Open durack1 opened 6 years ago
Yikes! This is very concerning. Until this is resolved consider using regrid2
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.
@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..
UV-CDAT/vcs#274
@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]
@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
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
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()
@durack1 algorithm for gap 0,(-180), 180-0 (westward grid)
For reference the output:
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!)
@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
@dnadeau4 @taylor13 I wonder if this is now fixed in cdms 3.1.2
?
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:![171116_to-ferret](https://user-images.githubusercontent.com/3229632/32909162-b45e941c-caba-11e7-9d5e-9dc1ecc1ba32.png)
The output from
![171116_esmfpat-ferret](https://user-images.githubusercontent.com/3229632/32909199-d029d6d4-caba-11e7-9482-7aa2fb13c3bf.png)
regridMethod='bilinear'
looks likeregridMethod='patch'
looks likeBut![171116_esmfcon-ferret](https://user-images.githubusercontent.com/3229632/32909220-ea08b480-caba-11e7-8528-c2ebc05c65b9.png)
regridMethod='conservative'
is problematic, with missing equatorial values, and a weird striping down the centreAn example script to generate these results is below:
Along with the required data attached below 171116_esmfExample.zip
@dnadeau4 @doutriaux1 @rokuingh any tips would be appreciated