Implement horizontal regridding of irregular grids

[mattiarighi]

The regridding of irregular grids to regular grids is availalbe in the regrid module but is not yet tested. As a test case, the sic variable (sea-ice concentration) will be used.

[ledm]

@valeriupredoi, I'd like to volunteer to test changes related to irregular grids. Let me know when you're ready for me to have a crack at it!

[valeriupredoi]

hey guys, I've just started looking into this just today -- I grabbed MPI-ESM-LR/MR files for sic variable and had a first looksee through the files. Data has shape = (time, len(i), len(j)) where (i, j) represent the indices along (dim1, dim2), in this case (lat, lon). The metadata comes with a set of (lat, lon) as well, as AuxCoords but len(j) differs from len(lon) (whereas len(i) = len(lat = 220). My question to the oceans and gridding experts @bjlittle @ledm and @mattiarighi is how to project the irregular grid onto a regular grid in the most precise way, without massive loss of generality. I found this that explains the relation between x,y (lat, lon) and p,q (i,j): http://polar.ncep.noaa.gov/waves/workshop/pdfs/wwws_2013_irregular_grids.pdf but the problem is that the indexing of (lat, lon) won't work if len(j) = len(lat). Any advice how to get a metric g=(i, j, lat, lon) to map to a regular (lat_rg, lon_rg) grid?

[tillku]

@valeriupredoi @ledm Generally, the question of how to correctly interpolate data from the irregular grids of some ocean models to regular (lat-lon) grids has been hotly debated in the ocean modelling community for a long time. Some people would like to do an interpolation (e.g. us a for a meaningful intercomparison of several CMIPi models), and other people object to it because any interpolation always means a loss of accuracy, and hence loss of conservation in (interpolated) tracer fields. But I guess you know this already. Back to the actual question. For the case of the ORCA configurations, at least the direct model output netcdf files and the mesh_grid netcdf file contain all the information about where exactly in lat-lon space every single grid vertex sits. From that, interpolation scripts (e.g. SCRIP) can calculate mapping weights for the actual interpolation - your metric g(i,j,lat,lon). I guess the problem you've got with the MPI data is that this mapping doesn't work because the information about the grid in the metadata doesn't match the grid that the actual data are on. Is that right? If so, I'd contact MPI to be honest. Perhaps it's a bug they've sorted out already.

[tillku]

one more thought @valeriupredoi : if you len(j) and len(lon) differ by only 1 or 2, it could be a wrap-around thing? in that case, if we assume len(j) is that bit larger, the first and last values are identical perhaps?

[koldunovn]

We have just discussed with @valeriupredoi different ways of interpolation. I will just live here several links that can be useful in this respect: Interpolation between grids with Basemap Interpolation between grids with cKDTree Interpolation between grids with pyresample

Several libraries that can do interpolation: Relativelly simple: scipy interp2d

More powerful, aware of geographcal coordinates: pyresample

Very powerful, but the API is not the most "pythonic" one: ESMpy

[bouweandela]

It would be best if the regridding could be done through iris, possibly using one of those packages. @bjlittle Can you say something about this?

[valeriupredoi]

hey guys, cheers muchly for the input -- great responses! I am thinking a simple interpolation will do for now, and then we can get to more complex solutions starting from this. @bouweandela iris has an interpolator: https://scitools.org.uk/iris/docs/v2.0/userguide/interpolation_and_regridding.html but am not sure how well it works for a case like ours, I have never used it myself. @bjlittle is the main man here :grin:

[koldunovn]

Sorry for double posting to some of you. Here is what is available now in iris: https://scitools.org.uk/iris/docs/latest/iris/iris/experimental/regrid.html

The functions are experimental for now.

I would still go for solutions like pyresample or especially ESMPy since they are more general. For them, it does not matter if the original data are on regular or the curvilinear grid (or unstructured mesh for that matter). So you write a function that covers all cases at once.

[mattiarighi]

Outcome of the workshop discsussion on this topic:

• the first goal would be to have a regridding package for dealing with 2D variables on irregular grids (like sic)
• Iris does not provide such a method (yet)
• the ESMF package is the preferred alternative and it is well established
• @koldunovn @zklaus and @valeriupredoi will work towards a first implementation within the ESMValTool .preprocessor._regrid.py module
• the order of operations of the preprocessor will not be changed for now

[zklaus]

This is implemented in the branch version2_irregular_2d_regridding. If someone wants to give it a spin, they are welcome, but at the moment automated tests are missing.

[mattiarighi]

Thanks @zklaus, this is great!

Could you open a pull request? This would allow us to see code changes and comment.

[mattiarighi]

Some specific model / variable combination may require a fix for the irregular regridding to work properly. This will be addressed in separate issues.

Thanks @zklaus for this essential new feature!

[valeriupredoi]

yes - big thanks from me as well @zklaus :beer:

[ledm]

Fantastic work @zklaus, I'm very glad that this is done. We can now use the main development branch of ESMValTool for evaluating ocean models. This is tremendous progress.

[zklaus]

Thanks guys! I'm happy too :)