Tutorials¶
The first few tutorials are stand-alone scripts that can be run from any Python interpreter.
Hello World¶
import ESMF # This call enables debug logging # esmpy = ESMF.Manager(debug=True) print "Hello ESMPy World from PET (processor) {0}!".format(ESMF.local_pet())
Create and Read from File¶
# This example demonstrates how to create ESMPy Grid, Mesh and Field objects from file. # The data files can be retrieved from the ESMF data repository by uncommenting the # following block of code: # # import os # DD = os.path.join(os.getcwd(), "examples/data") # if not os.path.isdir(DD): # os.makedirs(DD) # from ESMF.util.cache_data import cache_data_file # cache_data_file(os.path.join(DD, "so_Omon_GISS-E2.nc")) # cache_data_file(os.path.join(DD, "mpas_uniform_10242_dual_counterclockwise.nc")) import os import ESMF # This call enables debug logging # ESMF.Manager(debug=True) # set up the DATADIR DATADIR = os.path.join(os.getcwd(), "examples/data") # Create a uniform global latlon grid from a SCRIP formatted file grid = ESMF.Grid(filename=os.path.join(DATADIR, "so_Omon_GISS-E2.nc"), filetype=ESMF.FileFormat.GRIDSPEC) # Create a field on the centers of the grid field = ESMF.Field(grid, staggerloc=ESMF.StaggerLoc.CENTER)#, ndbounds=[2]) # field.read does not work if ESMF is built with MPIUNI if ESMF.api.constants._ESMF_COMM is not ESMF.api.constants._ESMF_COMM_MPIUNI: field.read(filename=os.path.join(DATADIR, "so_Omon_GISS-E2.nc"), variable="so") # create an ESMF formatted unstructured mesh with clockwise cells removed mesh = ESMF.Mesh(filename=os.path.join(DATADIR, "mpas_uniform_10242_dual_counterclockwise.nc"), filetype=ESMF.FileFormat.ESMFMESH) # create a field on the nodes of the mesh field = ESMF.Field(mesh, meshloc=ESMF.MeshLoc.NODE) if ESMF.local_pet() == 0: print "Grid, Mesh and Field all created/read from file successfully :)"
Regridding from Grid to Mesh¶
# This example demonstrates how to regrid between a Grid and a Mesh. # The data files can be retrieved from the ESMF data repository by uncommenting the # following block of code: # # import os # DD = os.path.join(os.getcwd(), "examples/data") # if not os.path.isdir(DD): # os.makedirs(DD) # from ESMF.util.cache_data import cache_data_file # cache_data_file(os.path.join(DD, "ll1deg_grid.nc")) # cache_data_file(os.path.join(DD, "mpas_uniform_10242_dual_counterclockwise.nc")) import ESMF import numpy # This call enables debug logging # esmpy = ESMF.Manager(debug=True) grid1 = "examples/data/ll1deg_grid.nc" grid2 = "examples/data/mpas_uniform_10242_dual_counterclockwise.nc" # Create a uniform global latlon grid from a SCRIP formatted file grid = ESMF.Grid(filename=grid1, filetype=ESMF.FileFormat.SCRIP) # NOTE: corners are needed for conservative regridding # grid = ESMF.Grid(filename=grid1, filetype=ESMF.FileFormat.SCRIP, # add_corner_stagger=True) # create a field on the center stagger locations of the source grid srcfield = ESMF.Field(grid, name='srcfield', staggerloc=ESMF.StaggerLoc.CENTER) # create an ESMF formatted unstructured mesh with clockwise cells removed mesh = ESMF.Mesh(filename=grid2, filetype=ESMF.FileFormat.ESMFMESH) # create a field on the nodes of the destination mesh dstfield = ESMF.Field(mesh, name='dstfield', meshloc=ESMF.MeshLoc.NODE) xctfield = ESMF.Field(mesh, name='xctfield', meshloc=ESMF.MeshLoc.NODE) # NOTE: Field must be built on elements of Mesh for conservative regridding # dstfield = ESMF.Field(mesh, name='dstfield', meshloc=ESMF.MeshLoc.ELEMENT) # xctfield = ESMF.Field(mesh, name='xctfield', meshloc=ESMF.MeshLoc.ELEMENT) # initialize the fields [lon,lat] = [0, 1] deg2rad = 3.14159/180 gridXCoord = srcfield.grid.get_coords(lon, ESMF.StaggerLoc.CENTER) gridYCoord = srcfield.grid.get_coords(lat, ESMF.StaggerLoc.CENTER) srcfield.data[...] = 10.0 + (gridXCoord * deg2rad) ** 2 + (gridYCoord * deg2rad) ** 2 gridXCoord = xctfield.grid.get_coords(lon, ESMF.StaggerLoc.CENTER) gridYCoord = xctfield.grid.get_coords(lat, ESMF.StaggerLoc.CENTER) xctfield.data[...] = 10.0 + (gridXCoord * deg2rad) ** 2 + (gridYCoord * deg2rad) ** 2 dstfield.data[...] = 1e20 # create an object to regrid data from the source to the destination field regrid = ESMF.Regrid(srcfield, dstfield, regrid_method=ESMF.RegridMethod.BILINEAR, unmapped_action=ESMF.UnmappedAction.ERROR) # do the regridding from source to destination field dstfield = regrid(srcfield, dstfield) # compute the mean relative error from operator import mul num_nodes = reduce(mul, xctfield.data.shape) relerr = 0 meanrelerr = 0 if num_nodes is not 0: ind = numpy.where((dstfield.data != 1e20) & (xctfield.data != 0))[0] relerr = numpy.sum(numpy.abs(dstfield.data[ind] - xctfield.data[ind]) / numpy.abs(xctfield.data[ind])) meanrelerr = relerr / num_nodes # handle the parallel case if ESMF.pet_count() > 1: try: from mpi4py import MPI except: raise ImportError comm = MPI.COMM_WORLD relerr = comm.reduce(relerr, op=MPI.SUM) num_nodes = comm.reduce(num_nodes, op=MPI.SUM) # output the results from one processor only if ESMF.local_pet() is 0: meanrelerr = relerr / num_nodes print "ESMPy Grid Mesh Regridding Example" print " interpolation mean relative error = {0}".format(meanrelerr) assert (meanrelerr < 2e-3)
Regridding from Grid to LocStream¶
# This example demonstrates how to regrid between a Grid and a LocStream. # The data files can be retrieved from the ESMF data repository by uncommenting the # following block of code: # # import os # DD = os.path.join(os.getcwd(), "examples/data") # if not os.path.isdir(DD): # os.makedirs(DD) # from ESMF.util.cache_data import cache_data_file # cache_data_file(os.path.join(DD, "ll1deg_grid.nc")) import ESMF import numpy # This call enables debug logging ESMF.Manager(debug=True) from ESMF.test.test_api.locstream_utilities import create_locstream_spherical_16, create_locstream_spherical_16_parallel coord_sys=ESMF.CoordSys.SPH_DEG domask=True if ESMF.pet_count() == 1: locstream = create_locstream_spherical_16(coord_sys=coord_sys, domask=domask) else: if ESMF.pet_count() is not 4: raise ValueError("processor count must be 4 or 1 for this example") else: locstream = create_locstream_spherical_16_parallel(coord_sys=coord_sys, domask=domask) grid1 = "examples/data/ll1deg_grid.nc" grid = ESMF.Grid(filename=grid1, filetype=ESMF.FileFormat.SCRIP) # create a field srcfield = ESMF.Field(grid, name='srcfield') dstfield = ESMF.Field(locstream, name='dstfield') xctfield = ESMF.Field(locstream, name='xctfield') # initialize the fields [x, y] = [0, 1] deg2rad = 3.14159/180 gridXCoord = srcfield.grid.get_coords(x) gridYCoord = srcfield.grid.get_coords(y) srcfield.data[...] = 10.0 + numpy.cos(gridXCoord * deg2rad) ** 2 + numpy.cos(2 * gridYCoord * deg2rad) gridXCoord = locstream["ESMF:Lon"] gridYCoord = locstream["ESMF:Lat"] if coord_sys == ESMF.CoordSys.SPH_DEG: xctfield.data[...] = 10.0 + numpy.cos(gridXCoord * deg2rad) ** 2 + numpy.cos(2 * gridYCoord * deg2rad) elif coord_sys == ESMF.CoordSys.SPH_RAD: xctfield.data[...] = 10.0 + numpy.cos(gridXCoord) ** 2 + numpy.cos(2 * gridYCoord) else: raise ValueError("coordsys value does not work in this example") dstfield.data[...] = 1e20 # create an object to regrid data from the source to the destination field dst_mask_values=None if domask: dst_mask_values=numpy.array([0]) regrid = ESMF.Regrid(srcfield, dstfield, regrid_method=ESMF.RegridMethod.BILINEAR, unmapped_action=ESMF.UnmappedAction.ERROR, dst_mask_values=dst_mask_values) # do the regridding from source to destination field dstfield = regrid(srcfield, dstfield, zero_region=ESMF.Region.SELECT) # compute the mean relative error from operator import mul num_nodes = reduce(mul, xctfield.data.shape) relerr = 0 meanrelerr = 0 dstfield = numpy.ravel(dstfield.data) xctfield = numpy.ravel(xctfield.data) if num_nodes is not 0: ind = numpy.where((dstfield != 1e20) & (xctfield != 0))[0] relerr = numpy.sum(numpy.abs(dstfield[ind] - xctfield[ind]) / numpy.abs(xctfield[ind])) meanrelerr = relerr / num_nodes # handle the parallel case if ESMF.pet_count() > 1: try: from mpi4py import MPI except: raise ImportError comm = MPI.COMM_WORLD relerr = comm.reduce(relerr, op=MPI.SUM) num_nodes = comm.reduce(num_nodes, op=MPI.SUM) # output the results from one processor only if ESMF.local_pet() is 0: meanrelerr = relerr / num_nodes print "ESMPy Grid LocStream Regridding Example" print " interpolation mean relative error = {0}".format(meanrelerr) assert (meanrelerr < 2e-2)
Regridding from Mesh to LocStream¶
# This example demonstrates how to regrid between a mesh and a locstream. import ESMF import numpy # This call enables debug logging # ESMF.Manager(debug=True) from ESMF.test.test_api.mesh_utilities import mesh_create_5, mesh_create_5_parallel from ESMF.test.test_api.locstream_utilities import create_locstream_16, create_locstream_16_parallel if ESMF.pet_count() == 1: mesh, _, _, _, _, _ = mesh_create_5() locstream = create_locstream_16() else: if ESMF.pet_count() is not 4: raise ValueError("processor count must be 4 or 1 for this example") else: mesh, _, _, _, _ = mesh_create_5_parallel() locstream = create_locstream_16_parallel() # create a field srcfield = ESMF.Field(mesh, name='srcfield')#, meshloc=ESMF.MeshLoc.ELEMENT) # create a field on the locstream dstfield = ESMF.Field(locstream, name='dstfield') xctfield = ESMF.Field(locstream, name='xctfield') # initialize the fields [x, y] = [0, 1] deg2rad = 3.14159/180 gridXCoord = srcfield.grid.get_coords(x) gridYCoord = srcfield.grid.get_coords(y) srcfield.data[...] = 10.0 + (gridXCoord * deg2rad) ** 2 + (gridYCoord * deg2rad) ** 2 gridXCoord = locstream["ESMF:X"] gridYCoord = locstream["ESMF:Y"] xctfield.data[...] = 10.0 + (gridXCoord * deg2rad) ** 2 + (gridYCoord * deg2rad) ** 2 dstfield.data[...] = 1e20 # create an object to regrid data from the source to the destination field # TODO: this example seems to fail occasionally with UnmappedAction.ERROR, probably due to a tolerance issue - ask Bob regrid = ESMF.Regrid(srcfield, dstfield, regrid_method=ESMF.RegridMethod.BILINEAR, unmapped_action=ESMF.UnmappedAction.IGNORE) # do the regridding from source to destination field dstfield = regrid(srcfield, dstfield) # compute the mean relative error from operator import mul num_nodes = reduce(mul, xctfield.data.shape) relerr = 0 meanrelerr = 0 if num_nodes is not 0: ind = numpy.where((dstfield.data != 1e20) & (xctfield.data != 0))[0] relerr = numpy.sum(numpy.abs(dstfield.data[ind] - xctfield.data[ind]) / numpy.abs(xctfield.data[ind])) meanrelerr = relerr / num_nodes # handle the parallel case if ESMF.pet_count() > 1: try: from mpi4py import MPI except: raise ImportError comm = MPI.COMM_WORLD relerr = comm.reduce(relerr, op=MPI.SUM) num_nodes = comm.reduce(num_nodes, op=MPI.SUM) # output the results from one processor only if ESMF.local_pet() is 0: meanrelerr = relerr / num_nodes print "ESMPy Grid Mesh Regridding Example" print " interpolation mean relative error = {0}".format(meanrelerr) assert (meanrelerr < 3e-5)
Regridding from LocStream to Grid¶
# This example demonstrates how to regrid between a LocStream and a Grid. # The data files can be retrieved from the ESMF data repository by uncommenting the # following block of code: # # import os # DD = os.path.join(os.getcwd(), "examples/data") # if not os.path.isdir(DD): # os.makedirs(DD) # from ESMF.util.cache_data import cache_data_file # cache_data_file(os.path.join(DD, "ll1deg_grid.nc")) import ESMF import numpy # This call enables debug logging ESMF.Manager(debug=True) grid1 = "examples/data/ll1deg_grid.nc" grid = ESMF.Grid(filename=grid1, filetype=ESMF.FileFormat.SCRIP) from ESMF.test.test_api.locstream_utilities import create_locstream_spherical_16, create_locstream_spherical_16_parallel coord_sys=ESMF.CoordSys.SPH_DEG domask=True if ESMF.pet_count() == 1: locstream = create_locstream_spherical_16(coord_sys=coord_sys, domask=domask) else: if ESMF.pet_count() is not 4: raise ValueError("processor count must be 4 or 1 for this example") else: locstream = create_locstream_spherical_16_parallel(coord_sys=coord_sys, domask=domask) # create a field srcfield = ESMF.Field(locstream, name='srcfield') dstfield = ESMF.Field(grid, name='dstfield') xctfield = ESMF.Field(grid, name='xctfield') # initialize the fields [x, y] = [0, 1] deg2rad = 3.14159/180 gridXCoord = locstream["ESMF:Lon"] gridYCoord = locstream["ESMF:Lat"] if coord_sys == ESMF.CoordSys.SPH_DEG: srcfield.data[...] = 10.0 + numpy.cos(gridXCoord * deg2rad) ** 2 + numpy.cos(2 * gridYCoord * deg2rad) elif coord_sys == ESMF.CoordSys.SPH_RAD: srcfield.data[...] = 10.0 + numpy.cos(gridXCoord) ** 2 + numpy.cos(2 * gridYCoord) else: raise ValueError("coordsys value does not apply in this example") gridXCoord = xctfield.grid.get_coords(x) gridYCoord = xctfield.grid.get_coords(y) xctfield.data[...] = 10.0 + numpy.cos(gridXCoord * deg2rad) ** 2 + numpy.cos(2 * gridYCoord * deg2rad) dstfield.data[...] = 1e20 # create an object to regrid data from the source to the destination field mask_values=None if domask: mask_values=numpy.array([0]) regrid = ESMF.Regrid(srcfield, dstfield, regrid_method=ESMF.RegridMethod.NEAREST_DTOS, unmapped_action=ESMF.UnmappedAction.ERROR, src_mask_values=mask_values) # do the regridding from source to destination field dstfield = regrid(srcfield, dstfield, zero_region=ESMF.Region.SELECT) # compute the mean relative error from operator import mul num_nodes = reduce(mul, xctfield.data.shape) relerr = 0 meanrelerr = 0 dstfield = numpy.ravel(dstfield.data) xctfield = numpy.ravel(xctfield.data) if num_nodes is not 0: ind = numpy.where((dstfield != 1e20) & (xctfield != 0))[0] relerr = numpy.sum(numpy.abs(dstfield[ind] - xctfield[ind]) / numpy.abs(xctfield[ind])) meanrelerr = relerr / num_nodes # handle the parallel case if ESMF.pet_count() > 1: try: from mpi4py import MPI except: raise ImportError comm = MPI.COMM_WORLD relerr = comm.reduce(relerr, op=MPI.SUM) num_nodes = comm.reduce(num_nodes, op=MPI.SUM) # output the results from one processor only if ESMF.local_pet() is 0: meanrelerr = relerr / num_nodes print "ESMPy LocStream Grid Regridding Example" print " interpolation mean relative error = {0}".format(meanrelerr) assert (meanrelerr < 9e-5)
Regridding Helper Functions¶
The following code snippets demonstrate how to build all of the pieces necessary to regrid data between Fields built on Grids, Meshes and LocStreams.
LocStream Create¶
def create_locstream_spherical_16(coord_sys=ESMF.CoordSys.SPH_DEG, domask=False):
"""
:param coord_sys: the coordinate system of the LocStream
:param domask: a boolean to tell whether or not to add a mask
:return: LocStream
"""
if ESMF.pet_count() is not 1:
raise ValueError("processor count must be 1 to use this function")
locstream = ESMF.LocStream(16, coord_sys=coord_sys)
deg_rad = pi
if coord_sys == ESMF.CoordSys.SPH_DEG:
deg_rad = 180
locstream["ESMF:Lon"] = [0.0, 0.5*deg_rad, 1.5*deg_rad, 2*deg_rad, 0.0, 0.5*deg_rad, 1.5*deg_rad, 2*deg_rad, 0.0, 0.5*deg_rad, 1.5*deg_rad, 2*deg_rad, 0.0, 0.5*deg_rad, 1.5*deg_rad, 2*deg_rad]
locstream["ESMF:Lat"] = [deg_rad/-2.0, deg_rad/-2.0, deg_rad/-2.0, deg_rad/-2.0, -0.25*deg_rad, -0.25*deg_rad, -0.25*deg_rad, -0.25*deg_rad, 0.25*deg_rad, 0.25*deg_rad, 0.25*deg_rad, 0.25*deg_rad, deg_rad/2.0, deg_rad/2.0, deg_rad/2.0, deg_rad/2.0]
if domask:
locstream["ESMF:Mask"] = np.array([1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=np.int32)
return locstream
LocStream Create Parallel¶
def create_locstream_spherical_16_parallel(coord_sys=ESMF.CoordSys.SPH_DEG, domask=False):
"""
:param coord_sys: the coordinate system of the LocStream
:param domask: a boolean to tell whether or not to add a mask
:return: LocStream
"""
if ESMF.pet_count() is not 4:
raise ValueError("processor count must be 4 to use this function")
deg_rad = pi
if coord_sys == ESMF.CoordSys.SPH_DEG:
deg_rad = 180.0
locstream = None
if ESMF.local_pet() is 0:
locstream = ESMF.LocStream(4, coord_sys=coord_sys)
locstream["ESMF:Lon"] = [0.0, 0.5*deg_rad, 0.0, 0.5*deg_rad]
locstream["ESMF:Lat"] = [deg_rad/-2.0, deg_rad/-2.0, -0.25*deg_rad, -0.25*deg_rad]
if domask:
locstream["ESMF:Mask"] = np.array([1, 0, 1, 1], dtype=np.int32)
elif ESMF.local_pet() is 1:
locstream = ESMF.LocStream(4, coord_sys=coord_sys)
locstream["ESMF:Lon"] = [1.5*deg_rad, 2*deg_rad, 1.5*deg_rad, 2*deg_rad]
locstream["ESMF:Lat"] = [deg_rad/-2.0, deg_rad/-2.0, -0.25*deg_rad, -0.25*deg_rad]
if domask:
locstream["ESMF:Mask"] = np.array([0, 1, 1, 1], dtype=np.int32)
elif ESMF.local_pet() is 2:
locstream = ESMF.LocStream(4, coord_sys=coord_sys)
locstream["ESMF:Lon"] = [0.0, 0.5*deg_rad, 0.0, 0.5*deg_rad]
locstream["ESMF:Lat"] = [0.25*deg_rad, 0.25*deg_rad, deg_rad/2.0, deg_rad/2.0]
if domask:
locstream["ESMF:Mask"] = np.array([1, 1, 1, 1], dtype=np.int32)
elif ESMF.local_pet() is 3:
locstream = ESMF.LocStream(4, coord_sys=coord_sys)
locstream["ESMF:Lon"] = [1.5*deg_rad, 2*deg_rad, 1.5*deg_rad, 2*deg_rad]
locstream["ESMF:Lat"] = [0.25*deg_rad, 0.25*deg_rad, deg_rad/2.0, deg_rad/2.0]
if domask:
locstream["ESMF:Mask"] = np.array([1, 1, 1, 1], dtype=np.int32)
return locstream
Create a 2D Grid¶
def grid_create(xdom, ydom, nx, ny, corners=False, domask=False, doarea=False, ctk=ESMF.TypeKind.R8): """ :param xdom: 2,1 list containing the x domain :param ydom: 2,1 list conatining the y domain :param nx: number of longitude values at cell centers :param ny: number of latitude values at cell centers :param corners: boolean to determine whether or not to add corner coordinates to this grid :param domask: boolean to determine whether to set an arbitrary mask or not :param doarea: boolean to determine whether to set an arbitrary area values or not :param ctk: the coordinate typekind :return: grid """ [x, y] = [0, 1] # Create arrays of center and corner values to emulate what would be read from a standard CF-like file # +1 because corners have one more than center xs = np.linspace(xdom[0], xdom[1], nx + 1) xcorner = np.array([xs[0:-1], xs[1::]]).T xcenter = (xcorner[:, 1] + xcorner[:, 0]) / 2 # +1 because corners have one more than center ys = np.linspace(ydom[0], ydom[1], ny + 1) ycorner = np.array([ys[0:-1], ys[1::]]).T ycenter = (ycorner[:, 1] + ycorner[:, 0]) / 2 # create a grid given the number of grid cells in each dimension, the center stagger location is allocated, the # Cartesian coordinate system and type of the coordinates are specified max_index = np.array([nx, ny]) grid = ESMF.Grid(max_index, staggerloc=[ESMF.StaggerLoc.CENTER], coord_sys=ESMF.CoordSys.CART, coord_typekind=ctk) # set the grid coordinates using numpy arrays, parallel case is handled using grid bounds gridXCenter = grid.get_coords(x) x_par = xcenter[grid.lower_bounds[ESMF.StaggerLoc.CENTER][x]:grid.upper_bounds[ESMF.StaggerLoc.CENTER][x]] gridXCenter[...] = x_par.reshape((x_par.size, 1)) gridYCenter = grid.get_coords(y) y_par = ycenter[grid.lower_bounds[ESMF.StaggerLoc.CENTER][y]:grid.upper_bounds[ESMF.StaggerLoc.CENTER][y]] gridYCenter[...] = y_par.reshape((1, y_par.size)) # create grid corners in a slightly different manner to account for the bounds format common in CF-like files if corners: grid.add_coords([ESMF.StaggerLoc.CORNER]) lbx = grid.lower_bounds[ESMF.StaggerLoc.CORNER][x] ubx = grid.upper_bounds[ESMF.StaggerLoc.CORNER][x] lby = grid.lower_bounds[ESMF.StaggerLoc.CORNER][y] uby = grid.upper_bounds[ESMF.StaggerLoc.CORNER][y] gridXCorner = grid.get_coords(x, staggerloc=ESMF.StaggerLoc.CORNER) for i0 in range(ubx - lbx - 1): gridXCorner[i0, :] = xcorner[i0+lbx, 0] gridXCorner[i0 + 1, :] = xcorner[i0+lbx, 1] gridYCorner = grid.get_coords(y, staggerloc=ESMF.StaggerLoc.CORNER) for i1 in range(uby - lby - 1): gridYCorner[:, i1] = ycorner[i1+lby, 0] gridYCorner[:, i1 + 1] = ycorner[i1+lby, 1] # add an arbitrary mask if domask: mask = grid.add_item(ESMF.GridItem.MASK) mask[:] = 1 mask[np.where((1.75 <= gridXCenter.data < 2.25) & (1.75 <= gridYCenter.data < 2.25))] = 0 # add arbitrary areas values if doarea: area = grid.add_item(ESMF.GridItem.AREA) area[:] = 5.0 return grid
Create a 3D Grid¶
def grid_create_3d(xdom, ydom, zdom, nx, ny, nz, corners=False, domask=False, doarea=False): """ :param xdom: 2,1 list containing the x domain :param ydom: 2,1 list conatining the y domain :param zdom: 2,1 list conatining the z domain :param nx: number of x values at cell centers :param ny: number of y values at cell centers :param nz: number of z values at cell centers :param corners: boolean to determine whether or not to add corner coordinates to this grid :param domask: boolean to determine whether to set an arbitrary mask or not :param doarea: boolean to determine whether to set an arbitrary area values or not :return: grid """ [x, y, z] = [0, 1, 2] # Create arrays of center and corner values to emulate what would be read from a standard CF-like file # +1 because corners have one more than center xs = np.linspace(xdom[0], xdom[1], nx + 1) xcorner = np.array([xs[0:-1], xs[1::]]).T xcenter = (xcorner[:, 1] + xcorner[:, 0]) / 2 # +1 because corners have one more than center ys = np.linspace(ydom[0], ydom[1], ny + 1) ycorner = np.array([ys[0:-1], ys[1::]]).T ycenter = (ycorner[:, 1] + ycorner[:, 0]) / 2 # +1 because corners have one more than center zs = np.linspace(zdom[0], zdom[1], nz + 1) zcorner = np.array([zs[0:-1], zs[1::]]).T zcenter = (zcorner[:, 1] + zcorner[:, 0]) / 2 # create a grid given the number of grid cells in each dimension, the center stagger location is allocated and the # Cartesian coordinate system is specified max_index = np.array([nx, ny, nz]) grid = ESMF.Grid(max_index, staggerloc=[ESMF.StaggerLoc.CENTER_VCENTER], coord_sys=ESMF.CoordSys.CART) # set the grid coordinates using numpy arrays, parallel case is handled using grid bounds gridXCenter = grid.get_coords(x) x_par = xcenter[grid.lower_bounds[ESMF.StaggerLoc.CENTER_VCENTER][x]:grid.upper_bounds[ESMF.StaggerLoc.CENTER_VCENTER][x]] gridXCenter[...] = x_par.reshape(x_par.size, 1, 1) gridYCenter = grid.get_coords(y) y_par = ycenter[grid.lower_bounds[ESMF.StaggerLoc.CENTER_VCENTER][y]:grid.upper_bounds[ESMF.StaggerLoc.CENTER_VCENTER][y]] gridYCenter[...] = y_par.reshape(1, y_par.size, 1) gridZCenter = grid.get_coords(z) z_par = zcenter[grid.lower_bounds[ESMF.StaggerLoc.CENTER_VCENTER][z]:grid.upper_bounds[ESMF.StaggerLoc.CENTER_VCENTER][z]] gridZCenter[...] = z_par.reshape(1, 1, z_par.size) # create grid corners in a slightly different manner to account for the bounds format common in CF-like files if corners: grid.add_coords([ESMF.StaggerLoc.CORNER_VFACE]) lbx = grid.lower_bounds[ESMF.StaggerLoc.CORNER_VFACE][x] ubx = grid.upper_bounds[ESMF.StaggerLoc.CORNER_VFACE][x] lby = grid.lower_bounds[ESMF.StaggerLoc.CORNER_VFACE][y] uby = grid.upper_bounds[ESMF.StaggerLoc.CORNER_VFACE][y] lbz = grid.lower_bounds[ESMF.StaggerLoc.CORNER_VFACE][z] ubz = grid.upper_bounds[ESMF.StaggerLoc.CORNER_VFACE][z] gridXCorner = grid.get_coords(x, staggerloc=ESMF.StaggerLoc.CORNER_VFACE) for i0 in range(ubx - lbx - 1): gridXCorner[i0, :, :] = xcorner[i0+lbx, 0] gridXCorner[i0 + 1, :, :] = xcorner[i0+lbx, 1] gridYCorner = grid.get_coords(y, staggerloc=ESMF.StaggerLoc.CORNER_VFACE) for i1 in range(uby - lby - 1): gridYCorner[:, i1, :] = ycorner[i1+lby, 0] gridYCorner[:, i1 + 1, :] = ycorner[i1+lby, 1] gridZCorner = grid.get_coords(z, staggerloc=ESMF.StaggerLoc.CORNER_VFACE) for i2 in range(ubz - lbz - 1): gridZCorner[:, :, i2] = zcorner[i2+lbz, 0] gridZCorner[:, :, i2 + 1] = zcorner[i2+lbz, 1] # add an arbitrary mask if domask: mask = grid.add_item(ESMF.GridItem.MASK) mask[:] = 1 mask[np.where((1.75 < gridXCenter.data < 2.25) & (1.75 < gridYCenter.data < 2.25) & (1.75 < gridZCenter.data < 2.25))] = 0 # add arbitrary areas values if doarea: area = grid.add_item(ESMF.GridItem.AREA) area[:] = 5.0 return grid
Create a periodic Grid¶
def grid_create_periodic(nlon, nlat, corners=False, domask=False): """ :param nlons: number of longitude values at cell centers :param nlats: number of latitude values at cell centers :param corners: boolean to determine whether or not to add corner coordinates to this grid :param domask: boolean to determine whether to set an arbitrary mask or not :return: grid """ [lon, lat] = [0, 1] # Create arrays of center and corner values to emulate what would be read from a standard CF-like file # +1 because corners have one more than center lons = np.linspace(-180, 180, nlon + 1) loncorner = np.array([lons[0:-1], lons[1::]]).T loncenter = (loncorner[:, 1] + loncorner[:, 0]) / 2 # +1 because corners have one more than center lats = np.linspace(-90, 90, nlat + 1) latcorner = np.array([lats[0:-1], lats[1::]]).T latcenter = (latcorner[:, 1] + latcorner[:, 0]) / 2 # create a grid given the number of grid cells in each dimension the center stagger location is allocated max_index = np.array([nlon, nlat]) grid = ESMF.Grid(max_index, num_peri_dims=1, staggerloc=[ESMF.StaggerLoc.CENTER]) # set the grid coordinates using numpy arrays, parallel case is handled using grid bounds gridXCenter = grid.get_coords(lon) lon_par = loncenter[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 = latcenter[grid.lower_bounds[ESMF.StaggerLoc.CENTER][lat]:grid.upper_bounds[ESMF.StaggerLoc.CENTER][lat]] gridYCenter[...] = lat_par.reshape((1, lat_par.size)) # create grid corners in a slightly different manner to account for the bounds format common in CF-like files if corners: grid.add_coords([ESMF.StaggerLoc.CORNER]) 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, :] = loncorner[i0+lbx, 0] gridXCorner[i0 + 1, :] = loncorner[i0+lbx, 1] gridYCorner = grid.get_coords(lat, staggerloc=ESMF.StaggerLoc.CORNER) for i1 in range(uby - lby - 1): gridYCorner[:, i1] = latcorner[i1+lby, 0] gridYCorner[:, i1 + 1] = latcorner[i1+lby, 1] # add an arbitrary mask if domask: mask = grid.add_item(ESMF.GridItem.MASK) mask[:] = 1 mask[np.where((1.75 <= gridXCenter.data < 2.25) & (1.75 <= gridYCenter.data < 2.25))] = 0 return grid
Create a 5 element Mesh¶
def mesh_create_5(): ''' PRECONDITIONS: None POSTCONDITIONS: A 5 element Mesh has been created. RETURN VALUES: \n Mesh :: mesh \n 4.0 31 ------ 32 ------ 33 | | 22 / | | 21 | / | | | / 23 | 2.0 21 ------ 22 ------ 23 | | | | 11 | 12 | | | | 0.0 11 ------ 12 ------ 13 0.0 2.0 4.0 Node Ids at corners Element Ids in centers Note: This mesh is not parallel, it can only be used in serial ''' # Two parametric dimensions, and two spatial dimensions mesh = ESMF.Mesh(parametric_dim=2, spatial_dim=2) num_node = 9 num_elem = 5 nodeId = np.array([11,12,13,21,22,23,31,32,33]) nodeCoord = np.array([0.0,0.0, # node 11 2.0,0.0, # node 12 4.0,0.0, # node 13 0.0,2.0, # node 21 2.0,2.0, # node 22 4.0,2.0, # node 23 0.0,4.0, # node 31 2.0,4.0, # node 32 4.0,4.0]) # node 33 nodeOwner = np.zeros(num_node) elemId = np.array([11,12,21,22,23]) elemType=np.array([ESMF.MeshElemType.QUAD, ESMF.MeshElemType.QUAD, ESMF.MeshElemType.QUAD, ESMF.MeshElemType.TRI, ESMF.MeshElemType.TRI]) elemConn=np.array([0,1,4,3, # element 11 1,2,5,4, # element 12 3,4,7,6, # element 21 4,8,7, # element 22 4,5,8]) # element 23 elemCoord = np.array([1.0, 1.0, 3.0, 1.0, 1.0, 3.0, 2.5, 3.5, 3.5, 2.5]) mesh.add_nodes(num_node,nodeId,nodeCoord,nodeOwner) mesh.add_elements(num_elem,elemId,elemType,elemConn, element_coords=elemCoord) return mesh, nodeCoord, nodeOwner, elemType, elemConn, elemCoord
Create a Field¶
Initialize an analytic Field¶
def initialize_field_grid_periodic(field): ''' PRECONDITIONS: A Field has been created as 'field' with a 'grid' where coordinates have been set on both the center and corner stagger locations. \n POSTCONDITIONS: The 'field' has been initialized to an analytic field.\n RETURN VALUES: \n Field :: field \n ''' DEG2RAD = 3.141592653589793/180.0 # get the coordinate pointers and set the coordinates [x,y] = [0, 1] gridXCoord = field.grid.get_coords(x, ESMF.StaggerLoc.CENTER) gridYCoord = field.grid.get_coords(y, ESMF.StaggerLoc.CENTER) field.data[:] = 2.0 + np.cos(DEG2RAD*gridXCoord)**2 * \ np.cos(2.0*DEG2RAD*(90.0 - gridYCoord)) return field
Run ESMPy regridding¶
Compute Field mass¶
def compute_mass_grid(valuefield, dofrac=False, fracfield=None, uninitval=422397696.): ''' PRECONDITIONS: 'fracfield' contains the fractions of each cell which contributed to a regridding operation involving 'valuefield. 'dofrac' is a boolean value that gives the option to not use the 'fracfield'.\n POSTCONDITIONS: The mass of the data field is computed.\n RETURN VALUES: float :: mass \n ''' mass = 0.0 areafield = ESMF.Field(valuefield.grid, name='areafield') areafield.get_area() ind = np.where(valuefield.data != uninitval) if dofrac: mass = np.sum(areafield.data[ind] * valuefield.data[ind] * fracfield.data[ind]) else: mass = np.sum(areafield.data[ind] * valuefield.data[ind]) return mass