def plot(self, ax, V=None, **kwargs):
        '''Plot the contours into matplotlib axis.

        Parameters
        ----------
        ax : matplotlib.Axes
            Axes to plot into
        V : array-like
            A list of contour values to plot. If not None, the internal contour
            values will be overriden during plotting, but not inside the
            object.
        kwargs : dict
            Keyword arguments to pass on to the ax.contour() method.
        '''
        if V is None:
            V = self.V
        d, X, Y = self.data.getData()
        # hack - add zero value to close contours
        d = np.hstack((d, np.zeros((d.shape[0], 1))))
        d = np.vstack((d, np.zeros((1, d.shape[1]))))
        dx = X[0, 1] - X[0, 0]
        dy = Y[1, 0] - Y[0, 0]
        x_longer = X[0, :].tolist()
        x_longer.append(X[0, -1] + dx)
        y_longer = Y[:, 0].tolist()
        y_longer.append(Y[-1, 0] + dy)
        x_interp, y_interp = np.meshgrid(
            np.linspace(x_longer[0], x_longer[-1],
                        len(x_longer) * self.upsample_factor),
            np.linspace(x_longer[0], y_longer[-1],
                        len(y_longer) * self.upsample_factor))
        spl = RectBivariateSpline(x_longer, y_longer, d.T)
        d_interp = spl.ev(x_interp, y_interp)
        ax.contour(x_interp, y_interp, d_interp, V, **kwargs)

def resample2d(i_data, i_s, i_e, i_i, o_s, o_e, o_i, kx=3, ky=3, s=0, 
               gauss_sig=0, median_boxcar_size=0, clip=True):
  '''
    Resample a square 2D input grid with extents defined by [i_s] and [i_e] with 
    increment [i_i] to a new 2D grid with extents defined by [o_s] and [o_e] 
    with increment [o_i].
    
    Returns a 2D resampled array, with options for smoothing (gaussian and 
    median) and clipping.
  '''
  
  # calculate bivariate spline, G, using input grid and data
  grid_pre_rebin = np.arange(i_s, i_e, i_i)
  G = RectBivariateSpline(grid_pre_rebin, grid_pre_rebin, i_data, kx=kx, ky=ky)

  # evaluate this spline at new points on output grid
  grid_x, grid_y = np.mgrid[o_s:o_e:o_i, o_s:o_e:o_i]

  data = G.ev(grid_x, grid_y)
  
  if gauss_sig != 0:
    data = gaussian_filter(data, gauss_sig)
    
  if median_boxcar_size != 0:
    data = median_filter(data, median_boxcar_size)
    
  if clip:
    input_max = np.max(i_data)
    input_min = np.min(i_data)
    
    data[np.where(data>input_max)] = input_max
    data[np.where(data<input_min)] = input_min

  return data

  def setModel(self,model,dx,think_positive=False):
    '''
    Set new model for the source.

    :param model: ``(n, n)``
      Numpy image array.
    :param dx: scalar
      Pixel size in microarcseconds.
    :param think_positive: (optional) bool
        Should we enforce that the source image has no negative pixel values?
    '''
    self.nx = int(ceil(model.shape[-1] * dx / self.dx))          # number of image pixels
    self.model = model                                           # source model
    self.model_dx = dx                                           # source model resolution

    # load source image that has size and resolution compatible with the screen.
    self.isrc = np.empty(2*(self.nx,))
    self.think_positive = think_positive

    M = self.model.shape[1]       # size of original image array
    f_img = RectBivariateSpline(self.model_dx/self.dx*(np.arange(M) - 0.5*(M-1)),\
                                self.model_dx/self.dx*(np.arange(M) - 0.5*(M-1)),\
                                self.model)

    xx_,yy_ = np.meshgrid((np.arange(self.nx) - 0.5*(self.nx-1)),\
                          (np.arange(self.nx) - 0.5*(self.nx-1)),indexing='xy')

      
    m  = f_img.ev(yy_.flatten(),xx_.flatten()).reshape(2*(self.nx,))
    self.isrc  = m * (self.dx/self.model_dx)**2     # rescale for change in pixel size

    if self.think_positive:
      self.isrc[self.isrc < 0] = 0

    if not self.live_dangerously: self._checkSanity()

def calcAnker(IS, inputPoints, rasterdata, gp):
    """
    """
    dhm = rasterdata['subraster']
    [Xa, Ya, Xe, Ye] = inputPoints
    # Letzte Koordinate in xi/yi entspricht nicht exakt den Endkoordinaten
    Xe_ = gp['xi'][-1]
    Ye_ = gp['yi'][-1]

    AnkA_dist = IS['d_Anker_A'][0]
    AnkE_dist = IS['d_Anker_E'][0]
    stueA_H = IS['HM_Anfang'][0]
    stueE_H = IS['HM_Ende_max'][0]

    # X- und Y-Koordinate der Geodaten im Projektionssystem berechnen
    dx = float(Xe - Xa)
    dy = float(Ye - Ya)
    if dx == 0:
        dx = 0.0001
    azimut = math.atan(dy/dx)
    if dx > 0:
        azimut += 2 * math.pi
    else:
        azimut += math.pi
    # X- und Y-Koordinaten der beiden Ankerpunkte am Boden
    AnkXa = Xa - AnkA_dist * math.cos(azimut)
    AnkYa = Ya - AnkA_dist * math.sin(azimut)
    AnkXe = Xe_ + AnkE_dist * math.cos(azimut)
    AnkYe = Ye_ + AnkE_dist * math.sin(azimut)

    # Linear Interpolation
    # Koordinatenarrays des DHMs
    coordX = gp['linspaces'][0]
    coordY = gp['linspaces'][1]
    # kx, ky bezeichnen grad der interpolation, 1=linear
    spline = RectBivariateSpline(-coordY, coordX, dhm, kx=1, ky=1)
    xi = np.array([AnkXa, Xa, Xe_, AnkXe])
    yi = np.array([AnkYa, Ya, Ye_, AnkYe])
    # Z-Koordinate der Anker für Anfangs- und Endpunkte
    zAnker = spline.ev(-yi, xi)     # Höhenangaben am Boden

    AnkA_z = stueA_H + 0.1*(zAnker[1] - zAnker[0])
    AnkE_z = stueE_H + 0.1*(zAnker[2] - zAnker[3])

    if AnkA_dist == 0:
        AnkA_z = 0.0
    if AnkE_dist == 0:
        AnkE_z = 0.0

    Ank = [AnkA_dist, AnkA_z, AnkE_dist, AnkE_z]

    # Ausdehnungen der Anker Felder, alles in [m]
    #Ank = [d_Anker_A, z_Anker_A * 0.1, d_Anker_E, z_Anker_E * 0.1]
    Laenge_Ankerseil = (AnkA_dist**2 + AnkA_z**2)**0.5 + \
                       (AnkE_dist**2 + AnkE_z**2)**0.5

    # Eventuell nicht nötig
    #IS['z_Anker_A'][0] = z_Anker_A
    #IS['z_Anker_E'][0] = z_Anker_E
    return [Ank, Laenge_Ankerseil, zAnker]

    def interpolate_individual(self, image):
        # unpacking
        ogridx, ogridy = self.ogrid
        ngridx, ngridy = self.ngrid

        f = RectBivariateSpline(ogridy, ogridx, image, kx=1, ky=1)
        return f.ev(ngridy.flatten(), ngridx.flatten()).reshape(ngridx.shape)

    def __init__(self, alpha, Re, cl, cd):
        """Setup CCAirfoil from raw airfoil data on a grid.
        Parameters
        ----------
        alpha : array_like (deg)
            angles of attack where airfoil data are defined
            (should be defined from -180 to +180 degrees)
        Re : array_like
            Reynolds numbers where airfoil data are defined
            (can be empty or of length one if not Reynolds number dependent)
        cl : array_like
            lift coefficient 2-D array with shape (alpha.size, Re.size)
            cl[i, j] is the lift coefficient at alpha[i] and Re[j]
        cd : array_like
            drag coefficient 2-D array with shape (alpha.size, Re.size)
            cd[i, j] is the drag coefficient at alpha[i] and Re[j]
        """

        alpha = np.radians(alpha)
        self.one_Re = False

        # special case if zero or one Reynolds number (need at least two for bivariate spline)
        if len(Re) < 2:
            Re = [1e1, 1e15]
            cl = np.c_[cl, cl]
            cd = np.c_[cd, cd]
            self.one_Re = True

        kx = min(len(alpha)-1, 3)
        ky = min(len(Re)-1, 3)

        # a small amount of smoothing is used to prevent spurious multiple solutions
        self.cl_spline = RectBivariateSpline(alpha, Re, cl, kx=kx, ky=ky, s=0.1)
        self.cd_spline = RectBivariateSpline(alpha, Re, cd, kx=kx, ky=ky, s=0.001)

def getStraightenWormInt(worm_img, skeleton, half_width = -1, cnt_widths  = np.zeros(0), width_resampling = 7, ang_smooth_win = 12, length_resampling = 49):
    '''
        Code to straighten the worm worms.
        worm_image - image containing the worm
        skeleton - worm skeleton
        half_width - half width of the worm, if it is -1 it would try to calculated from cnt_widths
        cnt_widths - contour widths used in case the half width is not given
        width_resampling - number of data points used in the intensity map along the worm width
        length_resampling - number of data points used in the intensity map along the worm length
        ang_smooth_win - window used to calculate the skeleton angles. 
            A small value will introduce noise, therefore obtaining bad perpendicular segments.
            A large value will over smooth the skeleton, therefore not capturing the correct shape.
        
    '''
    #if np.all(np.isnan(skeleton)):
    #    buff = np.empty((skeleton.shape[0], width_resampling))
    #    buff.fill(np.nan)
    #    return buff
    assert half_width>0 or cnt_widths.size>0
    assert not np.any(np.isnan(skeleton))
    
    if ang_smooth_win%2 == 1:
        ang_smooth_win += 1; 
    
    if skeleton.shape[0] != length_resampling:
        skeleton, _ = curvspace(np.ascontiguousarray(skeleton), length_resampling)
    
    skelX = skeleton[:,0];
    skelY = skeleton[:,1];
    
    assert np.max(skelX) < worm_img.shape[0]
    assert np.max(skelY) < worm_img.shape[1]
    assert np.min(skelY) >= 0
    assert np.min(skelY) >= 0
    
    #calculate smoothed angles
    skel_angles = angleSmoothed(skelX, skelY, ang_smooth_win)
    
    #%get the perpendicular angles to define line scans (orientation doesn't
    #%matter here so subtracting pi/2 should always work)
    perp_angles = skel_angles - np.pi/2;
    
    #%for each skeleton point get the coordinates for two line scans: one in the
    #%positive direction along perpAngles and one in the negative direction (use
    #%two that both start on skeleton so that the intensities are the same in
    #%the line scan)
    
    #resample the points along the worm width
    if half_width <= 0:
        half_width = (np.median(cnt_widths[10:-10])/2.) #add half a pixel to get part of the contour
    r_ind = np.linspace(-half_width, half_width, width_resampling)
    
    #create the grid of points to be interpolated (make use of numpy implicit broadcasting Nx1 + 1xM = NxM)
    grid_x = skelX + r_ind[:, np.newaxis]*np.cos(perp_angles);
    grid_y = skelY + r_ind[:, np.newaxis]*np.sin(perp_angles);
    
    
    f = RectBivariateSpline(np.arange(worm_img.shape[0]), np.arange(worm_img.shape[1]), worm_img)
    return f.ev(grid_y, grid_x) #return interpolated intensity map

def main():
    
    filenameEffArea='aeff_P7REP_ULTRACLEAN_V15_back.fits'
    directoryEffectiveArea='/Users/dspolyar/Documents/IRF/EffectiveArea/' 
    print pyfits.info( directoryEffectiveArea+filenameEffArea) 
    CTHETA_LO, CTHETA_HI, energyLow, energyHigh, EFFAREA = importEffectiveArea(directoryEffectiveArea+filenameEffArea)
    energylog, Ctheta=centeringDataAndConvertingToLog(energyHigh,energyLow,CTHETA_HI,CTHETA_LO)
    SplineEffectiveArea=RectBivariateSpline(Ctheta,energylog,EFFAREA)
    plotofEffectiveArea(SplineEffectiveArea,EFFAREA,energylog,Ctheta)
    print SplineEffectiveArea.ev(1.,5.)

def main():
    
    if len(sys.argv) != 2:
        print """
Usage: python img2spline.py [path_to_img] \n"""
        sys.exit(-1)
        
    else:
        path_to_img = sys.argv[1]
    
    
    img = Image.open(path_to_img).convert('L') # RGB -> [0..255]
    print "Image was opened and converted to grayscale."
    img.show()
    
    width, height = img.size
    
    data = np.array(list(img.getdata()), dtype=int)
    data = data.reshape((height, width))
    print "Data was extracted."
    
#    Start plotting original surface of image
    fig = plt.figure()
#    ax = fig.add_subplot(111, projection='3d')
    ax = fig.add_subplot(111)
    ax.invert_yaxis()
    x = range(0, width)
    y = range(0, height)
#    rev_y = range(height-1, -1, -1) # reverse y
    
    X, Y = np.meshgrid(x, y)
    print Y
    
    
    r_stride = 1 + width / 20 
    c_stride = 1 + height / 20
#    ax.plot_surface(X, Y, data, rstride=r_stride, cstride=c_stride)
    mappable = ax.pcolor(X, Y, data)
    plt.colorbar(mappable)
    ax.set_title("Original grayscale image")
    ax.set_xlabel('Width (px)')
    ax.set_ylabel('Height (px)')
    plt.draw()
#    Finish plotting original surface of image
    
    
#    2D Interpolation here
    spline = RectBivariateSpline(x, y, data)
    
    print spline.get_coeffs()
    
    
    
    
    plt.show()

  def scatter(self,move_pix=0,scale=1):
    '''
    Generate the scattered image which is stored in the ``iss`` member.

    :param move_pix: (optional) int 
      Number of pixels to roll the screen (for time evolution).
    :param scale: (optional) scalar
      Scale factor for gradient.  To simulate the scattering effect at another 
      wavelength this is (lambda_new/lambda_old)**2
    '''

    M = self.model.shape[-1]       # size of original image array
    N = self.nx                    # size of output image array

    #if not self.live_dangerously: self._checkSanity()

    # calculate phase gradient
    dphi_x,dphi_y = self._calculate_dphi(move_pix=move_pix)

    if scale != 1:
        dphi_x *= scale/sqrt(2.)
        dphi_y *= scale/sqrt(2.)

    xx_,yy = np.meshgrid((np.arange(N) - 0.5*(N-1)),\
                         (np.arange(N) - 0.5*(N-1)),indexing='xy')

    # check whether we care about PA of scattering kernel
    if self.pa != None:
      f_model = RectBivariateSpline(self.model_dx/self.dx*(np.arange(M) - 0.5*(M-1)),\
                                    self.model_dx/self.dx*(np.arange(M) - 0.5*(M-1)),\
                                    self.model)

      # apply rotation
      theta = -(90 * pi / 180) + np.radians(self.pa)     # rotate CW 90 deg, then CCW by PA
      xx_ += dphi_x
      yy  += dphi_y
      xx = cos(theta)*xx_ - sin(theta)*yy
      yy = sin(theta)*xx_ + cos(theta)*yy
      self.iss  = f_model.ev(yy.flatten(),xx.flatten()).reshape((self.nx,self.nx))

      # rotate back and clip for positive values for I
      if self.think_positive:
          self.iss  = clip(rotate(self.iss,-1*theta/np.pi*180,reshape=False),a_min=0,a_max=1e30) * (self.dx/self.model_dx)**2
      else:
          self.iss  = rotate(self.iss,-1*theta/np.pi*180,reshape=False) * (self.dx/self.model_dx)**2

    # otherwise do a faster lookup rather than the expensive interpolation.
    else:
      yyi = np.rint((yy+dphi_y+self.nx/2)).astype(np.int) % self.nx
      xxi = np.rint((xx_+dphi_x+self.nx/2)).astype(np.int) % self.nx
      if self.think_positive:
        self.iss = clip(self.isrc[yyi,xxi],a_min=0,a_max=1e30)
      else:
        self.iss = self.isrc[yyi,xxi]

 def __init__(self, x, y, z, kx=1, ky=1, xname=None, xunits=None,
              yname=None, yunits=None, zname=None, zunits=None):
     """Constructor.
     """
     if hasattr(z, '__call__'):
         _x, _y = numpy.meshgrid(y, x)
         z = z(_x, _y)
     xBivariateSplineBase.__init__(self, x, y, z, xname, xunits, yname,
                                   yunits, zname, zunits)
     RectBivariateSpline.__init__(self, x, y, z,
                                  bbox=[None, None, None, None],
                                  kx=kx, ky=ky, s=0)

 def x_sig( self, x, sigma ):
     eps_list, mu_q = self.spirrid_response
     eps_sig = InterpolatedUnivariateSpline( mu_q[0, :], eps_list[1] )
     if max( mu_q ) > sigma:
         pass
     else:
         raise ValueError( 'applied stress higher than the maximum in micromechanical evaluation of a CB' )
     eps = eps_sig( sigma )
     spline = RectBivariateSpline( eps_list[0], eps_list[1], mu_q )
     sigma_f = spline.ev( x, ones( len( x ) ) * eps ) / self.V_f
     sigma_m = ( sigma - sigma_f * self.V_f ) / self.V_m
     return sigma_m

def _approx(fmapnii, s=14.):
    """
    Slice-wise approximation of a smooth 2D bspline
    credits: http://scipython.com/book/chapter-8-scipy/examples/two-dimensional-interpolation-\
    with-scipyinterpolaterectbivariatespline/

    """
    from scipy.interpolate import RectBivariateSpline
    from builtins import str, bytes

    if isinstance(fmapnii, (str, bytes)):
        fmapnii = nb.load(fmapnii)

    if not isinstance(s, (tuple, list)):
        s = np.array([s] * 2)

    data = fmapnii.get_data()
    zooms = fmapnii.header.get_zooms()

    knot_decimate = np.floor(s / np.array(zooms)[:2]).astype(np.uint8)
    knot_space = np.array(zooms)[:2] * knot_decimate

    xmax = 0.5 * data.shape[0] * zooms[0]
    ymax = 0.5 * data.shape[1] * zooms[1]

    x = np.arange(-xmax, xmax, knot_space[0])
    y = np.arange(-ymax, ymax, knot_space[1])

    x2 = np.arange(-xmax, xmax, zooms[0])
    y2 = np.arange(-ymax, ymax, zooms[1])

    coeffs = []
    nslices = data.shape[-1]
    for k in range(nslices):
        data2d = data[..., k]
        data2dsubs = data2d[::knot_decimate[0], ::knot_decimate[1]]
        interp_spline = RectBivariateSpline(x, y, data2dsubs)

        data[..., k] = interp_spline(x2, y2)
        coeffs.append(interp_spline.get_coeffs().reshape(data2dsubs.shape))

    # Save smoothed data
    hdr = fmapnii.header.copy()
    caff = fmapnii.affine
    datanii = nb.Nifti1Image(data.astype(np.float32), caff, hdr)

    # Save bspline coeffs
    caff[0, 0] = knot_space[0]
    caff[1, 1] = knot_space[1]
    coeffnii = nb.Nifti1Image(np.stack(coeffs, axis=2), caff, hdr)

    return datanii, coeffnii

    def add_dem_2D(self, x, dem, y0=0., y1=np.infty, yref=None, kx=3, ky=1,
                   s=None):
        '''
        Add topography by vertically stretching the domain in the region [y0,
        y1] - points below y0 are kept fixed, points above y1 are moved as
        the DEM, points in between are interpolated.

        Usage: first call add_dem_2D for each boundary that is to be perturbed
            and finally call apply_dem to add the perturbation to the mesh
            coordinates.

        :param x: x coordinates of the DEM
        :type x: numpy array
        :param dem: the DEM
        :type dem: numpy array
        :param y0: vertical coordinate, at which the stretching begins
        :type y0: float
        :param y1: vertical coordinate, at which the stretching ends, can be
            infinity
        :type y1: float
        :param yref: vertical coordinate, at which the stretching ends
        :type yref: float
        :param kx: horizontal degree of the spline interpolation
        :type kx: integer
        :param ky: vertical degree of the spline interpolation
        :type ky: integer
        :param s: smoothing factor
        :type s: float
        '''

        if not self.ndim == 2:  # pragma: no cover
            raise ValueError('apply_dem_2D works on 2D meshes only')

        if yref is None:
            yref = self.points[:, 1].max()

        if y1 < np.infty:
            y = np.array([y0, yref, y1])
            d = np.c_[np.zeros(len(dem)), dem, np.zeros(len(dem))]
        else:
            y = np.array([y0, yref])
            d = np.c_[np.zeros(len(dem)), dem]

        xx, yy = np.meshgrid(x, y, indexing='ij')
        rbs = RectBivariateSpline(x, y, d, kx=kx, ky=ky, s=s)

        # add to topography
        if self.topography is None:
            self.topography = np.zeros_like(self.points[:, -1])

        self.points[:, 1] += rbs.ev(self.points[:, 0], self.points[:, 1])

def kde_histogram(events_x, events_y, xout=None, yout=None, bins=None):
    """ Histogram-based Kernel Density Estimation

    Parameters
    ----------
    events_x, events_y: 1D ndarray
        The input points for kernel density estimation. Input
        is flattened automatically.
    xout, yout: ndarray
        The coordinates at which the KDE should be computed.
        If set to none, input coordinates are used.
    bins: tuple (binsx, binsy)
        The number of bins to use for the histogram.

    Returns
    -------
    density: ndarray, same shape as `xout`
        The KDE for the points in (xout, yout)

    See Also
    --------
    `numpy.histogram2d`
    `scipy.interpolate.RectBivariateSpline`
    """
    valid_combi = ((xout is None and yout is None) or
                   (xout is not None and yout is not None)
                   )
    if not valid_combi:
        raise ValueError("Both `xout` and `yout` must be (un)set.")

    if yout is None and yout is None:
        xout = events_x
        yout = events_y

    if bins is None:
        bins = (max(5, bin_num_doane(events_x)),
                max(5, bin_num_doane(events_y)))

    # Compute the histogram
    hist2d, xedges, yedges = np.histogram2d(x=events_x,
                                            y=events_y,
                                            bins=bins,
                                            normed=True)
    xip = xedges[1:]-(xedges[1]-xedges[0])/2
    yip = yedges[1:]-(yedges[1]-yedges[0])/2

    estimator = RectBivariateSpline(x=xip, y=yip, z=hist2d)
    density = estimator.ev(xout, yout)
    density[density < 0] = 0

    return density.reshape(xout.shape)

def run_slimscat(isrc,idx,screenfile='screen.bin'):
  '''
  Scatter source image.

  :param isrc: source image
  :param idx: source pixel scale
  :param screenfile: screen file

  '''

  # read in screen parameters
  f = open(screenfile,'rb')
  hdrsize = struct.unpack('i',f.read(4))[0]
  nphi = int(np.sqrt(struct.unpack('i',f.read(4))[0]))
  dx   = struct.unpack('d',f.read(8))[0]
  f.seek(hdrsize)

  # filter image

  # check fov
  iny,inx = isrc.shape
  ny = int(np.floor(iny*idx/dx))
  nx = int(np.floor(inx*idx/dx))
  assert idx*max([iny,inx]) < nphi*dx

  # read in screen
  dphi_x = np.empty((ny,nx),dtype=np.float64)
  dphi_y = np.empty((ny,nx),dtype=np.float64)
  for i in range(ny):
    f.seek(hdrsize + i*8*1*nphi,0)
    dphi_x[i,:] = struct.unpack('{0:d}d'.format(nx),f.read(nx*8))
  for i in range(ny):
    f.seek(hdrsize + 8*nphi*nphi + i*8*1*nphi,0)
    dphi_y[i,:] = struct.unpack('{0:d}d'.format(nx),f.read(nx*8))

  f.close()

  # construct spline fit to source image
  f_isrc = RectBivariateSpline(idx/dx*(np.arange(iny) - 0.5*(iny-1)),\
                               idx/dx*(np.arange(inx) - 0.5*(inx-1)),\
                               isrc)

  # scatter pixel coordinates
  yy,_xx = np.meshgrid((np.arange(ny) - 0.5*(ny-1)),\
                       (np.arange(nx) - 0.5*(nx-1)),indexing='ij')
  _xx += dphi_x
  yy += dphi_y

  return f_isrc.ev(yy.flatten(),_xx.flatten()).reshape((ny,nx))

def interp_cross(x_in, y_in, data, x0, x1, y0, y1, rint=None, mode='linear'):
    u"""
    (x0,y0)-(x1,y1)に沿った直線上の点に内挿する
    """
    if np.ndim(x_in) == 2: x_in = x_in[0,:]
    if np.ndim(y_in) == 2: y_in = y_in[:,0]

    ndim = np.ndim(data)
    if ndim > 2:
        oldshape = data.shape        
        data     = data.reshape(-1, oldshape[-2], oldshape[-1])
    elif ndim < 2:
        raise ValueError, "input data dimension size must be larger than 2. input ndim is {}".format(ndim)
        
    # make output grid
    if rint is None:
        rint = max(np.diff(x_in).max(), np.diff(y_in).max())
    rmax  = np.hypot(x1-x0, y1-y0)
    nr    = rmax//rint + 1
    theta = np.arctan2(y1-y0, x1-x0)
    r_out = np.linspace(0, rmax, nr)
    x_out = r_out*np.cos(theta) + x0
    y_out = r_out*np.sin(theta) + y0

    xrev = yrev = 1
    if x0 > x1:
        x_out = x_out[::-1]
        xrev  = -1
    if y0 > y1:
        y_out = y_out[::-1]
        yrev  = -1
        
    # interpolate    
    if ndim == 2:
        if mode == 'linear':
            out = RectBivariateSpline(y_in, x_in, data, kx=1, ky=1)(y_out, x_out)[::xrev,::yrev].diagonal()
        elif mode == 'cubic':
            out = RectBivariateSpline(y_in, x_in, data, kx=3, ky=3)(y_out, x_out)[::xrev,::yrev].diagonal()
    else:
        zn, yn, xn = data.shape
        out = np.ma.empty((zn,nr), dtype=np.float32)
        for z in range(zn):
            if mode == 'linear':
                out[z,:] = RectBivariateSpline(y_in, x_in, data[z,:,:], kx=1, ky=1)(y_out, x_out)[::xrev,::yrev].diagonal()
            elif mode == 'cubic':
                out[z,:] = RectBivariateSpline(y_in, x_in, data[z,:,:], kx=3, ky=3)(y_out, x_out)[::xrev,::yrev].diagonal()
        out = out.reshape(oldshape[:-2] + (-1,))
                
    return x_out, y_out, r_out, out

def get_interpolated_pixel_color_rbspline(pts, s_im, size):
    """given pts in floats, linear interpolate pixel values nearby to get a good colour"""
    pts = clamp(pts, size)

    s_im = np.atleast_3d(s_im)
    ys,xs = size
    ycoords, xcoords = np.arange(ys), np.arange(xs)
    out = np.empty(pts.shape[1:] + (s_im.shape[-1],),dtype=s_im.dtype)
    
    pts_vec = pts.reshape((2,-1))
    out_vec = out.reshape((-1,s_im.shape[-1])) #flatten for easier vectorization
    for i in range(s_im.shape[-1]): #loop over color channels
        rbspline = RectBivariateSpline(ycoords, xcoords, s_im[...,i])
        out_vec[:,i] = rbspline.ev(pts_vec[0],pts_vec[1])
    return out

def interp3d_emergence(uplift, data, out_times, verbose=False):
    """Interpolate uplift surfaces (xyz data at a specific t) to data locations 
    (non-grid) and data times (between times calculated). 
    
    Uses progressive linear interpolations: first the uplift at each outputted 
    time is interpolated to the data locations in data.locs, then they are
    interpolated to the data times in each location.
    
    Parameters
    ----------
    uplift (array-like) - size (times, lon, lat) array of uplift surfaces
    data - data whose data.locs are the locations to interpolate to.
    out_times - the times for the first index of the uplift array (should be
            of uplift eventually, yes?)
    """
    time_start = time.clock()
    ##########################################
    # STUFF TO FIX HEEEEEEERE!!!!!!!!!!!!!  
    N = uplift[0].shape
    # TODO These should be gotten from somewhere, right? uplift.grid?? 
    X = np.linspace(0, 4900000, num=N[1], endpoint=True)
    Y = np.linspace(0, 4700000, num=N[0], endpoint=True)
    ##########################################
    
    # interp_data will be an array of size (N_output_times, N_locations)
    # for use in interpolating the calculated emergence to the locations and
    # times at which there are data in data
    interp_data = []
    # Interpolate the calculated uplift at each time on the Lat-Lon grid
    # to the data locations.
    for uplift_at_a_time in uplift:
        interp_func = RectBivariateSpline(X, Y, uplift_at_a_time.T)
        interp_data.append(interp_func.ev(data.locs[:,0], data.locs[:,1]))
    interp_data = np.array(interp_data).T
    
    calc_vector = []
    # Interpolate the calculated uplifted at each time and data location
    # to the times of the data location.
    for interp, loc in zip(interp_data, data):
        calc_vector.append(np.interp(loc['data_dict']['times'],
        out_times[::-1], interp[::-1]))
    
    # flatten the array    
    calc_vector = np.array([item for l in calc_vector for item in l])
    
    if verbose: print 'Interpolation time: {0}s'.format(time.clock()-time_start)

    return calc_vector

	def __init__(self, terrainAP=None, rasterDist=1.0, humusType=None):
		"""
		The humus layer has a triangular thickness distribution with parameters in dPar.
		RasterDist is the distance between each node on our raster with thicknesses given by the distribution.
		A large distance mimics a situation where the humuslayer is more or less spatially invariant,
		whilst a small distance mimics a situation where the thickness varies on a small length scale.
		"""
		if humusType=='2': dPar=[0.00,0.05,0.01]#triangular distribution [min,max,mode] meters
		elif humusType=='3': dPar=[0.05,0.15,0.10]
		elif humusType=='4': dPar=[0.15,0.30,0.22]
		else: raise Exception('HumusType is not correct.')
		xLim = terrainAP[1][0]
		yLim = terrainAP[3][1]
		self.xNumInv = int(xLim/rasterDist)
		self.yNumInv = int(yLim/rasterDist)
		z=[]
		self.pointDepth = 1
		x=np.linspace(0,xLim,self.xNumInv)
		y=np.linspace(0,yLim,self.yNumInv)
		#These forloops are in order to achieve a two dimensional array with thicknesses
		for i in range(self.xNumInv):
			ztemp=[]
			for j in range(self.yNumInv):
				ztemp.append(random.triangular(dPar[0],dPar[1],dPar[2]))
			z.append(ztemp)
		self.z=np.array(z)#Must be a numpy-array for the interpolation to work.
		self.interpolDepth=RectBivariateSpline(x,y,self.z)