-
Notifications
You must be signed in to change notification settings - Fork 68
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
ESMF regrid problem #1125
Comments
Can you also send us the latitudes? Or point me to the file. Thanks From: "Paul J. Durack" <[email protected]mailto:[email protected]> From: , Mark Hi Charles and Paul, Have you guys ever encountered the following issue? If I regrid using But when I regrid using My grids seem to be fine: SWCRE.getLongitude()[:] Destination longitudes: horiz_grid.getLongitude()[:] Any ideas? Mark — |
Destination latitudes: 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:
|
@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
|
Ryan scripts: # 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
The text was updated successfully, but these errors were encountered: