3

I have two data arrays that are on lat-lon grids. The first one, A, has the shape (89, 180). The second one, B, has the shape (94, 192). A's latitudes are in descending order from 88. to -88. & longitudes are in ascending order from 0. to 358. B's latitudes are in descending order from 88.54199982 to -88.54199982 & longitudes are in ascending order from 0. to 358.125.

I want to regrid/interpolate B's data onto A's coordinate system so that I can get both arrays the same size and calculate the spatial correlation between them. (I can also regrid/interpolate A's data onto B's coordinate system if that's easier.) I tried mpl_toolkits.basemap.interp(datain, xin, yin, xout, yout), but that requires xout & yout to be the same size. I also tried scipy.interpolate.griddata, but I can't figure out how it works and I'm not even sure that'll get me what I'm looking for...

Cebbie
  • 1,741
  • 6
  • 23
  • 37
  • how bout simple linear interpolation? Interpolate the x coordinates and the y separately and it should work. Basically scale any coordinate by the grid difference. e.g., x_new = x_old*(B_width/A_width). e.g., if B_width was 2*A_width then x_new = x_old*2. Obviously if B_width = A_width then there is no scaling because they are the exact same size. – AbstractDissonance Mar 01 '16 at 21:44
  • Now, once you have the intermediate value(since it may not be integral) you'll need to create a "blend" of the data points around that value to get the new data value. In fact, it works best to work backwards. E.g., for each point (i,j) in B, the data value is B[i,j]. To get this value, we need to look it up in A. But we must "scale" our i,j into A which won't be integral points so we can't use them directly to access our array. We can either simply round/floor/ceil or interpolate between the values around A using a weighted sum or other methods(cubic interpolation). – AbstractDissonance Mar 01 '16 at 21:47
  • Sorry, I don't really understand your second comment... any way you could answer in code form? – Cebbie Mar 01 '16 at 21:52
  • e.g., B[i,j] = A[floor(i*A_width/B_width), floor(j*A_height/B_height)] is a simple mapping. But it is inaccurate since some of the points can fall between "grid points" in A yet we always "clamp" the points to the grid points(e.g., if i = 3.9 we use i = 3 which may not represent the value well. Rounding would be best but linear interpolation between 3(*0.1) and 4(*0.9) would be best. Hence a weighted approach is better depending on your needs. – AbstractDissonance Mar 01 '16 at 21:54

2 Answers2

2

You may want to look at pyresample for this and other similar geographical interpolation problems. It provides multiple methods for interpolation, works well with lat/lon data, and incorporates basemap support. I suggest this package because you can also create AreaDefinition objects that define a domain using Proj4 definitions, then register data to the AreaDefinition.

For your specific problem, I would do the following (note, the interpolation step is not complete, see below):

from pyresample.geometry import SwathDefinition
from pyresample.kd_tree import resample_nearest

def interp_b_to_a(a, b):
    '''Take in two dictionaries of arrays and interpolate the second to the first.
    The dictionaries must contain the following keys: "data", "lats", "lons"
    whose values must be numpy arrays.
    '''
    def_a = SwathDefinition(lons=a['lons'], lats=a['lats'])
    def_b = SwathDefinition(lons=b['lons'], lats=b['lats'])

    interp_dat = resample_nearest(def_b, b['data'], def_a, ...)
    new_b = {'data':interp_dat,
             'lats':copy(a['lats']),
             'lons':copy(a['lons'])
            }
    return new_b

Note that the interpolation step where resample_nearest is called is not complete. You will also need to specify the radius_of_influence which is the search radius to use around each point in meters. This is dependent on the resolution of your data. You may also want to specify nprocs to speed things up and fill_value if you are using masked data.

Vorticity
  • 4,582
  • 4
  • 32
  • 49
  • Are def_a & def_b meant to have the same definition? – Cebbie Mar 01 '16 at 22:00
  • No, def_a would be the lats and lots you want to interpolate to. def_b would be the lats and lots you want to interpolate from. Editing the answer now. – Vorticity Mar 01 '16 at 22:04
  • 1
    I began to test it but didn't know what to add for the ... section. In the end I didn't use any of the solutions on this page because I realized I didn't need to interpolate lat-lon grids to each other for the problem I was trying to solve. Thanks for your answer though. – Cebbie Mar 04 '16 at 18:52
  • Just FYI, the ... indicated other keywords that are case specific as was discussed in the last paragraph. So, if you need this in the future, just look at the documentation for pyresample and look at the additional keywords. – Vorticity Mar 07 '16 at 17:48
0

Based on @Vorticity, I solved it editing to:

from pyresample.geometry import SwathDefinition
from pyresample.kd_tree import resample_nearest

def interp_b_to_a(a, b):
    def_a = SwathDefinition(lons=a['lons'], lats=a['lats'])
    def_b = SwathDefinition(lons=b['lons'], lats=b['lats'])

    interp_dat = resample_nearest(def_a, a['data'], def_b,radius_of_influence = 5000)
    new_b = {'data':interp_dat,
             'lats':b['lats'],
             'lons':b['lons']
            }
    return new_b
Javier A. Concha
  • 1
  • 1