Closed durack1 closed 9 years ago
Can you also send us the latitudes? Or point me to the file.
Thanks
From: "Paul J. Durack" notifications@github.com<mailto:notifications@github.com> Reply-To: UV-CDAT/uvcdat reply@reply.github.com<mailto:reply@reply.github.com> Date: Friday, March 6, 2015 at 4:51 PM To: UV-CDAT/uvcdat uvcdat@noreply.github.com<mailto:uvcdat@noreply.github.com> Subject: [uvcdat] ESMF regrid problem (#1125)
From: , Mark Date: Wednesday, March 4, 2015 at 11:40 AM To: "Doutriaux, Charles", "Paul J. Durack" Subject: CDAT regridder issue
Hi Charles and Paul,
Have you guys ever encountered the following issue?
If I regrid using SWCRE_grd=SWCRE.regrid(horiz_grid,regridTool="regrid2”), everything looks okay (see attached top right figure).
But when I regrid using SWCRE_grd=SWCRE.regrid(horiz_grid,regridTool="esmf",regridMethod = "linear”), things get messed up, but only in the NH. Note how things look okay in the SH. (See attached bottom right figure.)
My grids seem to be fine: Original longitudes:
SWCRE.getLongitude()[:] Out[39]: array([ 135.5, 136.5, 137.5, 138.5, 139.5, 140.5, 141.5, 142.5, 143.5, 144.5, 145.5, 146.5, 147.5, 148.5, 149.5, 150.5, 151.5, 152.5, 153.5, 154.5, 155.5, 156.5, 157.5, 158.5, 159.5, 160.5, 161.5, 162.5, 163.5, 164.5, 165.5, 166.5, 167.5, 168.5, 169.5, 170.5, 171.5, 172.5, 173.5, 174.5, 175.5, 176.5, 177.5, 178.5, 179.5, 180.5, 181.5, 182.5, 183.5, 184.5, 185.5, 186.5, 187.5, 188.5, 189.5, 190.5, 191.5, 192.5, 193.5, 194.5, 195.5, 196.5, 197.5, 198.5, 199.5, 200.5, 201.5, 202.5, 203.5, 204.5, 205.5, 206.5, 207.5, 208.5, 209.5, 210.5, 211.5, 212.5, 213.5, 214.5, 215.5, 216.5, 217.5, 218.5, 219.5, 220.5, 221.5, 222.5, 223.5, 224.5, 225.5, 226.5, 227.5, 228.5, 229.5, 230.5, 231.5, 232.5, 233.5, 234.5], dtype=float32)
Destination longitudes:
horiz_grid.getLongitude()[:] Out[38]: array([ 135., 137., 139., 141., 143., 145., 147., 149., 151., 153., 155., 157., 159., 161., 163., 165., 167., 169., 171., 173., 175., 177., 179., 181., 183., 185., 187., 189., 191., 193., 195., 197., 199., 201., 203., 205., 207., 209., 211., 213., 215., 217., 219., 221., 223., 225., 227., 229., 231., 233., 235.])
Any ideas?
Mark [swcre_regrid_2options]https://cloud.githubusercontent.com/assets/3229632/6539096/fce3ac0c-c420-11e4-8194-e27242995741.png
— Reply to this email directly or view it on GitHubhttps://github.com/UV-CDAT/uvcdat/issues/1125.
1) Here is what I am using (on feedback): $ which cdat /usr/local/uvcdat/2015-02-10/bin/cdat $ which uvcdat uvcdat: aliased to /usr/local/uvcdat/2013-06-05/bin/uvcdat $ which ipython /usr/local/uvcdat/2015-02-10/bin/ipython
2) Here are the latitudes of the two grids: Original latitudes: SWCRE.getLatitude()[:] Out[3]: array([-29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5, -20.5, -19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5], dtype=float32)
Destination latitudes: horiz_grid.getLatitude()[:] Out[4]: array([-29., -27., -25., -23., -21., -19., -17., -15., -13., -11., -9., -7., -5., -3., -1., 1., 3., 5., 7., 9., 11., 13., 15., 17., 19., 21., 23., 25., 27., 29.])
Thanks!
The fact that your uvcdat is different is worrisome (mixed versions) I will try to reproduce next week
I just changed to this uvcdat: /usr/local/uvcdat/2015-02-10/bin/uvcdat and the problem remains.
The following code reproduces the error:
import vcs,cdms2
import os,sys
f=cdms2.open("/home/doutriaux1/Mark/charles.nc")
s=f("swcre")
lon_dest = cdms2.createAxis([ 135., 137., 139., 141., 143., 145., 147., 149., 151.,
153., 155., 157., 159., 161., 163., 165., 167., 169.,
171., 173., 175., 177., 179., 181., 183., 185., 187.,
189., 191., 193., 195., 197., 199., 201., 203., 205.,
207., 209., 211., 213., 215., 217., 219., 221., 223.,
225., 227., 229., 231., 233., 235.])
lon_dest.designateLongitude()
lon_dest.id="longitude"
lon_dest.units="degrees_east"
lat_dest = cdms2.createAxis([-29., -27., -25., -23., -21., -19., -17., -15., -13., -11., -9.,
-7., -5., -3., -1., 1., 3., 5., 7., 9., 11., 13.,
15., 17., 19., 21., 23., 25., 27., 29.])
lat_dest.designateLatitude()
lat_dest.units="degrees_north"
dummy = cdms2.MV2.ones((len(lat_dest),len(lon_dest)))
dummy.setAxisList((lat_dest,lon_dest))
grid_dest=dummy.getGrid()
s.id="orig"
s_regrid2 = s.regrid(grid_dest,regridTool="regrid2")
s_regrid2.id="regrid2"
s_esmf_lin = s.regrid(grid_dest)
s_esmf_lin.id = "ESMF Linear"
s_esmf_con = s.regrid(grid_dest,regridTool="esmf",regridMethod="conservative")
s_esmf_lin.id = "ESMF Conservative"
import EzTemplate
x=vcs.init()
t=x.createtemplate()
t.blank()
t.data.priority=1
t.legend.priority=1
t.dataname.priority=1
t.dataname.y=t.dataname.y*.95
M=EzTemplate.Multi(template=t,x=x,rows=2,columns=2)
gm=x.createboxfill()
levels= vcs.mkscale(*vcs.minmax(s))
cols = vcs.getcolors(levels)
gm.boxfill_type = "custom"
gm.fillareacolor = cols
gm.levels = levels
x.plot(s,M.get(),gm)
x.plot(s_regrid2,M.get(),gm)
x.plot(s_esmf_lin,M.get(),gm)
x.plot(s_esmf_con,M.get(),gm)
x.png("esmf_issue_1125")
raw_input("Press enter")
And here is the output:
@dnadeau4 @doutriaux1 just pinging you both on this..
it's on my list for 2.4... keep your finger crossed....
@dnadeau4 this is the answer I got from Ryan OKuinghttons
Hi Charles,
I've got some good news and some bad news for you. The good news is
that I've found the issue, and the bad news is that it was something in
my ESMP reproducer script. This means that your problem is likely
caused by the wrapper layer between ESMP and UVCDAT. The mistake I made
in the reproducer was caused by creating a periodic Grid, instead of a
regional Grid. This is easy to fix by replacing:
grid = ESMP.ESMP_GridCreate1PeriDim(maxIndex)
with:
grid = ESMP.ESMP_GridCreateNoPeriDim(maxIndex)
After this change, and a few other minor mods to clean up the plots and
make the two reproducer codes more alike, the results are as expected.
I attached the ESMP and ESMPy reproducer codes, you can easily reproduce
this error by reversing the change I mentioned above. There is also one
minor error in the ESMPy script that I still need to fix, but I wanted
to get you a more satisfactory solution as soon as possible.
If I can make a guess, I would say that the error in the interface layer
(as well as the one I made in the reproducer) is due to the fact that
the latitude data in this example _appears_ to cross a pole (at 0
degrees). I suppose that the grid in this example is intended to have
the poles at -90 and 90 degrees. Either way, I think that the structure
of this data is throwing off the interface layer, causing it to create a
Grid with a periodic dimension. Again, this is pure speculation on my part.
Also, when talking to Ben today I realized that I forgot to mention in
my last mail that ESMPy now also interfaces with OCGIS using converters
between ESMPy Fields and OCGIS Variables. The current capability is
limited to IO and subsetting in serial, but we plan to expand that in
the future to include computations and parallel execution.
Thanks!
Ryan
Ryan scripts: uvcdat_error_esmp.py
# This is an ESMP reproducer of uvcdat bug 1125, esmf support request 3613723
import ESMP
import numpy as np
def make_index( nx, ny ):
return np.arange( nx*ny ).reshape( ny, nx )
def create_grid(lons, lats, lonbnds, latbnds):
maxIndex = np.array([len(lons), len(lats)], dtype=np.int32)
grid = ESMP.ESMP_GridCreateNoPeriDim(maxIndex)
## CENTERS
ESMP.ESMP_GridAddCoord(grid, staggerloc=ESMP.ESMP_STAGGERLOC_CENTER)
exLB_center, exUB_center = ESMP.ESMP_GridGetCoord(grid,
ESMP.ESMP_STAGGERLOC_CENTER)
# get the coordinate pointers and set the coordinates
[x, y] = [0, 1]
gridXCenter = ESMP.ESMP_GridGetCoordPtr(grid, x, ESMP.ESMP_STAGGERLOC_CENTER)
gridYCenter = ESMP.ESMP_GridGetCoordPtr(grid, y, ESMP.ESMP_STAGGERLOC_CENTER)
# make an array that holds indices from lower_bounds to upper_bounds
bnd2indX = np.arange(exLB_center[x], exUB_center[x], 1)
bnd2indY = np.arange(exLB_center[y], exUB_center[y], 1)
pts = make_index(exUB_center[x] - exLB_center[x],
exUB_center[y] - exLB_center[y])
for i0 in range(len(bnd2indX)):
gridXCenter[pts[:, i0]] = lons[i0]
for i1 in range(len(bnd2indY)):
gridYCenter[pts[i1, :]] = lats[i1]
## CORNERS
ESMP.ESMP_GridAddCoord(grid, staggerloc=ESMP.ESMP_STAGGERLOC_CORNER)
exLB_corner, exUB_corner = ESMP.ESMP_GridGetCoord(grid, ESMP.ESMP_STAGGERLOC_CORNER)
# get the coordinate pointers and set the coordinates
[x,y] = [0, 1]
gridXCorner = ESMP.ESMP_GridGetCoordPtr(grid, x, ESMP.ESMP_STAGGERLOC_CORNER)
gridYCorner = ESMP.ESMP_GridGetCoordPtr(grid, y, ESMP.ESMP_STAGGERLOC_CORNER)
# make an array that holds indices from lower_bounds to upper_bounds
bnd2indX = np.arange(exLB_corner[x], exUB_corner[x], 1)
bnd2indY = np.arange(exLB_corner[y], exUB_corner[y], 1)
pts = make_index(exUB_corner[x] - exLB_corner[x],
exUB_corner[y] - exLB_corner[y])
for i0 in range(len(bnd2indX)-1):
gridXCorner[pts[:, i0]] = lonbnds[i0, 0]
gridXCorner[pts[:, i0+1]] = lonbnds[i0, 1]
for i1 in range(len(bnd2indY)-1):
gridYCorner[pts[i1, :]] = latbnds[i1, 0]
gridYCorner[pts[i1+1, :]] = latbnds[i1, 1]
return grid
def build_analyticfield(field, grid):
# get the field pointer
fieldPtr = ESMP.ESMP_FieldGetPtr(field)
# get the grid bounds and coordinate pointers
exLB, exUB = ESMP.ESMP_GridGetCoord(grid, ESMP.ESMP_STAGGERLOC_CENTER)
# make an array that holds indices from lower_bounds to upper_bounds
[x,y] = [0, 1]
bnd2indX = np.arange(exLB[x], exUB[x], 1)
bnd2indY = np.arange(exLB[y], exUB[y], 1)
realdata = False
try:
import netCDF4 as nc
f = nc.Dataset('charles.nc')
swcre = f.variables['swcre']
swcre = swcre[:]
realdata = True
except:
raise ImportError('netCDF4 not available on this machine')
p = 0
for i1 in range(len(bnd2indX)):
for i0 in range(len(bnd2indY)):
if realdata:
fieldPtr[p] = swcre.flatten()[p]
else:
fieldPtr[p] = 2.0
p = p + 1
return field
def compute_mass(valuefield, areafield, fracfield, dofrac):
mass = 0.0
ESMP.ESMP_FieldRegridGetArea(areafield)
area = ESMP.ESMP_FieldGetPtr(areafield)
value = ESMP.ESMP_FieldGetPtr(valuefield)
frac = 0
if dofrac:
frac = ESMP.ESMP_FieldGetPtr(fracfield)
for i in range(valuefield.size):
if dofrac:
mass += area[i]*value[i]*frac[i]
else:
mass += area[i]*value[i]
return mass
def plot(srclons, srclats, srcfield, dstlons, dstlats, interpfield):
try:
import matplotlib
import matplotlib.pyplot as plt
except:
raise ImportError("matplotlib is not available on this machine")
exact = ESMP.ESMP_FieldGetPtr(srcfield)
exact = exact.reshape(len(srclats.flatten()), len(srclons.flatten()))
interp = ESMP.ESMP_FieldGetPtr(interpfield)
interp = interp.reshape(len(dstlats.flatten()), len(dstlons.flatten()))
fig = plt.figure(1, (15, 6))
fig.suptitle('ESMP Conservative Regridding', fontsize=14, fontweight='bold')
ax = fig.add_subplot(1, 2, 1)
im = ax.imshow(exact, vmin=-140, vmax=0, cmap='gist_ncar', aspect='auto',
extent=[min(srclons), max(srclons), min(srclats), max(srclats)])
ax.set_xbound(lower=min(srclons), upper=max(srclons))
ax.set_ybound(lower=min(srclats), upper=max(srclats))
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.set_title("Source Data")
ax = fig.add_subplot(1, 2, 2)
im = ax.imshow(interp, vmin=-140, vmax=0, cmap='gist_ncar', aspect='auto',
extent=[min(dstlons), max(dstlons), min(dstlats), max(dstlats)])
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.set_title("Conservative Regrid Solution")
uvcdat_error_esmpy.py
# This is an ESMPy reproducer of uvcdat bug 1125, esmf support request 3613723
import ESMF
import numpy as np
def create_grid_corners(lons, lats, lonbnds, latbnds):
[lon, lat] = [0, 1]
max_index = np.array([len(lons), len(lats)])
grid = ESMF.Grid(max_index,
staggerloc=[ESMF.StaggerLoc.CENTER, ESMF.StaggerLoc.CORNER])
gridXCenter = grid.get_coords(lon)
lon_par = lons[grid.lower_bounds[ESMF.StaggerLoc.CENTER][lon]:grid.upper_bounds[ESMF.StaggerLoc.CENTER][lon]]
gridXCenter[...] = lon_par.reshape((lon_par.size, 1))
gridYCenter = grid.get_coords(lat)
lat_par = lats[grid.lower_bounds[ESMF.StaggerLoc.CENTER][lat]:grid.upper_bounds[ESMF.StaggerLoc.CENTER][lat]]
gridYCenter[...] = lat_par.reshape((1, lat_par.size))
lbx = grid.lower_bounds[ESMF.StaggerLoc.CORNER][lon]
ubx = grid.upper_bounds[ESMF.StaggerLoc.CORNER][lon]
lby = grid.lower_bounds[ESMF.StaggerLoc.CORNER][lat]
uby = grid.upper_bounds[ESMF.StaggerLoc.CORNER][lat]
gridXCorner = grid.get_coords(lon, staggerloc=ESMF.StaggerLoc.CORNER)
for i0 in range(ubx - lbx - 1):
gridXCorner[i0, :] = lonbnds[i0, 0]
gridXCorner[i0 + 1, :] = lonbnds[i0, 1]
gridYCorner = grid.get_coords(lat, staggerloc=ESMF.StaggerLoc.CORNER)
for i1 in range(uby - lby - 1):
gridYCorner[:, i1] = latbnds[i1, 0]
gridYCorner[:, i1 + 1] = latbnds[i1, 1]
return grid
def initialize_field(field):
realdata = False
try:
import netCDF4 as nc
f = nc.Dataset('charles.nc')
swcre = f.variables['swcre']
swcre = swcre[:]
realdata = True
except:
raise ImportError('netCDF4 not available on this machine')
if realdata:
# transpose because uvcdat data is represented as lat/lon
field.data[...] = swcre.T
else:
field.data[...] = 2.0
return field
def compute_mass(valuefield, areafield, fracfield, dofrac):
mass = 0.0
areafield.get_area()
if dofrac:
mass = np.sum(areafield*valuefield*fracfield)
else:
mass = np.sum(areafield * valuefield)
return mass
def plot(srclons, srclats, srcfield, dstlons, dstlats, interpfield):
try:
import matplotlib
import matplotlib.pyplot as plt
except:
raise ImportError("matplotlib is not available on this machine")
fig = plt.figure(1, (15, 6))
fig.suptitle('ESMPy Conservative Regridding', fontsize=14, fontweight='bold')
ax = fig.add_subplot(1, 2, 1)
im = ax.imshow(srcfield.T, vmin=-140, vmax=0, cmap='gist_ncar', aspect='auto',
extent=[min(srclons), max(srclons), min(srclats), max(srclats)])
ax.set_xbound(lower=min(srclons), upper=max(srclons))
ax.set_ybound(lower=min(srclats), upper=max(srclats))
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.set_title("Source Data")
ax = fig.add_subplot(1, 2, 2)
im = ax.imshow(interpfield.T, vmin=-140, vmax=0, cmap='gist_ncar', aspect='auto',
extent=[min(dstlons), max(dstlons), min(dstlats), max(dstlats)])
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.set_title("Conservative Regrid Solution")
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.9, 0.1, 0.01, 0.8])
fig.colorbar(im, cax=cbar_ax)
plt.show()
##########################################################################################
# Start up ESMF, this call is only necessary to enable debug logging
esmpy = ESMF.Manager(logkind=ESMF.LogKind.MULTI, debug=True)
# Create a destination grid from a GRIDSPEC formatted file.
srcgrid = ESMF.Grid(filename="charles.nc",
filetype=ESMF.FileFormat.GRIDSPEC, add_corner_stagger=True, is_sphere=False)
try:
import netCDF4 as nc
except:
raise ImportError('netCDF4 not available on this machine')
f = nc.Dataset('charles.nc')
srclons = f.variables['lon'][:]
srclats = f.variables['lat'][:]
srclonbounds = f.variables['bounds_lon'][:]
srclatbounds = f.variables['bounds_lat'][:]
#srcgrid = create_grid_corners(srclons, srclats, srclonbounds, srclatbounds)
# original destination center from charles
# dstlons = numpy.array([135., 137., 139., 141., 143., 145., 147., 149., 151.,
# 153., 155., 157., 159., 161., 163., 165., 167., 169.,
# 171., 173., 175., 177., 179., 181., 183., 185., 187.,
# 189., 191., 193., 195., 197., 199., 201., 203., 205.,
# 207., 209., 211., 213., 215., 217., 219., 221., 223.,
# 225., 227., 229., 231., 233., 235.])
# dstlats = numpy.array([-29., -27., -25., -23., -21., -19., -17., -15., -13., -11., -9.,
# -7., -5., -3., -1., 1., 3., 5., 7., 9., 11., 13.,
# 15., 17., 19., 21., 23., 25., 27., 29.])
# create data slightly offset for a destination grid
dstlats = srclats - 0.5
dstlons = srclons - 0.5
dstlatbounds = srclatbounds - 0.5
dstlonbounds = srclonbounds - 0.5
dstgrid = create_grid_corners(dstlons, dstlats, dstlonbounds, dstlatbounds)
srcfield = ESMF.Field(srcgrid, "srcfield", staggerloc=ESMF.StaggerLoc.CENTER)
dstfield = ESMF.Field(dstgrid, "dstfield", staggerloc=ESMF.StaggerLoc.CENTER)
srcareafield = ESMF.Field(srcgrid, "srcfield", staggerloc=ESMF.StaggerLoc.CENTER)
dstareafield = ESMF.Field(dstgrid, "dstfield", staggerloc=ESMF.StaggerLoc.CENTER)
srcfracfield = ESMF.Field(srcgrid, "srcfield", staggerloc=ESMF.StaggerLoc.CENTER)
dstfracfield = ESMF.Field(dstgrid, "dstfield", staggerloc=ESMF.StaggerLoc.CENTER)
srcfield = initialize_field(srcfield)
# Regrid from source grid to destination grid.
regridSrc2Dst = ESMF.Regrid(srcfield, dstfield,
regrid_method=ESMF.RegridMethod.CONSERVE,
unmapped_action=ESMF.UnmappedAction.ERROR,
src_frac_field=srcfracfield,
dst_frac_field=dstfracfield)
dstfield = regridSrc2Dst(srcfield, dstfield)
srcmass = compute_mass(srcfield, srcareafield, srcfracfield, True)
dstmass = compute_mass(dstfield, dstareafield, 0, False)
print "Conservative error = {}".format(abs(srcmass-dstmass)/abs(srcmass))
plot(srclons, srclats, srcfield, dstlons, dstlats, dstfield)
print '\nUVCDAT example completed successfully.\n'
From:, Mark
Date: Wednesday, March 4, 2015 at 11:40 AM
To: "Doutriaux, Charles", "Paul J. Durack"
Subject: CDAT regridder issue
Hi Charles and Paul,
Have you guys ever encountered the following issue?
If I regrid using SWCRE_grd=SWCRE.regrid(horiz_grid,regridTool="regrid2”), everything looks okay (see attached top right figure).
But when I regrid using SWCRE_grd=SWCRE.regrid(horiz_grid,regridTool="esmf",regridMethod = "linear”), things get messed up, but only in the NH. Note how things look okay in the SH. (See attached bottom right figure.)
My grids seem to be fine: Original longitudes:
Destination longitudes:
Any ideas?
Mark