docs: def interpolate_data(self, data, times, resampled_times):
          """ Interpolates data feature
  
          :param data: Array in a shape of t x nobs, where nobs = h x w x n
          :type data: numpy.ndarray
          :param times: Array of reference times relative to the first timestamp
          :type times: numpy.array
          :param resampled_times: Array of reference times relative to the first timestamp in initial timestamp array.
          :type resampled_times: numpy.array
          :return: Array of interpolated values
          :rtype: numpy.ndarray
          """
          # pylint: disable=too-many-locals
          # get size of 2d array t x nobs
          nobs = data.shape[-1]
          if self.interpolate_pixel_wise:
              # initialise array of interpolated values
              new_data = data if self.resample_range is None else np.full(
                  (len(resampled_times), nobs),
                  np.nan, dtype=data.dtype)
  
              # Interpolate for each pixel, could be easily parallelized
              for obs in range(nobs):
                  valid = ~np.isnan(data[:, obs])
  
                  obs_interpolating_func = self.get_interpolation_function(times[valid], data[valid, obs])
  
                  new_data[:, obs] = obs_interpolating_func(resampled_times[:, np.newaxis])
  
              # return interpolated values
              return new_data
  
          # mask representing overlap between reference and resampled times
          time_mask = (resampled_times >= np.min(times)) & (resampled_times <= np.max(times))
  
          # define time values as linear monotonically increasing over the observations
          const = int(self.filling_factor * (np.max(times) - np.min(times)))
          temp_values = (times[:, np.newaxis] +
                         const * np.arange(nobs)[np.newaxis, :].astype(np.float64))
          res_temp_values = (resampled_times[:, np.newaxis] +
                             const * np.arange(nobs)[np.newaxis, :].astype(np.float64))
  
          # initialise array of interpolated values
          new_data = np.full((len(resampled_times), nobs), np.nan, dtype=data.dtype)
  
          # array defining index correspondence between reference times and resampled times
          ori2res = np.array([np.abs(resampled_times - o).argmin()
                              if np.min(resampled_times) <= o <= np.max(resampled_times) else None for o in times])
  
          # find NaNs that start or end a time-series
          row_nans, col_nans = np.where(self._get_start_end_nans(data))
          nan_row_res_indices = np.array([index for index in ori2res[row_nans] if index is not None], dtype=np.int32)
          nan_col_res_indices = np.array([index is not None for index in ori2res[row_nans]],
                                         dtype=np.bool)
          if nan_row_res_indices.size:
              # mask out from output values the starting/ending NaNs
              res_temp_values[nan_row_res_indices, col_nans[nan_col_res_indices]] = np.nan
          # if temporal values outside the reference dates are required (extrapolation) masked them to NaN
          res_temp_values[~time_mask, :] = np.nan
  
          # build 1d array for interpolation. Spline functions require monotonically increasing values of x,
          # so .T is used
          input_x = temp_values.T[~np.isnan(data).T]
          input_y = data.T[~np.isnan(data).T]
  
          # build interpolation function
          if len(input_x) > 1:
              interp_func = self.get_interpolation_function(input_x, input_y)
  
              # interpolate non-NaN values in resampled time values
              new_data[~np.isnan(res_temp_values)] = interp_func(res_temp_values[~np.isnan(res_temp_values)])
  
          # return interpolated values
          return new_data