Does interp() work on curvilinear grids (2D coordinates) ?

[JiaweiZhuang]

I am evaluating interp() against xESMF. Here's how xESMF would regrid the built-in 'rasm' data to a regular lat-lon grid.

Seems like interp() can convert rectilinear grids (1D coordinates) to curvilinear grids (2D coordinates), according to the last example. How about the reverse? Can it convert 2D coordinate to 1D, or to another 2D coordinate?

That's the test data:

dr = xr.tutorial.load_dataset('rasm')['Tair']
...
Coordinates:
  * time     (time) datetime64[ns] 1980-09-16T12:00:00 1980-10-17 ...
    xc       (y, x) float64 189.2 189.4 189.6 189.7 189.9 190.1 190.2 190.4 ...
    yc       (y, x) float64 16.53 16.78 17.02 17.27 17.51 17.76 18.0 18.25 ...

That's a simple destination grid:

lon = xr.DataArray(np.linspace(-180, 180, 120), dims=['lon'])
lat = xr.DataArray(np.linspace(-90, 90, 60), dims=['lat'])

I would expect a syntax like:

dr_out = dr.interp(xc=lon, yc=lat)

But I got ValueError: dimensions ['xc', 'yc'] do not exist, becauseinterp() only accepts dimensions (which must be 1D), not coordinates.

dr.interp(x=lon, y=lat) runs but the result is not correct. This is expected because x does not mean longitude in the original data.

@crusaderky @fujiisoup

[JiaweiZhuang]

One way I can think of to make interp() work on this example: Define a new coordinate system (i.e. two new coordinate variables) on the source curvilinear grid, and rewrite the destination coordinate using those new coordinate variables (not lat, lon anymore).

But this is absolutely too convoluted...

Updated: see Gridded with Scipy for a similar idea.

[fujiisoup]

Thanks, @JiaweiZhuang

Not yet.
interp() only works on N-dimensional regular grid.
Under the hood, we are just using scipy.interpolate.interp1d and interpn.

I am happy to see curvilinear interpolation in xarray if we could find a good general API for N-dimensional array.
Do you have any proposal?

For curvilinear interpolation, we may have some arbitrariness,
e.g.

dr_out = dr.interp(xc=lon)

the resultant dimension is not well determined. Maybe we need some limitation for the arguments.

[shoyer]

I think we could make dr.interp(xc=lon, yc=lat) work for the N-D -> M-D case by wrapping scipy.interpolate.griddata

[JiaweiZhuang]

Do you have any proposal?

I guess it is not an API design problem yet... The algorithm is not here since interpn doesn't deal with curvilinear grids.

I think we could make dr.interp(xc=lon, yc=lat) work for the N-D -> M-D case by wrapping scipy.interpolate.griddata

My concern with scipy.interpolate.griddata is that the performance might be miserable... griddata takes an arbitrary stream of data points in a D-dimensional space. It doesn't know if those source data points have a gridded/mesh structure. A curvilinear grid mesh needs to be flatten into a stream of points before passed to griddata(). Might not be too bad for nearest-neighbour search, but very inefficient for linear/bilinear method, where knowing the mesh structure beforehand can save a lot of computation.

Utilizingscipy.interpolate.griddata would be a nice feature, but it should probably be used for data point streams (more like a Pandas dataframe method?), not as a way to handle curvilinear grids.

PS: I have some broader concerns regarding interp vs xESMF: JiaweiZhuang/xESMF#24

[shoyer]

I'd like to figure out interfaces that make it possible for external, grid aware libraries to extend indexing and interpolation features in xarray. In particular, it would be nice to be able to associate a "grid index" used for caching computation that gets passed on in all xarray operations.

[fspaolo]

Great thread by @JiaweiZhuang! I just posted a question on stackoverflow about this exact problem. After hours navigating through the xarray documentation I finally found this opened issue... (since I'm new to xarray I thought the problem was my lack of understanding on how the package works)

It is surprising that a package targeting n-dimensional gridded datasets (particularly those from the geo/climate sciences) does not handle such a common task with spatial gridded data.

The problem on hand is this: I have two 3d arrays with different dimensions defined by 2d coordinates, all I want is to regrid the first cube onto the second. Is there a way to perform this operation with xarray?

This is what I've tried (which @JiaweiZhuang explained why it doesn't work):

da = xr.DataArray(cube, dims=['t', 'y', 'x'],
                  coords={'t': time,           
                          'xc': (['y', 'x'], X),  
                          'yc': (['y', 'x'], Y)})

da_interp = da.interp(x=x_new, y=y_new).interp(t=t_new)

Here x_new/y_new should be mapped onto xc/yc (physical coordinates), not onto x/y (logical coordinates).

[crusaderky]

I am not aware of a ND mesh interpolation algorithm. However, my package xarray_extras [1] offers highly optimized 1D interpolation on a ND hypercube, on any numerical coord (not just time). You may try applying it 3 times on each dimension in sequence and see if you get what you want - although performance won't be optimal.

[1] https://xarray-extras.readthedocs.io/en/latest/

Alternatively, if you do find the exact algorithm you want, but it's for numpy, then applying it to xarray is simple - just get DataArray.values -> apply function -> create new DataArray from the output.

[fspaolo]

@crusaderky I don't think we need a "proper" 3d interpolation in most cases (i.e. predicting each 3d grid node considering all dimensions simultaneously). If you see my example above, DataArray.interp(x=x_new, y=y_new).interp(t=t_new), I am first interpolating over the spatial coordinates and then over time. If instead I do DataArray.interp(x=x_new, y=y_new, t=t_new), the computing time is prohibited for large ndarrays. In fact, I think DataArray.interp() should figure this out and do this decomposition internally when calling DataArray.interp(x=x_new, y=y_new, t=t_new).

The main limitation here, however, is being able to interpolate over the spatial coordinates when these are defined as 2d arrays. I'll check your package... thanks!

[crusaderky]

@fspaolo I never tried using my algorithm to perform 2D interpolation, but this should work:

from xarray_extras.interpolate import splrep, splev

da = splev(x_new, splrep(da, 'x'))
da = splev(y_new, splrep(da, 'y'))
da = splev(t_new, splrep(da, 't'))

Add k=1 to downgrade from cubic to linear interpolation and get a speed boost.

You can play around with dask to increase performance by using all your CPUs (or more with dask distributed), although you have to remember that an original dim can't be broken on multiple chunks when you apply splrep to it:

from xarray_extras.interpolate import splrep, splev
da = da.chunk(t=TCHUNK)
da = splev(x_new, splrep(da, 'x'))
da = splev(y_new, splrep(da, 'y'))
da = da.chunk(x=SCHUNK, y=SCHUNK).chunk(t=-1)
da = splev(t_new, splrep(da, 't'))
da = da.compute()

where TCHUNK and SCHUNK are integers you'll have to play with. The rule of thumb is that you want to have your chunks 5~100 MBs each.

If you end up finding out that chunking along an interpolation dimension is important for you, it is possible to implement it with dask ghosting techniques, just painfully complicated.

[fspaolo]

@crusaderky It doesn't work. First, it tells me that if x_new/y_new are 2D, then I have to pass them as DataArray's... why do I have to do that? After doing that, it tells me there are conflicts between "overlapping" dimensions?! Ok, I then pass x_new/y_new as 1D arrays... and the result comes nonsense!

The only interpolation that works (both with DataArray.interp() and your splev()) is the one over the time dimension. This is because time is defined as a one-dimensional variable (t).

Why is it so hard to perform an interpolation with spatial coordinates defined with 2D variables?! I would think this is a pretty common operation on climate datasets...

[fspaolo]

@JiaweiZhuang your approach doesn't work either. After installing you package and the dependencies... and following the documentation, I got

Fatal error in MPI_Init_thread: Other MPI error, error stack:
MPIR_Init_thread(474)..............:
MPID_Init(190).....................: channel initialization failed
MPIDI_CH3_Init(89).................:
MPID_nem_init(320).................:
MPID_nem_tcp_init(173).............:
MPID_nem_tcp_get_business_card(420):
MPID_nem_tcp_init(379).............: gethostbyname failed, LMC-061769 (errno 1)
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=3191311
:
system msg for write_line failure : Bad file descriptor

An MPI error?!

What bothers me are statements like this "For usability and simplicity"...

[crusaderky]

@fspaolo 2d mesh interpolation and 1d interpolation with extra "free" dimensions are fundamentally different algorithms. Look up the scipy documentation on the various interpolation functions available.

I don't understand what you are trying to pass for x_new and y_new and it definitely doesn't sound right. Right now you have a 3d DataArray with dimensions (x, y, t) and 3 coords, each of which is a 1d numpy array (e.g. da.coords.x.values). If you want to rescale, you need to pass a 1d numpy array or array-like for x_new, and another separate 1d array for y_new. You are not doing that, as the error message you're receiving is saying that your x_new is a numpy array with 2 or more dimensions, which the algorithm doesn't know what to do with. It can accept a multi-dimensional DataArrays with brand new dimensions, but that does not sound like it's your case.

[fspaolo]

@crusaderky what I'm trying to do is what the title of this opened thread says "Does interp() work on curvilinear grids (2D coordinates)?" Working with (x, y, t) all 1D numpy arrays is an easy problem! But that's not the problem I'm dealing with.

Let me state exactly what my problem is. I have two data cubes, cube1 and cube2:

cube1 = xr.DataArray(cube1, dims=['t', 'y', 'x'],
                     coords={'time': ('t', t_cube1),           # 1D array
                             'lon': (['y', 'x'], X_cube1),     # 2D array!!!
                             'lat': (['y', 'x'], Y_cube1)})    # 2D array!!!

cube2 = xr.DataArray(cube2, dims=['y', 'x', 't'],
                     coords={'time': ('t', t_cube2),           # 1D array
                             'lon': (['y', 'x'], X_cube2),     # 2D array!!!
                             'lat': (['y', 'x'], Y_cube2)})    # 2D array!!!

All I want is to regrid cube1 onto cube2, or in other words, interpolate X_cube2/Y_cube2/t_cube2 onto cube1. Note that I cannot pass X_cube1/Y_cube1 as 1D arrays because they vary for every grid cell (this is a rotated grid). I can, however, pass X_cube2/Y_cube2 as 1D arrays (this is a standard lon/lat grid).

Also note the inverted time dimension between cube1 and cube2 (but that shouldn't be an issue).

So how to perform this operation... or am I missing something?

[shoyer]

So how to perform this operation... or am I missing something?

Sorry, i don't think there's an easy way to do this directly in xarray right now.

My concern with scipy.interpolate.griddata is that the performance might be miserable... griddata takes an arbitrary stream of data points in a D-dimensional space. It doesn't know if those source data points have a gridded/mesh structure. A curvilinear grid mesh needs to be flatten into a stream of points before passed to griddata(). Might not be too bad for nearest-neighbour search, but very inefficient for linear/bilinear method, where knowing the mesh structure beforehand can save a lot of computation.

Thinking a little more about this, I wonder if this the performance could actually be OK as long as the spatial grid is not too big, i.e., if we reuse the same grid many times for different variables/times.

In particular, SciPy's griddata either makes use of a scipy.spatial.KDTree (for nearest neighbor lookups) and scipy.spatial.Delaunay(for linear interpolation, on a triangular mesh). We could build these data structures once (and potentially even cache them in indexes on xarray objects), and likewise calculate the sparse interpolation coefficients once for repeated use.