Add a task for plotting the sea surface height vs. obs

mpas_analysis/config.default

@@ -212,6 +212,7 @@ ninoSubdirectory = Nino
 mhtSubdirectory = MHT
 meltSubdirectory = Melt
 soseSubdirectory = SOSE
+sshSubdirectory = SSH
 
 # first and last year of SST observational climatology (preferably one of the
 # two ranges given below)
@@ -582,6 +583,32 @@ seasons =  ['JFM', 'JAS', 'ANN']
 # comparison grid(s) ('latlon', 'antarctic') on which to plot analysis
 comparisonGrids = ['latlon']
 
+[climatologyMapSSH]
+## options related to plotting horizontally remapped climatologies of
+## sea surface height (SSH) against reference model results and observations
+
+# colormap for model/observations
+colormapNameResult = Spectral_r
+# color indices into colormapName for filled contours
+colormapIndicesResult = [0, 20, 40, 80, 100, 120, 160, 180, 200, 240, 255]
+# colormap levels/values for contour boundaries
+colorbarLevelsResult = [-240., -200., -160., -120., -80., -40., 0., 40., 80., 120.]
+
+# colormap for differences
+colormapNameDifference = RdBu_r
+# color indices into colormapName for filled contours
+colormapIndicesDifference = [0, 23, 46, 70, 93, 116, 128, 128, 139, 162, 185, 209, 232, 255]
+# colormap levels/values for contour boundaries
+colorbarLevelsDifference = [-100., -80., -60., -40., -20., -10., 0., 10., 20.,  40.,  60.,  80.,  100.]
+
+# Times for comparison times (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct,
+# Nov, Dec, JFM, AMJ, JAS, OND, ANN)
+seasons =  ['JFM', 'JAS', 'ANN']
+
+# comparison grid(s) ('latlon', 'antarctic') on which to plot analysis
+comparisonGrids = ['latlon']
+
+
 [climatologyMapAntarcticMelt]
 ## options related to plotting horizontally regridded maps of Antarctic
 ## sub-ice-shelf melt rates against reference model results and observations

mpas_analysis/ocean/__init__.py

@@ -1,6 +1,7 @@
 from mpas_analysis.ocean.climatology_map_sst import ClimatologyMapSST
 from mpas_analysis.ocean.climatology_map_mld import ClimatologyMapMLD
 from mpas_analysis.ocean.climatology_map_sss import ClimatologyMapSSS
+from mpas_analysis.ocean.climatology_map_ssh import ClimatologyMapSSH
 from mpas_analysis.ocean.climatology_map_antarctic_melt import \
     ClimatologyMapAntarcticMelt
 from mpas_analysis.ocean.climatology_map_sose import \

mpas_analysis/ocean/climatology_map_ssh.py

@@ -0,0 +1,252 @@
+from __future__ import absolute_import, division, print_function, \
+    unicode_literals
+
+import xarray as xr
+
+from mpas_analysis.shared import AnalysisTask
+
+from mpas_analysis.shared.io.utility import build_config_full_path
+
+from mpas_analysis.shared.climatology import RemapMpasClimatologySubtask, \
+    RemapObservedClimatologySubtask
+
+from mpas_analysis.ocean.plot_climatology_map_subtask import \
+    PlotClimatologyMapSubtask
+
+from mpas_analysis.shared.grid import LatLonGridDescriptor
+
+from mpas_analysis.shared.constants import constants
+
+
+class ClimatologyMapSSH(AnalysisTask):  # {{{
+    """
+    An analysis task for comparison of sea surface height (ssh) against
+    observations
+
+    Authors
+    -------
+    Xylar Asay-Davis
+    """
+    def __init__(self, config, mpasClimatologyTask,
+                 refConfig=None):  # {{{
+        """
+        Construct the analysis task.
+
+        Parameters
+        ----------
+        config :  ``MpasAnalysisConfigParser``
+            Configuration options
+
+        mpasClimatologyTask : ``MpasClimatologyTask``
+            The task that produced the climatology to be remapped and plotted
+
+        refConfig :  ``MpasAnalysisConfigParser``, optional
+            Configuration options for a reference run (if any)
+
+        Authors
+        -------
+        Xylar Asay-Davis
+        """
+        fieldName = 'ssh'
+        # call the constructor from the base class (AnalysisTask)
+        super(ClimatologyMapSSH, self).__init__(
+                config=config, taskName='climatologyMapSSH',
+                componentName='ocean',
+                tags=['climatology', 'horizontalMap', fieldName])
+
+        mpasFieldName = 'timeMonthly_avg_ssh'
+        iselValues = None
+
+        sectionName = self.taskName
+
+        # read in what seasons we want to plot
+        seasons = config.getExpression(sectionName, 'seasons')
+
+        if len(seasons) == 0:
+            raise ValueError('config section {} does not contain valid list '
+                             'of seasons'.format(sectionName))
+
+        comparisonGridNames = config.getExpression(sectionName,
+                                                   'comparisonGrids')
+
+        if len(comparisonGridNames) == 0:
+            raise ValueError('config section {} does not contain valid list '
+                             'of comparison grids'.format(sectionName))
+
+        # the variable mpasFieldName will be added to mpasClimatologyTask
+        # along with the seasons.
+        remapClimatologySubtask = RemapSSHClimatology(
+            mpasClimatologyTask=mpasClimatologyTask,
+            parentTask=self,
+            climatologyName=fieldName,
+            variableList=[mpasFieldName],
+            comparisonGridNames=comparisonGridNames,
+            seasons=seasons,
+            iselValues=iselValues)
+
+        if refConfig is None:
+
+            refTitleLabel = 'Observations (AVISO sea-level anomaly)'
+
+            observationsDirectory = build_config_full_path(
+                config, 'oceanObservations',
+                '{}Subdirectory'.format(fieldName))
+
+            obsFileName = \
+                "{}/zos_AVISO_L4_199210-201012.nc".format(
+                    observationsDirectory)
+            refFieldName = 'zos'
+            outFileLabel = 'sshAVISO'
+            galleryName = 'Observations: AVISO'
+
+            remapObservationsSubtask = RemapObservedSSHClimatology(
+                    parentTask=self, seasons=seasons, fileName=obsFileName,
+                    outFilePrefix=refFieldName,
+                    comparisonGridNames=comparisonGridNames)
+            self.add_subtask(remapObservationsSubtask)
+            diffTitleLabel = 'Model - Observations'
+
+        else:
+            remapObservationsSubtask = None
+            refRunName = refConfig.get('runs', 'mainRunName')
+            galleryName = None
+            refTitleLabel = 'Ref: {}'.format(refRunName)
+
+            refFieldName = mpasFieldName
+            outFileLabel = 'ssh'
+            diffTitleLabel = 'Main - Reference'
+
+        for comparisonGridName in comparisonGridNames:
+            for season in seasons:
+                # make a new subtask for this season and comparison grid
+                subtask = PlotClimatologyMapSubtask(self, season,
+                                                    comparisonGridName,
+                                                    remapClimatologySubtask,
+                                                    remapObservationsSubtask,
+                                                    refConfig)
+
+                subtask.set_plot_info(
+                        outFileLabel=outFileLabel,
+                        fieldNameInTitle='SSH',
+                        mpasFieldName=mpasFieldName,
+                        refFieldName=refFieldName,
+                        refTitleLabel=refTitleLabel,
+                        diffTitleLabel=diffTitleLabel,
+                        unitsLabel=r'cm',
+                        imageCaption='Mean Sea Surface Height',
+                        galleryGroup='Sea Surface Height',
+                        groupSubtitle=None,
+                        groupLink='ssh',
+                        galleryName=galleryName)
+
+                self.add_subtask(subtask)
+        # }}}
+    # }}}
+
+
+class RemapSSHClimatology(RemapMpasClimatologySubtask):  # {{{
+    """
+    Change units from m to cm
+
+    Authors
+    -------
+    Xylar Asay-Davis
+    """
+
+    def customize_masked_climatology(self, climatology):  # {{{
+        """
+        Mask the melt rates using ``landIceMask`` and rescale it to m/yr
+
+        Parameters
+        ----------
+        climatology : ``xarray.Dataset`` object
+            the climatology data set
+
+        Returns
+        -------
+        climatology : ``xarray.Dataset`` object
+            the modified climatology data set
+
+        Authors
+        -------
+        Xylar Asay-Davis
+        """
+
+        fieldName = self.variableList[0]
+
+        # scale the field to cm from m
+        climatology[fieldName] = constants.cm_per_m * climatology[fieldName]
+
+        return climatology  # }}}
+
+    # }}}
+
+
+class RemapObservedSSHClimatology(RemapObservedClimatologySubtask):  # {{{
+    """
+    A subtask for reading and remapping SSH observations
+
+    Authors
+    -------
+    Xylar Asay-Davis
+    """
+
+    def get_observation_descriptor(self, fileName):  # {{{
+        '''
+        get a MeshDescriptor for the observation grid
+
+        Parameters
+        ----------
+        fileName : str
+            observation file name describing the source grid
+
+        Returns
+        -------
+        obsDescriptor : ``MeshDescriptor``
+            The descriptor for the observation grid
+
+        Authors
+        -------
+        Xylar Asay-Davis
+        '''
+
+        # create a descriptor of the observation grid using the lat/lon
+        # coordinates
+        obsDescriptor = LatLonGridDescriptor.read(fileName=fileName,
+                                                  latVarName='lat',
+                                                  lonVarName='lon')
+        return obsDescriptor  # }}}
+
+    def build_observational_dataset(self, fileName):  # {{{
+        '''
+        read in the data sets for observations, and possibly rename some
+        variables and dimensions
+
+        Parameters
+        ----------
+        fileName : str
+            observation file name
+
+        Returns
+        -------
+        dsObs : ``xarray.Dataset``
+            The observational dataset
+
+        Authors
+        -------
+        Xylar Asay-Davis
+        '''
+
+        dsObs = xr.open_dataset(fileName)
+        dsObs.rename({'time': 'Time'}, inplace=True)
+        dsObs.coords['month'] = dsObs['Time.month']
+        dsObs.coords['year'] = dsObs['Time.year']
+
+        # scale the field to cm from m
+        dsObs['zos'] = constants.cm_per_m * dsObs['zos']
+
+        return dsObs  # }}}
+
+    # }}}
+
+# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python

mpas_analysis/shared/constants/constants.py

@@ -64,4 +64,7 @@
 # kilograms per gigatonne
 kg_per_GT = 1e12
 
+# cm per m
+cm_per_m = 100.
+
 # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python

run_mpas_analysis

@@ -85,6 +85,8 @@ def build_analysis_list(config, refConfig):  # {{{
                                             refConfig))
     analyses.append(ocean.ClimatologyMapSSS(config, oceanClimatolgyTask,
                                             refConfig))
+    analyses.append(ocean.ClimatologyMapSSH(config, oceanClimatolgyTask,
+                                            refConfig))
     analyses.append(ocean.ClimatologyMapSoseTemperature(
             config, oceanClimatolgyTask, refConfig))
     analyses.append(ocean.ClimatologyMapSoseSalinity(
@@ -107,6 +109,7 @@ def build_analysis_list(config, refConfig):  # {{{
                                         refConfig))
     analyses.append(ocean.MeridionalHeatTransport(config, oceanClimatolgyTask,
                                                   refConfig))
+
     analyses.append(ocean.StreamfunctionMOC(config, oceanClimatolgyTask))
     analyses.append(ocean.IndexNino34(config, oceanIndexTask, refConfig))