Extrapolation and masking

[RondeauG]

As discussed in Issue #93 (and earlier in #22), this adds the extrapolation functionalities of ESMPy 8.0.0 to xESMF. I also included masking functionalities that are compatible with ESMPy grids.

The changes are built on top of PR #91, since ESMP_FieldRegridStoreFile is outdated ("deprecated" even, based on a comment in the code) and does not include the new extrapolation functionalities. With #91, weights are kept in memory and the call to ESMF is instead made with ESMP_FieldRegridStore, which is up to date.

Masking is based on an automatic detection of a "mask" DataArray in the dataset.

[zxdawn]

Will this deal with #17? I'm curious of how to mask these diluted grids ...

[RondeauG]

@zxdawn I've not really explored conservative regridding yet. From what I see of #17, I think PR #82 is closer to what you're looking for, which is the addition of norm_type=ESMF.NormType.FRACAREA during the call to ESMPy.Regrid?

[RondeauG]

Once #82 is approved, I will replace my masking with his since I missed that a PR was already in the works for that. Both methods use a similar approach.

[trondkr]

@RondeauG I tried to apply your pull request 97 to the latest version of ESMF to experiment with extrapolation. The regridding to higher resolution works but I do not see that the extrapolation works and the extent of interpolation is limited to the extent of the original grid. Is there a trick to doing the extrapolation or some option I need to enable? My code is a copy of your test that looks like this regrid = xe.frontend.esmf_regrid_build(src_grid, dst_grid, 'bilinear', extrap='inverse_dist', extrap_num_pnts=3, extrap_exp=10). I would love to get the extrapolation to work :) Thanks!

[RondeauG]

Hi @trondkr . There are indeed a few tips and tricks scattered around the various issues threads. Your issue regarding the extent is very likely linked to #15. This is solved by using the add_matrix_NaNs function described in the thread.

# Temporary fix to make xESMF ignore grid cells outside of the domain
# https://github.com/JiaweiZhuang/xESMF/issues/15#issuecomment-609365496
def add_matrix_nans(r):
     """Deal with wrong boundary"""
     X = r.weights
     M = scipy.sparse.csr_matrix(X)
     num_nonzeros = np.diff(M.indptr)
     M[num_nonzeros == 0, 0] = np.NaN
     r.weights = scipy.sparse.coo_matrix(M)

     return r

As for your issue where the extrapolation does not seem to work, I have 2 possible fixes.

• First, your extrap_exp=10 seems high. Maybe that's the issue? I don't think this is it, but it should be fairly quick to test with extrap_exp=2.
• The other thing that comes to mind is that for the extrapolation to work, you really need to have proper masking to tell xESMF to ignore the pixels with NaNs. Otherwise, those NaNs will get interpolated/extrapolated. There's work in Masking with conservative_normed #82 to integrate masking in xESMF. Otherwise, my own PR Extrapolation and masking #97 uses a worse version. The main idea here is to feed xESMF a DataArray of 1's and 0's named "mask", with the 0's indicating pixels to ignore. In my own tests, for example, I used the land fraction sftlf to create the masks.

Here's the code I use locally to create my masks :

def _create_mask(ds, threshold):

    """ creates a SCRIP-compatible mask DataArray made of 0 (masked) and 1 (unmasked)

    ds : xr.Dataset()

    mask_in, mask_out: dict
      dict{variable: tuple(threshold, direction)}
      direction is either "smaller", "equal", "greater"
      NaNs will always be masked. If this is the only thing to mask, use np.nan as the threshold
      example: {"sftlf": (0, "equal")}

    :return mask : xr.DataArray()
    """

    v = list(threshold.keys())[0]
    t = threshold[v]

    if "time" in ds:
        ds = ds.isel(time=0)

    if t[1] in ["equal"]:
        mask = xr.where((xr.ufuncs.isnan(ds[v])) | (ds[v] == t[0]), 0, 1)
    elif t[1] in ["smaller"]:
        mask = xr.where((xr.ufuncs.isnan(ds[v])) | (ds[v] <= t[0]), 0, 1)
    elif t[1] in ["greater"]:
        mask = xr.where((xr.ufuncs.isnan(ds[v])) | (ds[v] >= t[0]), 0, 1)
    else:
        logging.error(t[1] + " not recognized.")
        raise ValueError

    return mask

And so, your xESMF code could look something like that :

# prepare the masks
src_grid["mask"] =  = _create_mask(src_sftlf, mask_conditions)
dst_grid["mask"] =  = _create_mask(dst_sftlf, mask_conditions)

# prepare the regridder
regridder = xe.Regridder(src_grid, dst_grid, 'bilinear', extrap='inverse_dist', extrap_num_pnts=3, extrap_exp=10)
regridder = add_matrix_nans(regridder)  # apply the grid extent workaround
da_out = regridder(src_grid[v], keep_attrs=True)

#(...) continue the code

[trondkr]

@RondeauG Thank you for the detailed help and explanation. I made a test script where I implemented your suggestions, but alas the extrapolation is still not doing what I hope it will do. I am interpolating and extrapolating (oxygen) from a 1/4 degree resolution grid to a 1/12th degree resolution grid (thetao). But the interpolation is only using the 1/4 degrees grid cells. I assume I am still doing something erroneous with regard to the masking. If you have time would you take a look at the short test script found on my Github here
The test script also includes two smaller src and dst netcdf files. Thank you very much!

[RondeauG]

@trondkr I believe that I have found the problem.

• _create_mask should be called before xe.Regridder, while in your script it comes after. (see my example above)
• Additionally, I did not account for variables other than time, lat, and lon in my _create_mask script. The mask variable should only have lat and lon as dimensions, but right now depth is still included. This line of code should be improved to include any variable that isn't lat or lon :

if "time" in ds:
        ds = ds.isel(time=0)

Alternatively, just to make sure that the extrapolation works (and it does, in my case), I just used src_ds["mask"] = _create_mask(src_ds, {"o2": (np.nan, "equal")}).isel(depth=0) and dst_ds["mask"] = _create_mask(dst_ds, {"o2": (np.nan, "equal")}).isel(depth=0)

Hope this helps!

[trondkr]

@RondeauG Thank you so much! That solved my problem and extrapolation works like a charm!

On a side note, I see that xesmf is moving over to pangeo-data but I do not see your pull request #97 in the current list of PR there? I certainly hope #97 will be included perhaps as part of #82 on pangeo-data. Thank you for your help.

[RondeauG]

Great! Happy to have been of assistance.

As for pangeo-data, I've seen a few of the discussions regarding the move, but I haven't had the time to look it up in details. If there's a need to re-submit a PR, I'll be happy to do so.

[yuvipanda]

Does pangeo-data/xESMF#1 override this?

[RondeauG]

Hi @yuvipanda. I tested the pangeo-data master around a month ago and as far as I can tell, this Issue has indeed been addressed by the pangeo-data fork: https://github.com/pangeo-data/xESMF

As development of xESMF will continue through pangeo-data, and since this issue was resolved in the pangeo-data fork, I think this issue can be closed.