diff --git a/README.md b/README.md
index 83f80ae..1a84652 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,12 @@
# MPAS Limited-Area - v2.1
-MPAS Limited-Area is a python tool that takes an MPAS global grid and produces
-a regional area grid given a region specifications. Regions can be specified in
-a number of different ways, making limited-area extensible and flexible.
+MPAS Limited-Area is a Python tool that takes MPAS global mesh files and
+can produce regional subsets meshes given a region specifications. The MPAS
+Limited-Area tool can create subsets MPAS mesh files that contain time variant,
+and time invariant fields, given a mesh connectivity fields are present.
+
+Regional specifications can be created in a number of ways, making limited-area
+extensible and flexible.
# Download and Installing
@@ -18,8 +22,8 @@ $ setenv PATH ${PATH}:/path/to/MPAS-Limited-Area
Then, run the `create_region` command-line script:
```Bash
$ # Running the script with no arguments will produce a usage output
-Usage: create_region [-h] [-v VERBOSE] points grid create_region
-create_region: error: the following arguments are required: points, grid
+usage: create_region [-h] [-v VERBOSE] points files [files ...]
+create_region: error: the following arguments are required: points, files
```
**Note**: It may be necessary to install the dependencies for this program if
@@ -51,20 +55,43 @@ to generate a help message, usage statement and a list of options. The command
line usage is:
``` bash
-$ create_region [options] points grid
+usage: create_region [-h] [-v VERBOSE] points files [files ...]
```
-Where `points` is a points specification file, and where `grid` is a global
-MPAS NetCDF grid file. You'll find example points file in the points-example
-directory within the docs directory of this repository and below in the *Points
-Syntax* section.
+Where `points` is a points specification file, and where `files` is any number
+of global MPAS NetCDF mesh files with the first containing the mesh
+connectivity information that will be used subset itself and the others.
+
+You can find example points file in the points-example directory within the
+docs directory of this repository and below in the *Points Syntax* section.
+
+## Example Subsetting Mesh Files
+
+MPAS-Limited-Area can subset single or multiple files:
+
+``` bash
+# Subset a mesh file with only connectivity information:
+create_region conus.circle.pts x1.10242.grid.nc
+
+# Subset a static file with time invariant fields:
+create_region conus.circle.pts x1.10242.static.nc
+
+# Subset an initial conditions files with time invariant and time variant
+# fields:
+create_region conus.circle.pts x1.10242.init.nc
+
+# Subset files that do not contain mesh connectivity (as well as the mesh that
+# contains mesh connectivity fields):
+create_region conus.circle.pts x1.10242.init.nc diag.2019-01-28_00.00.00.nc
+diag.2019-01-28_03.00.00.nc diag.2019-01-28-06.00.00.nc
+```
## Notes on creating regions from .grid.nc files
When creating region subsets from `grid.nc` files, it is important to take a
-few things into account. Because of the way the init_atmosphere model
+few things into account. Because of the way the init_atmosphere model
interpolates some static data, it is possible to create incorrect static
-feilds on a concave mesh, during the static interpolation of the
+fields on a concave mesh, during the static interpolation of the
init_atmosphere model.
Creating a mesh that is either **concave** in overall shape, or spans
@@ -115,7 +142,7 @@ appropriate value:
of the circle and ellipse. `Point` is not used with the channel
method.
-The value after `Name:` will be the name of the new regional mesh.
+The value after `Name:` will be the name of the new regional mesh.
If the `Name` within the pts file was set to be the following:
```
diff --git a/create_region b/create_region
index 789328b..079e039 100755
--- a/create_region
+++ b/create_region
@@ -5,8 +5,9 @@ import argparse
from limited_area.limited_area import LimitedArea
if __name__ == "__main__":
- description = ("Create a limited area MPAS grid from a global MPAS grid"
- " and a file that specifies the regional area boundary.")
+ description = ("Create regional subsets of global MPAS meshes using "
+ "a regional specification .pts file that specifies the "
+ "regions boundary.")
epilog = ("For more information please see: "
"https://github.com/MiCurry/MPAS-Limited-Area")
@@ -15,21 +16,23 @@ if __name__ == "__main__":
parser = argparse.ArgumentParser(description=description,
epilog=epilog)
- required = parser.add_argument_group('Required', "Limited Area requires a"
- " MPAS Grid as well as a"
- " file specifying the"
- " limited area.")
+ required_description = ("create_region requires a MPAS mesh file that "
+ "contains mesh connectivity fields and a file for "
+ "specifying a regional area (.pts file). ")
- options = parser.add_argument_group('Options', "Options to generating the"
- " MPAS limited area.")
+ required = parser.add_argument_group('Required', required_description)
+ options = parser.add_argument_group('Options')
required.add_argument('points',
help='Points file specifying the MPAS regional area',
type=str)
- required.add_argument('grid',
- help=('Global MPAS grid file. Either a grid.nc or'
- ' static.nc file'),
+ required.add_argument('files',
+ help=('Global MPAS file(s) to be subset. If multiple '
+ 'meshes are given, the first must contain mesh '
+ 'connectivity infromation that will be used to '
+ 'subset all files.'),
+ nargs='+',
type=str)
options.add_argument('-v', '--verbose',
@@ -45,13 +48,13 @@ if __name__ == "__main__":
print("DEBUG: DEBUG set to verbose level ", DEBUG, '\n')
if DEBUG > 1:
- print("DEBUG: Grid File: ", args.grid)
+ print("DEBUG: Mesh Files: ", args.files)
print("DEBUG: Points File: ", args.points)
kwargs = { 'DEBUG' : DEBUG }
- regional_area = LimitedArea(args.grid,
+ regional_area = LimitedArea(args.files,
args.points,
format='NETCDF3_64BIT_OFFSET',
**kwargs)
diff --git a/limited_area/limited_area.py b/limited_area/limited_area.py
index c6e5317..6f4f28d 100644
--- a/limited_area/limited_area.py
+++ b/limited_area/limited_area.py
@@ -17,7 +17,7 @@ class LimitedArea():
UNMARKED = 0
def __init__(self,
- mesh_file,
+ files,
region,
regionFormat='points',
format='NETCDF3_64BIT_OFFSET',
@@ -25,13 +25,19 @@ def __init__(self,
**kwargs):
""" Init function for Limited Area
- Check to see if mesh file exists and it is the correct type. Check to
- see that the region file exist and finally set the regionSpec to the
- requested regionFormat
+ Check to see if all mesh files that were passed in to files exist
+ and that they are all the correct type. Check to see if the first
+ (or only) file contains mesh connectivity. If it is, then load its
+ connectivity fields.
+ Add the mesh connectivity mesh to self.mesh, and then add all meshes
+ to the self.meshes attribute. Thus, we can use self.mesh to subset all
+ meshes in the self.meshes attribute in gen_region.
Keyword arguments:
- mesh_files -- Path to a valid MPAS Mesh file
+ files -- Path to valid MPAS mesh files. If multiple files are given, the first file
+ must contain mesh connectivity information, which will be used to subset
+ itself, and all following files.
region -- Path to pts file region specification
DEBUG -- Debug value used to turn on debug output, default == 0
@@ -39,19 +45,13 @@ def __init__(self,
is mark neighbor search
"""
+ self.meshes = []
+
# Keyword arguments
self._DEBUG_ = kwargs.get('DEBUG', 0)
self.boundary = kwargs.get('markNeighbors', 'search')
self.cdf_format = format
- # Check to see that all of the meshes exists and that they are
- # valid netCDF files.
- if os.path.isfile(mesh_file):
- self.mesh = MeshHandler(mesh_file, 'r', *args, **kwargs)
- else:
- print("ERROR: Mesh file was not found", mesh_file)
- sys.exit(-1)
-
# Check to see the points file exists and if it exists, then parse it
# and see that is is specified correctly!
self.region_file = region
@@ -64,10 +64,40 @@ def __init__(self,
elif self.boundary == 'search':
# Possibly faster for smaller regions
self.mark_neighbors = self._mark_neighbors_search
-
-
+
+ # Check to see that all given mesh files, simply exist
+ for mesh in files:
+ if not os.path.isfile(mesh):
+ print("ERROR: Mesh file was not found", mesh_file)
+ sys.exit(-1)
+
+ if len(files) == 1:
+ self.mesh = MeshHandler(files[0], 'r', *args, **kwargs)
+ if self.mesh.check_grid():
+ self.mesh.load_vars()
+ else:
+ print("ERROR:", self.mesh.fname, "did not contain needed mesh connectivity information")
+ sys.exit(-1)
+ self.meshes.append(self.mesh)
+
+ else:
+ self.mesh = MeshHandler(files.pop(0), 'r', *args, **kwargs)
+ if self.mesh.check_grid():
+ # Load the mesh connectivity variables needed to subset this mesh, and all the
+ # other ones
+ self.mesh.load_vars()
+ self.meshes.append(self.mesh)
+ else:
+ print("ERROR:", self.mesh.fname, "did not contain needed mesh connectivity information")
+ print("ERROR: The first file must contain mesh connectivity information")
+ sys.exit(-1)
+
+ for mesh in files:
+ self.meshes.append(MeshHandler(mesh, 'r', *args, **kwargs))
+
def gen_region(self, *args, **kwargs):
- """ Generate the boundary region of the given region for the given mesh(es). """
+ """ Generate the boundary region of the specified region
+ and subset meshes in self.meshes """
# Call the regionSpec to generate `name, in_point, boundaries`
name, inPoint, boundaries= self.regionSpec.gen_spec(self.region_file, **kwargs)
@@ -80,6 +110,7 @@ def gen_region(self, *args, **kwargs):
# For each mesh, create a regional mesh and save it
print('\n')
+ # TODO: Update this print statement to be more consistent to what is happening
print('Creating a regional mesh of ', self.mesh.fname)
# Mark boundaries
@@ -137,30 +168,40 @@ def gen_region(self, *args, **kwargs):
**kwargs)
- # Subset the grid into a new region:
- print('Subsetting mesh fields into the specified region mesh...')
- regionFname = self.create_regional_fname(name, self.mesh)
- regionalMesh = self.mesh.subset_fields(regionFname,
- bdyMaskCell,
- bdyMaskEdge,
- bdyMaskVertex,
- inside=self.INSIDE,
- unmarked=self.UNMARKED,
- format=self.cdf_format,
- *args,
- **kwargs)
-
- print('Copying global attributes...')
- self.mesh.copy_global_attributes(regionalMesh)
+ # Create subsets of all the meshes in self.meshes
+ print('Subsetting meshes...')
+ for mesh in self.meshes:
+ print("\nSubsetting:", mesh.fname)
+
+ regionFname = self.create_regional_fname(name, mesh.fname)
+ regionalMesh = mesh.subset_fields(regionFname,
+ bdyMaskCell,
+ bdyMaskEdge,
+ bdyMaskVertex,
+ self.INSIDE,
+ self.UNMARKED,
+ self.mesh,
+ format=self.cdf_format,
+ *args,
+ **kwargs)
+ print("Copying global attributes... ")
+ regionalMesh.copy_global_attributes(self.mesh)
+ print("Create a regional mesh:", regionFname)
+
+ if mesh.check_grid():
+ # Save the regional mesh that contains graph connectivity to create the regional
+ # graph partition file below
+ regionalMeshConn = regionalMesh
+ else:
+ regionalMesh.mesh.close()
+ mesh.mesh.close()
- print("Created a regional mesh: ", regionFname)
print('Creating graph partition file...', end=' '); sys.stdout.flush()
- graphFname = regionalMesh.create_graph_file(self.create_partiton_fname(name, self.mesh,))
+ graphFname = regionalMeshConn.create_graph_file(self.create_partiton_fname(name, self.mesh,))
print(graphFname)
- self.mesh.mesh.close()
- regionalMesh.mesh.close()
+ regionalMeshConn.mesh.close()
return regionFname, graphFname
@@ -169,26 +210,10 @@ def create_partiton_fname(self, name, mesh, **kwargs):
return name+'.graph.info'
- def create_regional_fname(self, name, mesh, **kwargs):
- """ Generate the filename for the regional mesh file. Depending on what "type" of MPAS file
- we are using, either static, grid or init, try to rename the region file as that type (i.e.
- x1.2562.static.nc becomes name.static.nc).
-
- If a file name is ambiguous, or the file name does not contain: static, init, or grid,
- rename the region file to be region. """
- # Static files
- if 'static' in mesh.fname and not ('grid' in mesh.fname or 'init' in mesh.fname):
- meshType = 'static'
- # Grid files
- elif 'grid' in mesh.fname and not ('static' in mesh.fname or 'init' in mesh.fname):
- meshType = 'grid'
- # Initialization Data
- elif 'init' in mesh.fname and not ('static' in mesh.fname or 'grid' in mesh.fname):
- meshType = 'init'
- else:
- meshType = 'region'
-
- return name+'.'+meshType+'.nc'
+ def create_regional_fname(self, regionName, meshFileName, **kwargs):
+ """ Create the regional file name by prepending the regional name
+ (specified by Name: ) in the points file, to the meshFileName. """
+ return regionName+'.'+meshFileName
# Mark_neighbors_search - Faster for smaller regions ??
diff --git a/limited_area/mesh.py b/limited_area/mesh.py
index 76c31ee..b7e0512 100644
--- a/limited_area/mesh.py
+++ b/limited_area/mesh.py
@@ -11,6 +11,7 @@
class MeshHandler:
""" Handle the operations related to NetCDF/MPAS grids. """
+
def __init__(self, fname, mode, format='NETCDF3_64BIT_OFFSET', *args, **kwargs):
""" Open fname with mode, for either reading, or creating
@@ -28,8 +29,7 @@ def __init__(self, fname, mode, format='NETCDF3_64BIT_OFFSET', *args, **kwargs):
self.fname = fname
if mode == 'r':
- if self.check_file(fname):
- self._load_vars()
+ if self.load_file(fname):
return
else:
sys.exit(-1)
@@ -48,15 +48,12 @@ def create_file(self, fname, mode, format):
sys.exit(-1)
- def check_file(self, fname):
- """ Check to see that fname exists and it is a valid NetCDF file """
+ def load_file(self, fname):
+ """ Load fname using the netCDF4 Dataset class. If fname is not a valid NetCDF file, the
+ program will exit. """
if os.path.isfile(fname):
try:
- mesh = open(fname, 'rb')
- nc_bytes = mesh.read()
- mesh.close()
-
- self.mesh = Dataset(fname, 'r', memory=nc_bytes)
+ self.mesh = Dataset(fname, 'r')
return True
except OSError as E:
print("ERROR: ", E)
@@ -66,7 +63,36 @@ def check_file(self, fname):
print("ERROR: This file did not exist!")
return False
- def _load_vars(self):
+ def check_grid(self):
+ """ Check to see if this mesh contains all of the needed mesh connectivity infromation for
+ subsetting. If it does, return True, and if not, return False. """
+ mesh_dims = ['nCells',
+ 'nEdges',
+ 'maxEdges',
+ 'nVertices',
+ 'vertexDegree']
+
+ mesh_vars = ['latCell',
+ 'lonCell',
+ 'nEdgesOnCell',
+ 'cellsOnCell',
+ 'cellsOnEdge',
+ 'cellsOnVertex',
+ 'indexToCellID',
+ 'indexToEdgeID',
+ 'indexToVertexID',]
+
+ for dim in mesh_dims:
+ if not dim in self.mesh.dimensions:
+ return False
+
+ for var in mesh_vars:
+ if not var in self.mesh.variables:
+ return False
+
+ return True
+
+ def load_vars(self):
""" Pre-load variables to avoid multiple, unnecessary IO calls
Pulling variables from a netCDF4 interface like the following:
@@ -75,7 +101,7 @@ def _load_vars(self):
I/O calls.
"""
if self._DEBUG_ > 2:
- print("DEBUG: In Load Vars")
+ print("DEBUG:", self.fname, "in Load Vars")
# Dimensions
self.nCells = self.mesh.dimensions['nCells'].size
@@ -100,16 +126,16 @@ def _load_vars(self):
# Attributes
self.sphere_radius = self.mesh.sphere_radius
- self.variables = { 'latCells' : self.latCells,
- 'lonCells' : self.lonCells,
- 'nEdgesOnCell' : self.nEdgesOnCell,
- 'cellsOnCell' : self.cellsOnCell,
- 'cellsOnEdge' : self.cellsOnEdge,
- 'cellsOnVertex' : self.cellsOnVertex,
- 'indexToCellID' : self.indexToCellIDs,
- 'indexToEdgeID' : self.indexToEdgeIDs,
+ self.variables = { 'latCells' : self.latCells,
+ 'lonCells' : self.lonCells,
+ 'nEdgesOnCell' : self.nEdgesOnCell,
+ 'cellsOnCell' : self.cellsOnCell,
+ 'cellsOnEdge' : self.cellsOnEdge,
+ 'cellsOnVertex' : self.cellsOnVertex,
+ 'indexToCellID' : self.indexToCellIDs,
+ 'indexToEdgeID' : self.indexToEdgeIDs,
'indexToVertexID' : self.indexToVertexIDs
- }
+ }
def nearest_cell(self, lat, lon):
@@ -194,10 +220,11 @@ def subset_fields(self,
bdyMaskVertex,
inside,
unmarked,
+ mesh,
format='NETCDF3_64BIT_OFFSET',
*args,
**kwargs):
- """ Subset the current mesh and return a new regional mesh with
+ """ Subset the mesh in self.mesh using mesh and return a new regional mesh with
subsetted fields
regionalFname -- Desired filename for the regional subset
@@ -210,6 +237,7 @@ def subset_fields(self,
unmarked -- The integer value that was used to mark cells,
edges, vertices as being 'outside' of the regional
mesh.
+ mesh -- The mesh connectivity to use to subset self.mesh
"""
# Don't pass on DEBUG to the regional mess - tone down output
@@ -228,9 +256,9 @@ def subset_fields(self,
indexingFields['edgesOnVertex'] = bdyMaskEdge
indexingFields['cellsOnVertex'] = bdyMaskCell
- glbBdyCellIDs = self.indexToCellIDs[np.where(bdyMaskCell != unmarked)] - 1
- glbBdyEdgeIDs = self.indexToEdgeIDs[np.where(bdyMaskEdge != unmarked)] - 1
- glbBdyVertexIDs = self.indexToVertexIDs[np.where(bdyMaskVertex != unmarked)] - 1
+ glbBdyCellIDs = mesh.indexToCellIDs[np.where(bdyMaskCell != unmarked)] - 1
+ glbBdyEdgeIDs = mesh.indexToEdgeIDs[np.where(bdyMaskEdge != unmarked)] - 1
+ glbBdyVertexIDs = mesh.indexToVertexIDs[np.where(bdyMaskVertex != unmarked)] - 1
if self._DEBUG_ > 0:
@@ -242,11 +270,11 @@ def subset_fields(self,
# len(bdyIndexToCellIDs) == nCells, then the specification was probably not
# specified correctly
force = False
- if len(glbBdyCellIDs) == self.nCells and not force:
+ if len(glbBdyCellIDs) == mesh.nCells and not force:
print("ERROR: The number of Cells in the specified region ",
"(", len(glbBdyCellIDs), ")")
print("ERROR: appears to be equal number of cells in the global mesh",
- "(", self.nCells, ")")
+ "(", mesh.nCells, ")")
print("ERROR: which means there was perhaps a problem in specifying the")
print("ERROR: region. Please insure your region specification is correct")
sys.exit(-1)
@@ -275,17 +303,18 @@ def subset_fields(self,
# Make boundary Mask's between 0 and the number of specified relaxation
# layers
- region.mesh.createVariable('bdyMaskCell', 'i4', ('nCells',))
- region.mesh.createVariable('bdyMaskEdge', 'i4', ('nEdges',))
- region.mesh.createVariable('bdyMaskVertex', 'i4', ('nVertices'))
+ if self.check_grid():
+ region.mesh.createVariable('bdyMaskCell', 'i4', ('nCells',))
+ region.mesh.createVariable('bdyMaskEdge', 'i4', ('nEdges',))
+ region.mesh.createVariable('bdyMaskVertex', 'i4', ('nVertices'))
- region.mesh.variables['bdyMaskCell'][:] = bdyMaskCell[bdyMaskCell != 0] - 1
- region.mesh.variables['bdyMaskEdge'][:] = bdyMaskEdge[bdyMaskEdge != 0] - 1
- region.mesh.variables['bdyMaskVertex'][:] = bdyMaskVertex[bdyMaskVertex != 0] - 1
+ region.mesh.variables['bdyMaskCell'][:] = bdyMaskCell[bdyMaskCell != 0] - 1
+ region.mesh.variables['bdyMaskEdge'][:] = bdyMaskEdge[bdyMaskEdge != 0] - 1
+ region.mesh.variables['bdyMaskVertex'][:] = bdyMaskVertex[bdyMaskVertex != 0] - 1
- scan(bdyMaskCell)
- scan(bdyMaskEdge)
- scan(bdyMaskVertex)
+ scan(bdyMaskCell)
+ scan(bdyMaskEdge)
+ scan(bdyMaskVertex)
# Variables - Create Variables
for var in self.mesh.variables:
@@ -310,8 +339,8 @@ def subset_fields(self,
continue
print("Copying variable ", var, "...", end=' ', sep=''); sys.stdout.flush()
- if var in self.variables:
- arrTemp = self.variables[var] # Use the pre-loaded variable if possible
+ if var in mesh.variables:
+ arrTemp = mesh.variables[var] # Use the pre-loaded variable if possible
else:
arrTemp = self.mesh.variables[var][:] # Else, read it from disk
@@ -355,10 +384,12 @@ def subset_fields(self,
return region
- def copy_global_attributes(self, region):
- """ Copy the global attributes into the regional mesh, but not 'np' """
- region.mesh.on_a_sphere = self.mesh.on_a_sphere
- region.mesh.sphere_radius = self.mesh.sphere_radius
+ def copy_global_attributes(self, mesh):
+ """ Copy the global attributes from mesh, onto self """
+ if self._DEBUG_ > 1:
+ print("DEBUG: Copying global attributes from", mesh.fname, "to", self.fname)
+ self.mesh.on_a_sphere = mesh.mesh.on_a_sphere
+ self.mesh.sphere_radius = mesh.mesh.sphere_radius
def scan(arr):