Update 2 coastal-alteration tools to be functions

The functions are in mpas_tools.ocean.coastline_alteration.
The scripts in ocean/coastline_alteration still work the same as
before by calling the functions.

mpas_tools/ocean/coastline_alteration.py

@@ -0,0 +1,69 @@
+from __future__ import absolute_import, division, print_function, \
+    unicode_literals
+
+import numpy
+
+
+def add_critical_land_blockages(dsMask, dsBlockages):
+    '''
+    Parameters
+    ----------
+    dsMask : `xarray.Dataset`
+        The mask to which critical blockages should be added
+    dsBlockage : `xarray.Dataset`
+        The transect masks defining critical land regions that should block
+        ocean flow (e.g. the Antarctic Peninsula)
+
+    Returns
+    -------
+    dsMask : `xarray.Dataset`
+        The mask with critical blockages included
+    '''
+
+    dsMask = dsMask.copy()
+
+    nTransects = dsBlockages.sizes['nTransects']
+    for transectIndex in range(nTransects):
+        dsMask.regionCellMasks[:, 0] = numpy.maximum(
+            dsBlockages.transectCellMasks[:, transectIndex],
+            dsMask.regionCellMasks[:, 0])
+
+    return dsMask
+
+
+def widen_transect_edge_masks(dsMask, dsMesh, latitude_threshold=43.0):
+    '''
+    Parameters
+    ----------
+    dsMask : `xarray.Dataset`
+        The mask to which critical blockages should be added
+    dsMesh : `xarray.Dataset`
+        The transect masks defining critical land regions that should block
+        ocean flow (e.g. the Antarctic Peninsula)
+    latitude_threshold : float
+        Minimum latitude, degrees, for transect widening
+
+    Returns
+    -------
+    dsMask : `xarray.Dataset`
+        The mask with critical blockages included
+    '''
+    latitude_threshold_radians = numpy.deg2rad(latitude_threshold)
+
+    dsMask = dsMask.copy()
+
+    maxEdges = dsMesh.sizes['maxEdges']
+
+    latMask = numpy.abs(dsMesh.latEdge) > latitude_threshold_radians
+
+    edgeMask = numpy.logical_and(
+        latMask, dsMask.transectEdgeMasks == 1)
+    for iEdge in range(maxEdges):
+        eoc = dsMesh.edgesOnCell[:, iEdge]-1
+        mask = numpy.logical_and(eoc >= 0,
+                                 edgeMask[eoc])
+        # cells with a neighboring transect edge should be masked to 1
+        dsMask['transectCellMasks'] = dsMask.transectCellMasks.where(
+            numpy.logical_not(mask), 1.)
+
+    return dsMask

ocean/coastline_alteration/add_critical_land_blockages_to_mask.py

@@ -11,50 +11,32 @@
 from __future__ import absolute_import, division, print_function, \
     unicode_literals
 
-import os
-import shutil
-from netCDF4 import Dataset
-import numpy as np
+import xarray
 import argparse
 
+from mpas_tools.ocean.coastline_alteration import add_critical_land_blockages
 
-def removeFile(fileName):
-    try:
-        os.remove(fileName)
-    except OSError:
-        pass
-
-
-parser = \
-    argparse.ArgumentParser(description=__doc__,
-                            formatter_class=argparse.RawTextHelpFormatter)
-parser.add_argument("-f", "--input_mask_file", dest="input_mask_filename",
-                    help="Mask file that includes cell and edge masks.",
-                    metavar="INPUTMASKFILE", required=True)
-parser.add_argument("-o", "--output_mask_file", dest="output_mask_filename",
-                    help="Mask file that includes cell and edge masks.",
-                    metavar="OUTPUTMASKFILE", required=True)
-parser.add_argument("-b", "--blockage_file", dest="blockage_file",
-                    help="Masks for each transect identifying critical land"
-                         "blockage.", metavar="BLOCKFILE",
-                    required=True)
-args = parser.parse_args()
-
-removeFile(args.output_mask_filename)
-shutil.copyfile(args.input_mask_filename, args.output_mask_filename)
-
-outMaskFile = Dataset(args.output_mask_filename, "r+")
-nRegions = len(outMaskFile.dimensions["nRegions"])
-regionCellMasks = outMaskFile.variables["regionCellMasks"]
-
-blockageFile = Dataset(args.blockage_file, "r+")
-nTransects = len(blockageFile.dimensions["nTransects"])
-transectCellMasks = blockageFile.variables["transectCellMasks"]
-for transectIndex in range(nTransects):
-    # make sure the regionCellMasks for the first region is 1 anywhere a
-    # transectCellMask is 1
-    regionCellMasks[:, 0] = np.maximum(transectCellMasks[:, transectIndex],
-                                       regionCellMasks[:, 0])
-
-blockageFile.close()
-outMaskFile.close()
+
+if __name__ == '__main__':
+    parser = \
+        argparse.ArgumentParser(description=__doc__,
+                                formatter_class=argparse.RawTextHelpFormatter)
+    parser.add_argument("-f", "--input_mask_file", dest="input_mask_filename",
+                        help="Mask file that includes cell and edge masks.",
+                        metavar="INPUTMASKFILE", required=True)
+    parser.add_argument("-o", "--output_mask_file",
+                        dest="output_mask_filename",
+                        help="Mask file that includes cell and edge masks.",
+                        metavar="OUTPUTMASKFILE", required=True)
+    parser.add_argument("-b", "--blockage_file", dest="blockage_file",
+                        help="Masks for each transect identifying critical "
+                             "land blockage.", metavar="BLOCKFILE",
+                        required=True)
+    args = parser.parse_args()
+
+    dsMask = xarray.open_dataset(args.input_mask_filename)
+
+    dsBlockages = xarray.open_dataset(args.blockage_file)
+
+    dsMask = add_critical_land_blockages(dsMask, dsBlockages)
+    dsMask.to_netcdf(args.output_mask_filename)

ocean/coastline_alteration/mpas_tools

@@ -0,0 +1 @@
+../../mpas_tools/

ocean/coastline_alteration/widen_transect_edge_masks.py

@@ -10,49 +10,38 @@
 from __future__ import absolute_import, division, print_function, \
     unicode_literals
 
-
-import numpy as np
-from netCDF4 import Dataset
 import argparse
-
-parser = \
-    argparse.ArgumentParser(description=__doc__,
-                            formatter_class=argparse.RawTextHelpFormatter)
-parser.add_argument("-f", "--mask_file", dest="mask_filename",
-                    help="Mask file with cell and edge transect masks.",
-                    metavar="MASKFILE",
-                    required=True)
-parser.add_argument("-m", "--mesh_file", dest="mesh_filename",
-                    help="MPAS Mesh filename.", metavar="MESHFILE",
-                    required=True)
-parser.add_argument("-l", "--latitude_threshold", dest="latitude_threshold",
-                    help="Minimum latitude, degrees, for transect widening.",
-                    required=False, type=float, default=43.0)
-args = parser.parse_args()
-
-latitude_threshold_radians = args.latitude_threshold * 3.1415 / 180.
-
-# Obtain mesh variables
-meshFile = Dataset(args.mesh_filename, "r")
-nEdges = len(meshFile.dimensions["nEdges"])
-cellsOnEdge = meshFile.variables["cellsOnEdge"][:, :]
-latEdge = meshFile.variables["latEdge"][:]
-meshFile.close()
-
-# Obtain transect mask variables
-maskFile = Dataset(args.mask_filename, "a")
-nTransects = len(maskFile.dimensions["nTransects"])
-transectCellMasks = maskFile.variables["transectCellMasks"][:, :]
-transectEdgeMasks = maskFile.variables["transectEdgeMasks"][:, :]
-
-print("widen_transect_edge_masks.py: Widening transects to two cells wide")
-for iEdge in range(nEdges):
-    if abs(latEdge[iEdge]) > latitude_threshold_radians:
-        for iTransect in range(nTransects):
-            if transectEdgeMasks[iEdge, iTransect] == 1:
-                maskFile['transectCellMasks'][cellsOnEdge[iEdge, 0] -
-                                              1, iTransect] = 1
-                maskFile['transectCellMasks'][cellsOnEdge[iEdge, 1] -
-                                              1, iTransect] = 1
-
-maskFile.close()
+import xarray
+
+from mpas_tools.ocean.coastline_alteration import widen_transect_edge_masks
+
+
+if __name__ == '__main__':
+
+    parser = \
+        argparse.ArgumentParser(description=__doc__,
+                                formatter_class=argparse.RawTextHelpFormatter)
+    parser.add_argument("-f", "--mask_file", dest="mask_filename",
+                        help="Mask file with cell and edge transect masks.",
+                        metavar="MASKFILE",
+                        required=True)
+    parser.add_argument("-m", "--mesh_file", dest="mesh_filename",
+                        help="MPAS Mesh filename.", metavar="MESHFILE",
+                        required=True)
+    parser.add_argument("-o", "--out_file", dest="out_filename",
+                        help="Output mask file,different from input filename.",
+                        metavar="MASKFILE",
+                        required=True)
+    parser.add_argument("-l", "--latitude_threshold",
+                        dest="latitude_threshold",
+                        help="Minimum latitude, degrees, for transect "
+                             "widening.",
+                        required=False, type=float, default=43.0)
+    args = parser.parse_args()
+
+    dsMask = xarray.open_dataset(args.mask_filename)
+
+    dsMesh = xarray.open_dataset(args.mesh_filename)
+
+    dsMask = widen_transect_edge_masks(dsMask, dsMesh, args.latitude_threshold)
+    dsMask.to_netcdf(args.out_filename)