docs: def interpolate_gridded_scalar(self, x, y, c, order=1, pad=1, offset=0):
          """Interpolate gridded scalar C to points x,y.
          
          Parameters
          ----------
          x, y : array-like
              Points at which to interpolate
          c : array-like
              The scalar, assumed to be defined on the grid.
          order : int
              Order of interpolation
          pad : int
              Number of pad cells added
          offset : int
              ??? 
              
          Returns
          -------
          ci : array-like
              The interpolated scalar
          """
          
          ## no longer necessary because we accept pre-padded arrays
          # assert c.shape == (self.Ny, self.Nx), 'Shape of c needs to be (Ny,Nx)'
          
          # first pad the array to deal with the boundaries
          # (map_coordinates can't seem to deal with this by itself)
          # pad twice so cubic interpolation can be used
          if pad > 0:
              cp = self._pad_field(c, pad=pad)
          else:
              cp = c
          # now the shape is Nx+2, Nx+2
          i = (x - self.xmin)/self.Lx*self.Nx + pad + offset - 0.5
          j = (y - self.ymin)/self.Ly*self.Ny + pad + offset - 0.5
          
          # for some reason this still does not work with high precision near the boundaries
          return scipy.ndimage.map_coordinates(cp, [j,i],
                  mode='constant', order=order, cval=np.nan)