Conservative remapping yields ESMC_FieldRegridStoreFile

I'm working on conservative remapping of air-sea CO2 fluxes on ocean grids. The goal is to remap from an arbitrary curvilinear grid to a 180x360 rectilinear grid. I've followed advice from https://github.com/JiaweiZhuang/xESMF/issues/14 and https://github.com/JiaweiZhuang/xESMF/issues/32, but am getting stuck.

As mentioned in those issues, typical ocean output has lat/lon bounds on (N, M, 4) dimensions, where 4 are the grid cell corners. I created a function following @JiaweiZhuang's advice (https://github.com/JiaweiZhuang/xESMF/issues/14#issuecomment-369686779) to convert from that convention to (N+1, M+1):

def _unravel(new_bounds, vertex_bounds, M, N):
    """Unravel vertices to (N+1, M+1) dimensions"""
    new_bounds[0:N, 0:M] = vertex_bounds[:, :, 0]

    # fill in missing row
    new_bounds[N, 0:M] = vertex_bounds[N-1, :, 1]
    # fill in missing column
    new_bounds[0:N, M] = vertex_bounds[:, M-1, 2]
    # fill in remaining element
    new_bounds[N, M] = vertex_bounds[N-1, M-1, 3]
    return new_bounds

def compress_vertices(ds, lat_bnds_name='lat_b', lon_bnds_name='lon_b'):
    """Convert from (N, M, nvertices) to (N+1, M+1) for grid cell corners

    Args:
        ds (xarray Dataset): Dataset with native ocean grid and (N, M, 4) 
                                          lat/lon bounds.
        lat_bnds_name (optional str): Name of lat bounds variable.
        lon_bnds_name (optional str): Name of lon bounds variable
    """
    M = ds.x.size
    N = ds.y.size

    # create arrays for 2D bounds info
    lat_b = np.zeros((N+1, M+1))
    lon_b = np.zeros((N+1, M+1))

    # unravel nvertices to 2D style
    lat_b = _unravel(lat_b, ds[lat_bnds_name], M, N)
    lon_b = _unravel(lon_b, ds[lon_bnds_name], M, N)

    # get rid of old coordinates
    del ds[lat_bnds_name], ds[lon_bnds_name]
    ds = ds.squeeze()

    # assign new coordinates
    ds.coords['lat_b'] = (('y_b', 'x_b'), lat_b)
    ds.coords['lon_b'] = (('y_b', 'x_b'), lon_b)
    return ds

When I feed this into the xESMF regridder with 'conservative', I get the following error:

ValueError: ESMC_FieldRegridStoreFile() failed with rc = 506.

The log file isn't very enlightening here:

20190813 094016.108 ERROR            PET0 ESMF_Regrid.F90:360 ESMF_RegridStore Invalid argument  - Internal subroutine call returned Error
20190813 094016.109 ERROR            PET0 ESMF_FieldRegrid.F90:1292 ESMF_FieldRegridStoreNX Invalid argument  - Internal subroutine call returned Error
20190813 094016.109 ERROR            PET0 ESMF_Field_C.F90:1168 f_esmf_regridstorefile Invalid argument  - Internal subroutine call returned Error

However, my modified dataset with lat/lon bounds works with bilinear remapping. I've gotten conservative remapping to work with NCL with the native corners, but am looking to use xESMF/python exclusively.

Any help here is appreciated. Sorry for duplicate issues on conservative remapping. I'd be happy to help implement some utility functions to help out with this for future releases, since it seems to be a persistent issue.

Here is a notebook with my attempt: https://nbviewer.jupyter.org/gist/bradyrx/421627385666eefdb0a20567c2da9976

Here is a dropbox with the notebook and the native CESM2 output: https://www.dropbox.com/sh/lef7i2v88sfke3d/AACNPgjIrXBV35ejLYTRFvqxa?dl=0

Thanks for the complete code!

The log file isn't very enlightening here:

You can make the ESMF log file more informative by setting

import ESMF
ESMF.Manager(debug=True)

The log shows:

20190813 223533.907 INFO             PET0 Running with ESMF Version 7.1.0r
20190813 223544.261 ERROR            PET0 ~~~~~~~~~~~~~~~~~~~~ Degenerate Element Detected ~~~~~~~~~~~~~~~~~~~~
20190813 223544.261 ERROR            PET0   degenerate elem. id=122561
20190813 223544.261 ERROR            PET0 
20190813 223544.261 ERROR            PET0   degenerate elem. coords (lon [-180 to 180], lat [-90 to 90]) (x,y,z)
20190813 223544.261 ERROR            PET0   -----------------------------------------------------------------
20190813 223544.261 ERROR            PET0     0  (-39.999603,  71.954628)  (0.237299, -0.199115, 0.950812)
20190813 223544.261 ERROR            PET0     1  (-39.548401,  71.959839)  (0.238793, -0.197185, 0.950840)
20190813 223544.261 ERROR            PET0     2  (-39.097107,  71.965042)  (0.240272, -0.195243, 0.950868)
20190813 223544.261 ERROR            PET0     3  (-39.548401,  71.959839)  (0.238793, -0.197185, 0.950840)
20190813 223544.261 ERROR            PET0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
20190813 223544.261 ERROR            PET0 ESMCI_Mesh_Regrid_Glue.C:161 c_esmc_regrid_create() Invalid argument  - - Src contains a cell that has corners close enough that the cell collapses to a line or point
20190813 223544.261 ERROR            PET0 ESMCI_Mesh_Regrid_Glue.C:512 c_esmc_regrid_create() Invalid argument  - Internal subroutine call returned Error
...

It says that the source grid contains a cell that has corners close enough that the cell collapses to a line or point. The location is around (lon=-40, lat=72). You can overlay this point with your input grid by:

plt.scatter(ds['lon_b'], ds['lat_b'], s=0.1, alpha=0.2)  # plot the input grid mesh
plt.plot(-40+360,  72, 'ro')  # plot the problematic point
plt.xlabel('lon')
plt.ylabel('lat')

which shows:

Because the grid cell near the red "pole" has very close corners (you can see two corners are exactly (-39.548401, 71.959839), from the log), ESMF thinks that it is a triangle instead of a quadrilateral (i.e. "degenerated") and throws out this error.

I think this error can be skipped by setting ignore_degenerate=True when calling ESMF.Regrid, as suggested by uturuncoglu/RegESM#15 (not configurable in xesmf yet, will test it). Alternatively, you can manually remove such degenerated cells.

JiaweiZhuang / xESMF

Conservative remapping yields ESMC_FieldRegridStoreFile #60