[FIXED] Basemap interpolation alternative - regridding data

Issue

I'm moving from basemap to cartopy given basemap is going to be phased out. I've previously used the basemap.interp functionality to interpolate data, e.g. say I have data at 1 degree resolution (180x360), I would run the following to interpolate to 0.5 degrees.

import numpy as np
from mpl_toolkits import basemap

Old_Lon = np.linspace(-180,180,360)
Old_Lat = np.linspace(-90,90,180)
New_Lon = np.linspace(-180,180,720)
New_Lat = np.linspace(-90,90,360)

New_Lon,New_Lat = np.meshgrid(New_Lon,New_Lat)

New_Data = basemap.interp(Old_Data,Old_Lon,Old_Lat,New_Lon,New_Lat,order=0)

order gives me options to choose from nearest neighbour, bi-linear etc. Is there an alternative that does this in as simple way? I've seen scipy has interpolation but I'm not sure how to apply it. Any help would be appreciated!

Solution

I eventually decided to take the raw code from Basemap and make it into a standalone function - I'll be recommending it to the cartopy guys to implement it as its a useful feature. Posting here as could be useful to someone else:

def Interp(datain,xin,yin,xout,yout,interpolation='NearestNeighbour'):

    """
       Interpolates a 2D array onto a new grid (only works for linear grids), 
       with the Lat/Lon inputs of the old and new grid. Can perfom nearest
       neighbour interpolation or bilinear interpolation (of order 1)'

       This is an extract from the basemap module (truncated)
    """

    # Mesh Coordinates so that they are both 2D arrays
    xout,yout = np.meshgrid(xout,yout)

   # compute grid coordinates of output grid.
    delx = xin[1:]-xin[0:-1]
    dely = yin[1:]-yin[0:-1]

    xcoords = (len(xin)-1)*(xout-xin[0])/(xin[-1]-xin[0])
    ycoords = (len(yin)-1)*(yout-yin[0])/(yin[-1]-yin[0])


    xcoords = np.clip(xcoords,0,len(xin)-1)
    ycoords = np.clip(ycoords,0,len(yin)-1)

    # Interpolate to output grid using nearest neighbour
    if interpolation == 'NearestNeighbour':
        xcoordsi = np.around(xcoords).astype(np.int32)
        ycoordsi = np.around(ycoords).astype(np.int32)
        dataout = datain[ycoordsi,xcoordsi]

    # Interpolate to output grid using bilinear interpolation.
    elif interpolation == 'Bilinear':
        xi = xcoords.astype(np.int32)
        yi = ycoords.astype(np.int32)
        xip1 = xi+1
        yip1 = yi+1
        xip1 = np.clip(xip1,0,len(xin)-1)
        yip1 = np.clip(yip1,0,len(yin)-1)
        delx = xcoords-xi.astype(np.float32)
        dely = ycoords-yi.astype(np.float32)
        dataout = (1.-delx)*(1.-dely)*datain[yi,xi] + \
                  delx*dely*datain[yip1,xip1] + \
                  (1.-delx)*dely*datain[yip1,xi] + \
                  delx*(1.-dely)*datain[yi,xip1]

    return dataout

--

Answered By - J W