WIP: Update MOC mask files and task to handle IndoPacific, Pacific and Indian regions

configs/edison/config.20180209.GMPAS-IAF.T62_oRRS30to10v3wLI.cori-knl.afterSalinityFix

@@ -3,7 +3,7 @@
 ## compared against
 
 # mainRunName is a name that identifies the simulation being analyzed.
-mainRunName = 20180129.DECKv1b_piControl.ne30_oEC.edison
+mainRunName = 20180209.GMPAS-IAF.T62_oRRS30to10v3wLI.cori-knl.afterSalinityFix
 
 # preprocessedReferenceRunName is the name of a reference run that has been
 # preprocessed to compare against (or None to turn off comparison).  Reference

mpas_analysis/config.default

@@ -482,6 +482,7 @@ plotVerticalSection = False
 # for MHT vertical section plots
 movingAveragePoints = 1
 
+
 [streamfunctionMOC]
 ## options related to plotting the streamfunction of the meridional overturning
 ## circulation (MOC)
@@ -492,30 +493,8 @@ movingAveragePoints = 1
 includeBolus = True
 
 # Region names for basin MOC calculation.
-# Supported options are Atlantic and IndoPacific
-regionNames = ['Atlantic']
-
-# Size of latitude bins over which MOC streamfunction is integrated
-latBinSizeGlobal = 1.
-latBinSizeAtlantic = 0.5
-latBinSizeIndoPacific = 0.5
-
-# colormap for model results
-colormapNameGlobal = RdYlBu_r
-colormapNameAtlantic = RdYlBu_r
-colormapNameIndoPacific = RdYlBu_r
-# colormap indices for contour color
-colormapIndicesGlobal = [0, 40, 80, 110, 140, 170, 200, 230, 255]
-colormapIndicesAtlantic = [0, 40, 80, 110, 140, 170, 200, 230, 255]
-colormapIndicesIndoPacific = [0, 40, 80, 110, 140, 170, 200, 230, 255]
-# colorbar levels/values for contour boundaries
-colorbarLevelsGlobal = [-20, -10, -5, -2, 2, 5, 10, 20, 30, 40]
-colorbarLevelsAtlantic = [-10, -5, -2, 0, 5, 8, 10, 14, 18, 22]
-colorbarLevelsIndoPacific = [-10, -5, -2, 0, 5, 8, 10, 14, 18, 22]
-# contour line levels
-contourLevelsGlobal = np.arange(-25.1, 35.1, 10)
-contourLevelsAtlantic = np.arange(-8, 20.1, 2)
-contourLevelsIndoPacific = np.arange(-8, 20.1, 2)
+# Supported options are Atlantic, Indian, Pacific and IndoPacific
+regionNames = ['Atlantic', 'IndoPacific']
 
 # Number of points over which to compute moving average for
 # MOC timeseries (e.g., for monthly output, movingAveragePoints=12
@@ -537,6 +516,76 @@ movingAveragePointsClimatological = 1
 # yearStrideXTicks = 1
 
 
+[streamfunctionMOCGlobal]
+# Size of latitude bins over which MOC streamfunction is integrated
+latBinSize = 1.
+
+# colormap for model results
+colormapName = RdYlBu_r
+# colormap indices for contour color
+colormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255]
+# colorbar levels/values for contour boundaries
+colorbarLevels = [-20, -10, -5, -2, 2, 5, 10, 20, 30, 40]
+# contour line levels
+contourLevels = np.arange(-25.1, 35.1, 10)
+
+
+[streamfunctionMOCAtlantic]
+# Size of latitude bins over which MOC streamfunction is integrated
+latBinSize = 0.5
+
+# colormap for model results
+colormapName = RdYlBu_r
+# colormap indices for contour color
+colormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255]
+# colorbar levels/values for contour boundaries
+colorbarLevels = [-10, -5, -2, 0, 5, 8, 10, 14, 18, 22]
+# contour line levels
+contourLevels = np.arange(-8, 20.1, 2)
+
+
+[streamfunctionMOCIndoPacific]
+# Size of latitude bins over which MOC streamfunction is integrated
+latBinSize = 0.5
+
+# colormap for model results
+colormapName = RdYlBu_r
+# colormap indices for contour color
+colormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255]
+# colorbar levels/values for contour boundaries
+colorbarLevels = [-10, -5, -2, 0, 5, 8, 10, 14, 18, 22]
+# contour line levels
+contourLevels = np.arange(-8, 20.1, 2)
+
+
+[streamfunctionMOCIndian]
+# Size of latitude bins over which MOC streamfunction is integrated
+latBinSize = 0.5
+
+# colormap for model results
+colormapName = RdYlBu_r
+# colormap indices for contour color
+colormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255]
+# colorbar levels/values for contour boundaries
+colorbarLevels = [-10, -5, -2, 0, 5, 8, 10, 14, 18, 22]
+# contour line levels
+contourLevels = np.arange(-8, 20.1, 2)
+
+
+[streamfunctionMOCPacific]
+# Size of latitude bins over which MOC streamfunction is integrated
+latBinSize = 0.5
+
+# colormap for model results
+colormapName = RdYlBu_r
+# colormap indices for contour color
+colormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255]
+# colorbar levels/values for contour boundaries
+colorbarLevels = [-10, -5, -2, 0, 5, 8, 10, 14, 18, 22]
+# contour line levels
+contourLevels = np.arange(-8, 20.1, 2)
+
+
 [climatologyMapSST]
 ## options related to plotting horizontally remapped climatologies of
 ## sea surface temperature (SST) against reference model results and

mpas_analysis/ocean/streamfunction_moc.py

@@ -220,8 +220,9 @@ def run_task(self):  # {{{
             x = self.lat[region]
             y = self.depth
             z = self.moc[region]
-            plot_vertical_section(config, x, y, z, self.sectionName,
-                                  suffix=region, colorbarLabel=colorbarLabel,
+            sectionName = '{}{}'.format(self.sectionName, region)
+            plot_vertical_section(config, x, y, z, sectionName,
+                                  colorbarLabel=colorbarLabel,
                                   title=title, xlabel=xLabel, ylabel=yLabel,
                                   fileout=figureName,
                                   N=movingAveragePointsClimatological)
@@ -354,18 +355,34 @@ def _compute_moc_climo_postprocess(self):  # {{{
         mpasMeshName = config.get('input', 'mpasMeshName')
         regionMaskDirectory = config.get('regions', 'regionMaskDirectory')
 
-        regionMaskFile = '{}/{}_SingleRegionAtlanticWTransportTransects_' \
-                         'masks.nc'.format(regionMaskDirectory, mpasMeshName)
+        regionMaskFile = '{}/{}_MOCBasinsAndTransectMasks.nc'.format(
+                regionMaskDirectory, mpasMeshName)
 
         if not os.path.exists(regionMaskFile):
             raise IOError('Regional masking file {} for MOC calculation '
                           'does not exist'.format(regionMaskFile))
-        iRegion = 0
+
+        ncFileRegional = netCDF4.Dataset(regionMaskFile, mode='r')
         self.dictRegion = {}
         for region in self.regionNames:
             self.logger.info('\n  Reading region and transect mask for '
                              '{}...'.format(region))
-            ncFileRegional = netCDF4.Dataset(regionMaskFile, mode='r')
+
+            nRegions = ncFileRegional.dimensions['nRegions'].size
+            found = False
+            for iRegion in range(nRegions):
+                name = ncFileRegional.variables['regionNames'][iRegion]
+                if type(name) == np.ma.MaskedArray:
+                    name = name.filled()
+                name = ''.join([char.decode('UTF-8') for char in name]).strip()
+                if name == '{}_MOC'.format(region):
+                    found = True
+                    break
+
+            if not found:
+                raise ValueError('Region {} not found in MOC region '
+                                 'masks'.format(region))
+
             maxEdgesInTransect = \
                 ncFileRegional.dimensions['maxEdgesInTransect'].size
             transectEdgeMaskSigns = \
@@ -374,7 +391,6 @@ def _compute_moc_climo_postprocess(self):  # {{{
                 ncFileRegional.variables['transectEdgeGlobalIDs'][iRegion, :]
             regionCellMask = \
                 ncFileRegional.variables['regionCellMasks'][:, iRegion]
-            ncFileRegional.close()
             iRegion += 1
 
             indRegion = np.where(regionCellMask == 1)
@@ -384,6 +400,8 @@ def _compute_moc_climo_postprocess(self):  # {{{
                 'maxEdgesInTransect': maxEdgesInTransect,
                 'transectEdgeMaskSigns': transectEdgeMaskSigns,
                 'transectEdgeGlobalIDs': transectEdgeGlobalIDs}
+
+        ncFileRegional.close()
         # Add Global regionCellMask=1 everywhere to make the algorithm
         # for the global moc similar to that of the regional moc
 
@@ -456,9 +474,8 @@ def _compute_moc_climo_postprocess(self):  # {{{
                                                          horizontalVel)
 
                 regionCellMask = self.dictRegion[region]['cellMask']
-                latBinSize = \
-                    config.getExpression(self.sectionName,
-                                         'latBinSize{}'.format(region))
+                sectionName = '{}{}'.format(self.sectionName, region)
+                latBinSize = config.getfloat(sectionName, 'latBinSize')
                 if region == 'Global':
                     latBins = np.arange(-90.0, 90.1, latBinSize)
                 else:

preprocess_masks/make_moc_masks.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+
+'''
+Make a mask file for a given mesh from geometric features describing
+MOC regions.
+
+The -m flag is used to specify the name of the ACME mesh to which the
+masks should be applied.
+
+Requires:
+    * a local link to the MPAS mask creator MpasMaskCreator.x
+    * a local link to the MOC southern boundary extractor tool
+      moc_southern_boundary_extractor.py
+    * a local link to a mesh file named <mesh_name>_mesh.nc describing the
+      desired mesh
+    * the region file MOCBasins.geojson produced by running
+      ./driver_scripts/setup_ocean_region_groups.py in the geometric_features
+      repo
+
+Produces:
+    * <mesh_name>_MOCBasinsAndTransectMasks.nc, the mask file
+
+Author: Xylar Asay-Davis
+'''
+
+import subprocess
+import argparse
+
+
+parser = \
+    argparse.ArgumentParser(description=__doc__,
+                            formatter_class=argparse.RawTextHelpFormatter)
+parser.add_argument('-m', '--mesh_name', dest='mesh_name',
+                    help='The ACME name of the mesh', metavar='MESH_NAME',
+                    required=True)
+args = parser.parse_args()
+
+meshFileName = '{}_mesh.nc'.format(args.mesh_name)
+maskFileName = '{}_MOCBasinsAndTransectMasks.nc'.format(args.mesh_name)
+regionFileName = 'MOCBasins.geojson'
+
+tempRegionMaskFile = 'tempRegionMasks.nc'
+subprocess.check_call(['./MpasMaskCreator.x', meshFileName,
+                       tempRegionMaskFile, '-f', regionFileName])
+
+subprocess.check_call(['./moc_southern_boundary_extractor.py',
+                       '-f', tempRegionMaskFile,
+                       '-m', meshFileName,
+                       '-o', maskFileName])
+
+subprocess.check_call(['rm', tempRegionMaskFile])
+
+