Add NetCDF C version of mesh conversion tools

conda_package/mpas_tools/logging.py

@@ -3,7 +3,7 @@
 import subprocess
 
 
-def check_call(args, logger, log_command=True, timeout=None, **kwargs):
+def check_call(args, logger=None, log_command=True, timeout=None, **kwargs):
     """
     Call the given subprocess with logging to the given logger.
 
@@ -13,7 +13,7 @@ def check_call(args, logger, log_command=True, timeout=None, **kwargs):
         A list or string of argument to the subprocess.  If ``args`` is a
         string, you must pass ``shell=True`` as one of the ``kwargs``.
 
-    logger : logging.Logger
+    logger : logging.Logger, optional
         The logger to write output to
 
     log_command : bool, optional
@@ -36,25 +36,28 @@ def check_call(args, logger, log_command=True, timeout=None, **kwargs):
         print_args = args
     else:
         print_args = ' '.join(args)
-    if log_command:
-        logger.info(f'Running: {print_args}')
-
-    process = subprocess.Popen(args, stdout=subprocess.PIPE,
-                               stderr=subprocess.PIPE, **kwargs)
-    stdout, stderr = process.communicate(timeout=timeout)
-
-    if stdout:
-        stdout = stdout.decode('utf-8')
-        for line in stdout.split('\n'):
-            logger.info(line)
-    if stderr:
-        stderr = stderr.decode('utf-8')
-        for line in stderr.split('\n'):
-            logger.error(line)
-
-    if process.returncode != 0:
-        raise subprocess.CalledProcessError(process.returncode,
-                                            print_args)
+
+    # make a logger if there isn't already one
+    with LoggingContext(print_args, logger=logger) as logger:
+        if log_command:
+            logger.info(f'Running: {print_args}')
+
+        process = subprocess.Popen(args, stdout=subprocess.PIPE,
+                                   stderr=subprocess.PIPE, **kwargs)
+        stdout, stderr = process.communicate(timeout=timeout)
+
+        if stdout:
+            stdout = stdout.decode('utf-8')
+            for line in stdout.split('\n'):
+                logger.info(line)
+        if stderr:
+            stderr = stderr.decode('utf-8')
+            for line in stderr.split('\n'):
+                logger.error(line)
+
+        if process.returncode != 0:
+            raise subprocess.CalledProcessError(process.returncode,
+                                                print_args)
 
 
 class LoggingContext(object):

conda_package/mpas_tools/mesh/conversion.py

@@ -1,10 +1,11 @@
 import os
 import xarray
-import subprocess
 from tempfile import TemporaryDirectory
 import shutil
 
+import mpas_tools.io
 from mpas_tools.io import write_netcdf
+from mpas_tools.logging import check_call
 
 
 def convert(dsIn, graphInfoFileName=None, logger=None, dir=None):
@@ -48,9 +49,9 @@ def convert(dsIn, graphInfoFileName=None, logger=None, dir=None):
 
         outDir = os.path.dirname(outFileName)
 
-        _call_subprocess(['MpasMeshConverter.x', inFileName, outFileName],
+        check_call(['MpasMeshConverter.x', inFileName, outFileName],
                          logger)
-        
+
         dsOut = xarray.open_dataset(outFileName)
         dsOut.load()
 
@@ -144,8 +145,8 @@ def cull(dsIn, dsMask=None, dsInverse=None, dsPreserve=None,
 
         outDir = os.path.dirname(outFileName)
 
-        _call_subprocess(args, logger)
-        
+        check_call(args=args, logger=logger)
+
         dsOut = xarray.open_dataset(outFileName)
         dsOut.load()
 
@@ -156,10 +157,10 @@ def cull(dsIn, dsMask=None, dsInverse=None, dsPreserve=None,
     return dsOut
 
 
-def mask(dsMesh, fcMask=None, fcSeed=None, logger=None, dir=None):
+def mask(dsMesh, fcMask=None, logger=None, dir=None, cores=1):
     """
-    Use ``MpasMaskCreator.x`` to create a set of region masks either from
-    mask feature collections or from seed points to be used to flood fill
+    Use ``compute_mpas_region_masks`` to create a set of region masks either
+    from mask feature collections
 
     Parameters
     ----------
@@ -169,17 +170,16 @@ def mask(dsMesh, fcMask=None, fcSeed=None, logger=None, dir=None):
     fcMask : geometric_features.FeatureCollection, optional
         A feature collection containing features to use to create the mask
 
-    fcSeed : geometric_features.FeatureCollection, optional
-        A feature collection with points to use a seeds for a flood fill that
-        will create a mask of all cells connected to the seed points
-
     logger : logging.Logger, optional
         A logger for the output if not stdout
 
     dir : str, optional
         A directory in which a temporary directory will be added with files
         produced during mask creation and then deleted upon completion.
 
+    cores : int, optional
+        The number of cores to use for python multiprocessing
+
     Returns
     -------
     dsMask : xarray.Dataset
@@ -189,48 +189,26 @@ def mask(dsMesh, fcMask=None, fcSeed=None, logger=None, dir=None):
     if dir is not None:
         dir = os.path.abspath(dir)
     with TemporaryDirectory(dir=dir) as tempdir:
-        inFileName = '{}/mesh_in.nc'.format(tempdir)
+        inFileName = f'{tempdir}/mesh_in.nc'
         write_netcdf(dsMesh, inFileName)
-        outFileName = '{}/mesh_out.nc'.format(tempdir)
-
-        args = ['MpasMaskCreator.x', inFileName, outFileName]
-
-        if fcMask is not None:
-            fileName = '{}/mask.geojson'.format(tempdir)
-            fcMask.to_geojson(fileName)
-            args.extend(['-f', fileName])
-
-        if fcSeed is not None:
-            fileName = '{}/seed.geojson'.format(tempdir)
-            fcSeed.to_geojson(fileName)
-            args.extend(['-s', fileName])
-
-        _call_subprocess(args, logger)
+        outFileName = f'{tempdir}/mask_out.nc'
+
+        geojsonFileName = f'{tempdir}/mask.geojson'
+        fcMask.to_geojson(geojsonFileName)
+        args = ['compute_mpas_region_masks',
+                '-m', inFileName,
+                '-o', outFileName,
+                '-g', geojsonFileName,
+                '-t', 'cell',
+                '--process_count', f'{cores}',
+                '--format', mpas_tools.io.default_format,
+                ]
+        if mpas_tools.io.default_engine is not None:
+            args.extend(['--engine', mpas_tools.io.default_engine])
+
+        check_call(args=args, logger=logger)
 
         dsOut = xarray.open_dataset(outFileName)
         dsOut.load()
 
     return dsOut
-
-
-def _call_subprocess(args, logger):
-    """Call the given subprocess and send the output to the logger"""
-    if logger is None:
-        subprocess.check_call(args)
-    else:
-        process = subprocess.Popen(args, stdout=subprocess.PIPE,
-                                   stderr=subprocess.PIPE)
-        stdout, stderr = process.communicate()
-
-        if stdout:
-            stdout = stdout.decode('utf-8')
-            for line in stdout.split('\n'):
-                logger.info(line)
-        if stderr:
-            stderr = stderr.decode('utf-8')
-            for line in stderr.split('\n'):
-                logger.error(line)
-
-        if process.returncode != 0:
-            raise subprocess.CalledProcessError(process.returncode,
-                                                ' '.join(args))

conda_package/mpas_tools/mesh/creation/build_mesh.py

@@ -6,8 +6,7 @@
 import cartopy.crs as ccrs
 import cartopy
 
-from mpas_tools.mesh.conversion import convert
-from mpas_tools.io import write_netcdf
+from mpas_tools.logging import check_call
 
 from mpas_tools.mesh.creation.jigsaw_driver import jigsaw_driver
 from mpas_tools.mesh.creation.jigsaw_to_netcdf import jigsaw_to_netcdf
@@ -101,9 +100,10 @@ def build_spherical_mesh(cellWidth, lon, lat, earth_radius,
                          sphere_radius=earth_radius)
 
         logger.info('Step 3. Convert from triangles to MPAS mesh')
-        write_netcdf(convert(xarray.open_dataset('mesh_triangles.nc'), dir=dir,
-                             logger=logger),
-                     out_filename)
+        args = ['MpasMeshConverter.x',
+                'mesh_triangles.nc',
+                out_filename]
+        check_call(args=args, logger=logger)
 
 
 def build_planar_mesh(cellWidth, x, y, geom_points, geom_edges,
@@ -152,7 +152,7 @@ def build_planar_mesh(cellWidth, x, y, geom_points, geom_edges,
                          output_name='mesh_triangles.nc', on_sphere=False)
 
         logger.info('Step 3. Convert from triangles to MPAS mesh')
-        write_netcdf(convert(xarray.open_dataset('mesh_triangles.nc'),
-                             logger=logger),
-                     out_filename)
-
+        args = ['MpasMeshConverter.x',
+                'mesh_triangles.nc',
+                out_filename]
+        check_call(args=args, logger=logger)

conda_package/mpas_tools/mesh/mask.py

@@ -38,7 +38,7 @@ def compute_mpas_region_masks(dsMesh, fcMask, maskTypes=('cell', 'vertex'),
     maskTypes : tuple of {'cell', 'edge', 'vertex'}, optional
         Which type(s) of masks to make.  Masks are created based on whether
         the latitude and longitude associated with each of these locations
-        (e.g. ``dsMesh.latCell`` and ``dsMesh.lonCell`` for ``'cells'``) are
+        (e.g. ``dsMesh.latCell`` and ``dsMesh.lonCell`` for ``'cell'``) are
         inside or outside of the regions in ``fcMask``.
 
     logger : logging.Logger, optional
@@ -213,7 +213,7 @@ def compute_mpas_transect_masks(dsMesh, fcMask, earthRadius,
     maskTypes : tuple of {'cell', 'edge', 'vertex'}, optional
         Which type(s) of masks to make.  Masks are created based on whether
         the latitude and longitude associated with each of these locations
-        (e.g. ``dsMesh.latCell`` and ``dsMesh.lonCell`` for ``'cells'``) are
+        (e.g. ``dsMesh.latCell`` and ``dsMesh.lonCell`` for ``'cell'``) are
         inside or outside of the transects in ``fcMask``.
 
     logger : logging.Logger, optional

conda_package/mpas_tools/viz/paraview_extractor.py

@@ -94,6 +94,7 @@
 from progressbar import ProgressBar, Percentage, Bar, ETA
 from mpas_tools.conversion import mask, cull
 from mpas_tools.io import write_netcdf
+from mpas_tools.logging import check_call
 
 
 def extract_vtk(filename_pattern, variable_list='all', dimension_list=None,
@@ -2146,12 +2147,14 @@ def _cull_files(fc_region_mask, temp_dir, mesh_filename, time_file_names,
         print('Making a region mask file')
         ds_mask = mask(dsMesh=ds_mesh, fcMask=fc_region_mask, logger=logger,
                        dir=temp_dir)
-        write_netcdf(ds_mask, '{}/mask.nc'.format(temp_dir))
+        write_netcdf(ds_mask, f'{temp_dir}/mask.nc')
         print('Cropping mesh to region')
         out_mesh_filename = '{}/mesh.nc'.format(temp_dir)
-        ds_culled = cull(dsIn=ds_mesh, dsInverse=ds_mask, logger=logger,
-                         dir=temp_dir)
-        write_netcdf(ds_culled, out_mesh_filename)
+        args = ['MpasCellCuller.x',
+                mesh_filename,
+                out_mesh_filename,
+                '-i', f'{temp_dir}/mask.nc']
+        check_call(args=args, logger=logger)
 
         region_masks = dict()
         cell_mask = ds_mask.regionCellMasks.sum(dim='nRegions') > 0
@@ -2201,5 +2204,3 @@ def _cull_files(fc_region_mask, temp_dir, mesh_filename, time_file_names,
     handler.close()
 
     return out_mesh_filename, out_time_file_names
-
-# vim: set expandtab:

conda_package/recipe/build.sh

@@ -20,7 +20,7 @@ cmake --build .
 cmake --install .
 
 # build and install mesh conversion tools
-cd ${SRC_DIR}/conda_package/mesh_tools/mesh_conversion_tools
+cd ${SRC_DIR}/conda_package/mesh_tools/mesh_conversion_tools_netcdf_c
 mkdir build
 cd build
 cmake \

conda_package/recipe/meta.yaml

@@ -102,7 +102,7 @@ test:
     - translate_planar_grid -f 'periodic_mesh_10x20_1km.nc' -d 'periodic_mesh_20x40_1km.nc'
     - MpasMeshConverter.x mesh_tools/mesh_conversion_tools/test/mesh.QU.1920km.151026.nc mesh.nc
     - MpasCellCuller.x mesh.nc culled_mesh.nc -m mesh_tools/mesh_conversion_tools/test/land_mask_final.nc
-    - MpasMaskCreator.x mesh.nc arctic_mask.nc -f mesh_tools/mesh_conversion_tools/test/Arctic_Ocean.geojson
+    # - MpasMaskCreator.x mesh.nc arctic_mask.nc -f mesh_tools/mesh_conversion_tools/test/Arctic_Ocean.geojson
     - planar_hex --nx=30 --ny=20 --dc=1000. --npx --npy --outFileName='nonperiodic_mesh_30x20_1km.nc'
     - MpasCellCuller.x nonperiodic_mesh_30x20_1km.nc culled_nonperiodic_mesh_30x20_1km.nc
     - python -m pytest conda_package/mpas_tools/tests

mesh_tools/mesh_conversion_tools_netcdf_c/CMakeLists.txt

@@ -0,0 +1,15 @@
+cmake_minimum_required (VERSION 3.0.2)
+project (mesh_conversion_tools)
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+
+add_executable (MpasMeshConverter.x mpas_mesh_converter.cpp)
+target_link_libraries (MpasMeshConverter.x netcdf)
+
+add_executable (MpasCellCuller.x mpas_cell_culler.cpp)
+target_link_libraries (MpasCellCuller.x netcdf)
+
+#add_executable (MpasMaskCreator.x mpas_mask_creator.cpp jsoncpp.cpp)
+#target_link_libraries (MpasMaskCreator.x netcdf)
+
+install (TARGETS MpasMeshConverter.x MpasCellCuller.x DESTINATION bin)