ESMF conservative regrid problems

[durack1]

I've just been trying to regrid a standard lat/lon dataset, and am getting very weird results using the conservative regrid option from ESMF.

The input data looks like:

The output from regridMethod='bilinear' looks like

regridMethod='patch' looks like

But regridMethod='conservative' is problematic, with missing equatorial values, and a weird striping down the centre

An example script to generate these results is below:

import cdms2 as cdm

dataFile = '171116_to.nc'
fH = cdm.open(dataFile)
to = fH('to')
gridFile = '171116_woa13.nc'
fG = cdm.open(gridFile)
t = fG('temperature_oa_mean')
woaGrid = t.getGrid() ; # Get WOA target grid
# Try regrid methods
diagBilin = {}
esmfBilin = to.regrid(woaGrid,regridTool='ESMF',regridMethod='bilinear',
                            srcMaskValues=to.missing,
                            dstMaskValues=1e20,diag=diagBilin)
diagCon = {}
esmfCon = to.regrid(woaGrid,regridTool='ESMF',regridMethod='conservative',
                            srcMaskValues=to.missing,
                            dstMaskValues=1e20,diag=diagCon)
diagPat = {}
esmfPat = to.regrid(woaGrid,regridTool='ESMF',regridMethod='patch',
                            srcMaskValues=to.missing,
                            dstMaskValues=1e20,diag=diagPat)

Along with the required data attached below
171116_esmfExample.zip

@dnadeau4 @doutriaux1 @rokuingh any tips would be appreciated

[gleckler1]

Yikes! This is very concerning. Until this is resolved consider using regrid2

[taylor13]

I'm guessing ESMF can easily handle a case where you go from one latxlon grid to another with (or without) masking, so it's either a problem with CDAT or perhaps the input is incorrectly configured (or wrong). Does the destination grid input have (or need to have) a "ghost" longitude, padding the eastern-most and/or western-most grid cells? Might account for the problem at 180 deg. Harder, perhaps, to explain the problem at the equator.

[durack1]

@taylor13 the longitude values are

Out[13]: 
array([   0.5,    1.5,    2.5,    3.5,    4.5,    5.5,    6.5,    7.5,
...
         -7.5,   -6.5,   -5.5,   -4.5,   -3.5,   -2.5,   -1.5,   -0.5], dtype=float32)

So these seem fine to me..

CDAT/vcs#274

[dnadeau4]

@durack1 You don't use the mask right. You should have an array of mask values you want to activate. Since your mask is "True" or "False" you need to set the mask to 1.

srcMaskValuesArr = numpy.array([1], dtype=numpy.int32)

[srcMaskValues]

• Mask information can be set in the Grid (see 31.3.16) or Mesh (see 33.3.9) upon which the srcField is built. The srcMaskValues argument specifies the values in that mask information which indicate a source point should be masked out. In other words, a location is masked if and only if the value for that location in the mask information matches one of the values listed in srcMaskValues. If srcMaskValues is not specified, no masking will occur.

[durack1]

@dnadeau4 is the URL of where you obtained the info above available? A major issue for users is that none of the documentation for ESMF is obviously available through inline documentation or web docs

[dnadeau4]

Thanks @durack1. It was an email between Ben and I. Here is the reference I found.
http://www.earthsystemmodeling.org/esmf_releases/last_built/ESMF_refdoc/node5.html#SECTION050366000000000000000

[dnadeau4]

Using your example, I was able to convert the code in native format. I did not have time to convert to conservative more.

import ESMF
import cdms2
import cdat_info
import numpy
import vcs
dataFile = '171116_to.nc'
fH = cdms2.open(dataFile)
so = fH('to')[:,:]

gridFile = '171116_woa13.nc'
fG = cdms2.open(gridFile)
clt = fG('temperature_oa_mean')[:,:]
numpy.save("sodata", so.data[:], allow_pickle=False)
#numpy.save("somask", so.mask[:], allow_pickle=False)
numpy.save("solat", so.getLatitude()[:], allow_pickle=False)
numpy.save("solon", so.getLongitude()[:], allow_pickle=False)
sodata = numpy.load("sodata.npy")
#somask = numpy.load("somask.npy")
solat = numpy.load("solat.npy")
solon = numpy.load("solon.npy")
so = numpy.ma.masked_array(sodata, mask=False)

clt = cdms2.open(cdat_info.get_sampledata_path() + '/clt.nc')('clt')[0, :, :] 
numpy.save("cltdata", clt.data[:], allow_pickle=False)
numpy.save("cltlat", clt.getLatitude()[:], allow_pickle=False)
numpy.save("cltlon", clt.getLongitude()[:], allow_pickle=False)
cltdata = numpy.load("cltdata.npy")
cltlat = numpy.load("cltlat.npy")
cltlon = numpy.load("cltlon.npy")
clt = numpy.ma.masked_array(cltdata, mask=False)
#periodic_dim=periodic_dim, pole_dim=pole_dim,

print(so.shape)
maxIndex = numpy.array(so.shape, dtype=numpy.int32)
srcGrid = ESMF.Grid(max_index=maxIndex, num_peri_dims=1, periodic_dim=1, pole_dim=0, staggerloc=[ESMF.StaggerLoc.CENTER], coord_sys=ESMF.CoordSys.SPH_DEG)
print(clt.shape)
maxIndex = numpy.array(clt.shape, dtype=numpy.int32)
dstGrid = ESMF.Grid(max_index=maxIndex, num_peri_dims=1, periodic_dim=1, pole_dim=0, staggerloc=[ESMF.StaggerLoc.CENTER], coord_sys=ESMF.CoordSys.SPH_DEG)



ny, nx = so.shape
y = solat
#x = clt.getGrid().getLongitude()
x = solon
yy = numpy.outer(y, numpy.ones((nx,), numpy.float32))
xx = numpy.outer(numpy.ones((ny,), numpy.float32), x)

ptr = srcGrid.get_coords(coord_dim=0, staggerloc=ESMF.StaggerLoc.CENTER)
#ptr[:] =so.getGrid().getLongitude()
ptr[:] =xx
ptr = srcGrid.get_coords(coord_dim=1, staggerloc=ESMF.StaggerLoc.CENTER)
#ptr[:] = so.getGrid().getLatitude()
ptr[:] = yy

ny, nx = clt.shape
#y = clt.getGrid().getLatitude()
y = cltlat
#x = clt.getGrid().getLongitude()
x = cltlon
yy = numpy.outer(y, numpy.ones((nx,), numpy.float32))
xx = numpy.outer(numpy.ones((ny,), numpy.float32), x)

ptr = dstGrid.get_coords(coord_dim=0, staggerloc=ESMF.StaggerLoc.CENTER)
ptr[:] = xx
ptr = dstGrid.get_coords(coord_dim=1, staggerloc=ESMF.StaggerLoc.CENTER)
ptr[:] = yy

mask = numpy.zeros(so.shape, numpy.int32)
mask[:] = numpy.ma.masked_where((so == 1.0e20), so).mask
srcGrid.add_item(item=ESMF.GridItem.MASK, staggerloc=ESMF.StaggerLoc.CENTER)
maskPtr = srcGrid.get_item( item=ESMF.GridItem.MASK, staggerloc=ESMF.StaggerLoc.CENTER)
lo = srcGrid.lower_bounds[ESMF.StaggerLoc.CENTER]
hi = srcGrid.upper_bounds[ESMF.StaggerLoc.CENTER]
slab = tuple([slice(lo[i], hi[i], None) for i in range(2)])
print(slab)
maskPtr[:] = mask[slab]

srcFld = ESMF.Field( grid=srcGrid, name="srcFld", typekind=ESMF.TypeKind.R8, staggerloc=ESMF.StaggerLoc.CENTER)
ptr = srcFld.data
ptr[:] = numpy.array(so[:])

dstFld = ESMF.Field( grid=dstGrid, name='dstFld', typekind=ESMF.TypeKind.R8, staggerloc=ESMF.StaggerLoc.CENTER)
ptr = dstFld.data
# set to missing_value
#ptr[:] = numpy.array(1.e20 * numpy.ones(clt.shape, numpy.float32))
ptr[:] = numpy.array(0 * numpy.ones(clt.shape, numpy.float64))

srcFracField= ESMF.Field( grid=srcGrid, name='src_cell_area_fractions', typekind=ESMF.TypeKind.R8, staggerloc=ESMF.StaggerLoc.CENTER)
dataPtr = srcFracField.data
dataPtr[:] = 1.0

dstFracField= ESMF.Field( grid=dstGrid, name='dst_cell_area_fractions', typekind=ESMF.TypeKind.R8, staggerloc=ESMF.StaggerLoc.CENTER)
dataPtr = dstFracField.data
dataPtr[:] = 1.0

srcMaskValuesArr = numpy.array([1], dtype=numpy.int32)
dstMaskValuesArr = numpy.array([1], dtype=numpy.int32)

regridHandle = ESMF.Regrid(srcFld, dstFld, 
                           src_frac_field=srcFracField,
                           dst_frac_field=dstFracField,
                           src_mask_values=srcMaskValuesArr,
                           dst_mask_values=dstMaskValuesArr,
                           regrid_method = ESMF.RegridMethod.BILINEAR,
                           unmapped_action = ESMF.UnmappedAction.IGNORE,
                           ignore_degenerate= False) 

regridHandle(srcfield=srcFld, dstfield=dstFld, zero_region=ESMF.Region.SELECT)

#numpy.save("dstFld", dstFld.data[:],allow_pickle=False)
print dstFld.data.min()
print dstFld.data.max()
import pdb
pdb.set_trace()
ppp=vcs.init()
b=ppp.createboxfill()
b.missing="grey"
b.boxfill_type="custom"
b.levels = range(0,30,1)
b.fillareacolors= vcs.getcolors(b.levels)
mask = numpy.ma.masked_where((dstFld.data == 1.0e20), dstFld.data).mask
out=cdms2.createVariable(dstFld.data,mask=mask)
ppp.plot(out,b)
ppp.interact()

[dnadeau4]

@durack1 algorithm for gap 0,(-180), 180-0 (westward grid)

For reference the output:

[dnadeau4]

esmfissue.zip
longitude bounds are wrong for some reasons. -180 and 180 is missing.....
[-179., -178.],
[ 0., -179.],
[ 179., 0.],
[ 178., 179.],

That would mess up the regridder.

I got these 2 cases to work!!!

x = numpy.arange(359.5, -0.5,-1)
tolonbnds = numpy.zeros((x.size, 2))
tolonbnds[:, 1] = x - 0.5  ## note  need to reverse order (left to right) [360, 359]
tolonbnds[:, 0] = x + 0.5
(no lines everything works!)


x = numpy.arange(0.5,360,1)
tolonbnds = numpy.zeros((x.size, 2))
tolonbnds[:, 0] = x - 0.5   ## note order (right to left)  is [359 360]
tolonbnds[:, 1] = x + 0.5
(no line everything is fine!)

[durack1]

@dnadeau4 @doutriaux1 @gleckler1 @taylor13 I’ll reopen this to double check (regenerate the outputs above) to make sure there’s not a conservative issue that has crept in here

[durack1]

@dnadeau4 @taylor13 I wonder if this is now fixed in cdms 3.1.2?