Add ability to subset multiple files in a single call

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<a name="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<a name="concave">
 
 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:
 ```

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)

limited_area/limited_area.py

@@ -17,41 +17,41 @@ class LimitedArea():
     UNMARKED = 0
 
     def __init__(self,
-                 mesh_file,
+                 files,
                  region,
                  regionFormat='points',
                  format='NETCDF3_64BIT_OFFSET',
                  *args,
                  **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
         markNeighbors -- Algorithm choice for choosing relaxation layers - Default
                          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 ??