Regridding with 3D Mask

[christophrenkl]

Hello, I have raised this question in the pangeo forum and was recommended to address it here:

I am currently trying to create initial conditions for a high-resolution, regional ocean model (based on ROMS) from the GLORYS12v1 reanalysis dataset. Since the model grids and vertical coordinate systems are different I am interpolating the three-dimensional (longitude, latitude, depth) GLORYS data horizontally onto the model grid using the xESMF module which generally works like a charm.

To account for the differences in the land/sea masks between my model and the reanalysis, I am using the extrapolation feature. This requires a mask variable which can only be 2D because of the underlying ESMF design.

However, the land/sea mask in GLORYS12v1 changes with depth, so that I have to apply the regridding for each vertical level separately. Is there any preferred way to do that?

Suppose that ds is a xarray.Dataset() with dimensions depth, lon, and lat containing multiple variables including a 3D mask. Off the top of my hat, I can think of two options:

1. Iterate over depth dimension, build the xesmf.Regridder() for each level and manually create an output dataset ds_int with the interpolated data. I’m afraid that this could be quite slow.
2. Write a dummy function that takes 2D fields as input arguments and applies the regridding for an individual model level. I would then apply this function to all depth levels using the xarray.apply_ufunc() method.

Do you have any other suggestions, tips, or tricks how to best approach this (in the best case computationally efficiently)?

Many thanks in advance!

[raphaeldussin]

Hi @christophrenkl this is a feature that, as an oceanographer, would like very much to have as well. AKAIK the ESMF library does not have the 3D interpolation available so xESMF cannot produce this natively. However, I'm thinking that there could be a solution:

When going from the surface to the deep ocean, we are only going to remove ocean points in the source and target grids. Starting from the surface regridder, we can modify the regridding matrix for the successive levels to reflect those changes. This would only necessitate to compute it once and then only modify the points that changed from ocean to land. Such an algo could look like:

1. creating the regridder for the surface level
2. copy the regridding sparse matrix from the regridder
3. loop: edit it for the various levels by taking out new target land points and zeroing out new source land points and renormalizing the weights
4. loop: Apply updated sparse matrix to new depth level

I don't have the bandwith to take this on right now but If you want to give it a try, I'd be happy to help.

[trondkr]

Hi @christophrenkl,

Just to chime in with some additional info here. I have a python toolbox called model2roms (https://github.com/trondkr/model2roms) that does a lot of what you are asking. It creates initial, clim and bry files for running ROMS using ESMF for interpolation. Model2roms uses GLORYS12V1 as the default dataset to create the model forcing. The toolbox uses ESMF which could be switched to xESMF instead to make it somewhat easier to read.

However, I do not take into account depth-varying mask when doing the interpolation so this is something that could be added. The interpolation is done first horizontally using ESMF and then vertically (linear interpolation implemented in fortran) to get to sigma coordinates. Given time it would be great to convert the interpolation routines  to only use xESMF, with masking, and output directly to sigma coordinates.

If you do implement depth-varying masking it would be great if you would send a pull request so I could integrate it into model2roms. The toolbox could definitely benefit from the help of others :)

Cheers, Trond On Jun 23, 2021, 12:57 PM -0700, raphael dussin @.***>, wrote:

Hi @christophrenkl this is a feature that, as an oceanographer, would like very much to have as well. AKAIK the ESMF library does not have the 3D interpolation available so xESMF cannot produce this natively. However, I'm thinking that there could be a solution: When going from the surface to the deep ocean, we are only going to remove ocean points in the source and target grids. Starting from the surface regridder, we can modify the regridding matrix for the successive levels to reflect those changes. This would only necessitate to compute it once and then only modify the points that changed from ocean to land. Such an algo could look like:

1. creating the regridder for the surface level
2. copy the regridding sparse matrix from the regridder
3. loop: edit it for the various levels by taking out new target land points and zeroing out new source land points and renormalizing the weights
4. loop: Apply updated sparse matrix to new depth level

I don't have the bandwith to take this on right now but If you want to give it a try, I'd be happy to help. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[stefraynaud]

Hi @christophrenkl

PR #65 is proposed to handle scattered missing values, where 2D are just a particular case. It basically allows the extrapolation to not be contaminated by NaNs, and the output mask depends on the level of contamination. However, it does not extrapolate from the scattered mask, which can be solved by applying a kind of spatial filter to the coastal area.

[christophrenkl]

Thank you for your input so far! I will try to carve out some time to have a more detailed look at the proposed approaches.