@shoyer and @crusaderky That's right, that is how I was actually dealing with this problem prior trying xarray ... by flattening the grid coordinates and performing either gridding (with scipy's griddata) or interpolation (with scipy's map_coordinate) ... instead of performing proper regridding (from cube to cube without having to flatten anything).

As a rule of thumb, any fancy algorithm should first exist for numpy-only data and then potentially it can be wrapped by the xarray library.

This is important information.

For the record, here is so far what I found to be best performant:

``` import xarray as xr from scipy.interpolate import griddata

Here x/y are dummy 1D coords that wont be used.

da1 = xr.DataArray(cube1, [('t', t_cube1) , ('y', range(cube1.shape[1])), ('x', range(cube1.shape[2]))])

Regrid t_cube1 onto t_cube2 first since time will always map 1 to 1 between cubes.

This operation is very fast.

print('regridding in time ...') cube1 = da1.interp(t=t_cube2).values

Regrid each 2D field (X_cube1/Y_cube1 onto X_cube2/Y_cube2), one at a time

print('regridding in space ...') cube3 = np.full_like(cube2, np.nan) for k in range(t_cube2.shape[0]): print('regridding:', k) cube3[:,:,k] = griddata((X_cube1.ravel(), Y_cube1.ravel()), cube1[k,:,:].ravel(), (X_cube2, Y_cube2), fill_value=np.nan, method='linear') ```

Performance is not that bad... for ~150 time steps and ~1500 nodes in x and y it takes less than 10 min.

I think this can be sped up by computing the interpolation weights between grids in the first iteration and cache them (I think xESMF does this).

@fspaolo where does that huge number come from? I thought you said you have 1500 nodes in total.

Not in total, I meant ~1500 on each x/y dimension (1500 x 1500). To be precise:

``` print(cube1.shape) # (t, y, x) (306, 240, 262)

print(cube2.shape) # (y, x, t) (1452, 1836, 104) ```

Did you select a single point on the t dimension before you applied bisplrep?

Yes.

Also, (pardon the ignorance, I never dealt with geographical data) what kind of information does having your lat and lon being bidimensional convey? Does it imply lat[i, j] < lat[i +1, j] and lon[i, j] < lon[i, j+1] for any possible (i, j)?

It does in my case, but cube1 is a rotated grid, meaning that lat is not the same along the x-axis neither lon is the same along the y-axis, while cube2 is a standard lon/lat grid so I can represent it simply by a 1D lon array (x-axis) and a 1D lat array (y-axis). To have them both with 1D coordinates I would have to either rotate cube2 and work with 1D rotated lon/lat, or unrotate cube1 and work with 1D standard lon/lat... but this gets me to the same problem as in the interpolation above since I have to transform (re-project) every 2D grid in the cube.

@crusaderky I also did the above using bisplrep and bispev, and it seems it cannot handle a grid of this size:

``` print(len(x), len(y), len(z)) 62880 62880 62880

print(len(x_new), len(y_new)) 2665872 2665872

Traceback (most recent call last): File "cubesmb.py", line 201, in <module> z_new = bisplev(x_new, y_new, tck) File "/Users/paolofer/anaconda3/lib/python3.7/site-packages/scipy/interpolate/_fitpack_impl.py", line 1047, in bisplev z, ier = _fitpack._bispev(tx, ty, c, kx, ky, x, y, dx, dy) RuntimeError: Cannot produce output of size 2665872x2665872 (size too large) ```

So griddata looks like my best option...

@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?

@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 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...

@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!

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).