SciTools / iris-esmf-regrid

A collection of structured and unstructured ESMF regridding schemes for Iris.
https://iris-esmf-regrid.readthedocs.io/en/latest
BSD 3-Clause "New" or "Revised" License
20 stars 17 forks source link

Unexpected masked data is introduced in bilinear mesh to grid #391

Open jrackham-mo opened 3 months ago

jrackham-mo commented 3 months ago

šŸ› Bug Report

When regridding from low resolution UGRID data (C4 cubed sphere) to a higher resolution grid (e.g. n216e), masked data points are introduced in the result, despite there being no masked data in the source.

How To Reproduce

Steps to reproduce the behaviour:

  1. Sample C4 data is available here. (see screenshots section for plots)
  2. Sample n216e target grid is available here.

Example script:

>>> from iris.experimental.ugrid import PARSE_UGRID_ON_LOAD
>>> import iris
>>> iris.__version__
'3.9.0'
>>> import esmf_regrid
>>> esmf_regrid.__version__
'0.9.0'
>>> with PARSE_UGRID_ON_LOAD.context():
...     source = iris.load_cube("data_C4.nc", "sample_data")
... 
>>> target = iris.load_cube("non_ugrid_data.nc", "land_area_fraction")
>>> bilinear_regridder = esmf_regrid.schemes.ESMFBilinearRegridder(source, target)
>>> result = bilinear_regridder(source)
>>> import numpy as np
>>> np.ma.is_masked(result.data)
True
>>> np.ma.is_masked(source.data)
False
>>> iris.save(result, "regridded_bilinear_C4_n96e.nc")

Expected behaviour

Expect there to be no masked points in the regridded result. See actual results in screenshots section.

Screenshots

Click to expand this section... ### Source data: C4 cubed sphere ![sample_data_C4](https://github.com/user-attachments/assets/acb3b022-8d45-46b5-81b2-0bbdd1200077) ![sample_data_C4_PlateCarree](https://github.com/user-attachments/assets/9e1dc052-5a55-430a-b9cd-9c0c0afd5c07) ### Regridded data: n96e ![bilinear_regrid_result](https://github.com/user-attachments/assets/f486f34e-d7c0-4fc8-baf8-1cbbc18f0aef) ![bilinear_regrid_result_PlateCarree](https://github.com/user-attachments/assets/ae034ea4-6201-478b-9463-6c0d7f02087f)

Environment

Additional context

We have a use case in UG-ANTS for regridding from a low resolution mesh to a high resolution grid. We would expect there to be no masked data introduced by regridding, or at least for a warning/error to be raised if this does occur.

stephenworsley commented 3 months ago

It looks like this is due to these points lying somewhere in the cracks between bilinear cells. Since these points are considered out of bounds by ESMF, it should be possible to fill these in by passing ESMF a suitable extrapolation method. After a bit of experimentation, I think this would be the NEAREST_IDAVG method with the ESMF Regrid arguments extrap_num_src_pnts=2 and extrap_dist_exponent=1 in order to best approximate the behaviour of bilinear regridding at these points (i.e. taking the weighted average of the two nearest points). Below is a sample of the results different extrapolation methods and arguments passed to ESMF. image

This was generated by hacking iris-esmf-regrid to allow extra arguments to be passed to ESMF with the below code:

esmf_args_1 = {"extrap_method": esmpy.ExtrapMethod.NEAREST_IDAVG,
               "extrap_num_src_pnts": 2,
               "extrap_dist_exponent": 1}
esmf_args_2 = {"extrap_method": esmpy.ExtrapMethod.NEAREST_IDAVG}
esmf_args_3 = {"extrap_method": esmpy.ExtrapMethod.CREEP_FILL,
               "extrap_num_levels": 1}

bilinear_regridder_0 = esmf_regrid.schemes.ESMFBilinearRegridder(source, target)
result_0 = bilinear_regridder_0(source)
bilinear_regridder_1 = esmf_regrid.schemes.ESMFBilinearRegridder(source, target, esmf_args=esmf_args_1)
result_1 = bilinear_regridder_1(source)
bilinear_regridder_2 = esmf_regrid.schemes.ESMFBilinearRegridder(source, target, esmf_args=esmf_args_2)
result_2 = bilinear_regridder_2(source)
bilinear_regridder_3 = esmf_regrid.schemes.ESMFBilinearRegridder(source, target, esmf_args=esmf_args_3)
result_3 = bilinear_regridder_3(source)

plt.subplot(221)
qplt.pcolor(result_0[30:40, 35:45])
plt.title("default")
plt.subplot(222)
qplt.pcolor(result_1[30:40, 35:45])
plt.title("NEAREST_IDAVG (+ args)")
plt.subplot(223)
qplt.pcolor(result_2[30:40, 35:45])
plt.title("NEAREST_IDAVG (no args)")
plt.subplot(224)
qplt.pcolor(result_3[30:40, 35:45])
plt.title("CREEP_FILL")
plt.show()

There are two potential options to solve this:

Changing the default behaviour may have the unintended consequence that out of bounds behaviour changes when dealing with local grids, and there is a more general application for passing ESMF arguments so #375 may be the way to go for solving this. With that said, passing the above suggsted kwargs is a bit cumbersome and non-obvious, so there's still a good case to make for at least providing some sort of shortcut to choose these arguments via the API.