diff --git a/README.md b/README.md index fa3326b2d..973dfa3cd 100644 --- a/README.md +++ b/README.md @@ -27,7 +27,7 @@ available: * scipy * matplotlib * netCDF4 - * xarray >= 0.9.1 + * xarray >= 0.10.0 * dask * bottleneck * basemap diff --git a/configs/anvil/job_script.anvil.bash b/configs/anvil/job_script.anvil.bash index d79ff7cfa..36f3d2e9a 100644 --- a/configs/anvil/job_script.anvil.bash +++ b/configs/anvil/job_script.anvil.bash @@ -13,7 +13,7 @@ cd $PBS_O_WORKDIR # needed to prevent interference with acme-unified unset LD_LIBRARY_PATH -soft add +acme-unified-1.1.1-nox +soft add +e3sm-unified-1.1.2-nox # MPAS/ACME job to be analyzed, including paths to simulation data and # observations. Change this name and path as needed diff --git a/configs/cori/job_script.cori-haswell.bash b/configs/cori/job_script.cori-haswell.bash index f40bd3d59..ec4d6b79c 100644 --- a/configs/cori/job_script.cori-haswell.bash +++ b/configs/cori/job_script.cori-haswell.bash @@ -22,8 +22,7 @@ export OMP_NUM_THREADS=1 module unload python python/base module use /global/project/projectdirs/acme/software/modulefiles/all -module load python/anaconda-2.7-acme -export PATH=/global/homes/z/zender/bin_cori:${PATH} +module load e3sm-unified/1.1.2 # MPAS/ACME job to be analyzed, including paths to simulation data and # observations. Change this name and path as needed diff --git a/configs/cori/job_script.cori-knl.bash b/configs/cori/job_script.cori-knl.bash index c21ef7b61..212af40cc 100644 --- a/configs/cori/job_script.cori-knl.bash +++ b/configs/cori/job_script.cori-knl.bash @@ -22,8 +22,7 @@ export OMP_NUM_THREADS=1 module unload python python/base module use /global/project/projectdirs/acme/software/modulefiles/all -module load python/anaconda-2.7-acme -export PATH=/global/homes/z/zender/bin_cori:${PATH} +module load e3sm-unified/1.1.2 # MPAS/ACME job to be analyzed, including paths to simulation data and # observations. Change this name and path as needed diff --git a/configs/edison/job_script.edison.bash b/configs/edison/job_script.edison.bash index 300731b0b..f1098dc6e 100644 --- a/configs/edison/job_script.edison.bash +++ b/configs/edison/job_script.edison.bash @@ -22,7 +22,7 @@ export OMP_NUM_THREADS=1 module unload python python/base module use /global/project/projectdirs/acme/software/modulefiles/all -module load python/anaconda-2.7-acme +module load e3sm-unified/1.1.2 # MPAS/ACME job to be analyzed, including paths to simulation data and # observations. Change this name and path as needed diff --git a/configs/lanl/job_script.lanl.bash b/configs/lanl/job_script.lanl.bash index 428522253..f42ffb377 100644 --- a/configs/lanl/job_script.lanl.bash +++ b/configs/lanl/job_script.lanl.bash @@ -16,7 +16,7 @@ export OMP_NUM_THREADS=1 module unload python module use /usr/projects/climate/SHARED_CLIMATE/modulefiles/all/ -module load python/anaconda-2.7-climate +module load e3sm-unified/1.1.2 # MPAS/ACME job to be analyzed, including paths to simulation data and # observations. Change this name and path as needed diff --git a/configs/olcf/job_script.olcf.bash b/configs/olcf/job_script.olcf.bash index a3aed1262..42114f0e7 100644 --- a/configs/olcf/job_script.olcf.bash +++ b/configs/olcf/job_script.olcf.bash @@ -16,8 +16,8 @@ cd $PBS_O_WORKDIR module unload python -module use /ccs/proj/cli115/pwolfram/modulefiles/all -module load python/anaconda-2.7-climate +module use /ccs/proj/cli127/software/modulefiles/all +module load e3sm-unified/1.1.2 # MPAS/ACME job to be analyzed, including paths to simulation data and # observations. Change this name and path as needed diff --git a/configs/theta/job_script.theta.bash b/configs/theta/job_script.theta.bash index cdb06d9f4..0b681a7e9 100755 --- a/configs/theta/job_script.theta.bash +++ b/configs/theta/job_script.theta.bash @@ -8,7 +8,7 @@ export OMP_NUM_THREADS=1 module unload python module use /projects/OceanClimate/modulefiles/all -module load python/anaconda-2.7-acme +module load e3sm-unified/1.1.2 # MPAS/ACME job to be analyzed, including paths to simulation data and # observations. Change this name and path as needed diff --git a/docs/api.rst b/docs/api.rst index 3dd80e432..f5e9bbe01 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -55,9 +55,18 @@ Ocean tasks IndexNino34 MeridionalHeatTransport StreamfunctionMOC - TimeSeriesOHC + TimeSeriesOHCAnomaly + TimeSeriesTemperatureAnomaly + TimeSeriesSalinityAnomaly TimeSeriesSST +.. currentmodule:: mpas_analysis.ocean.compute_anomaly_subtask + +.. autosummary:: + :toctree: generated/ + + ComputeAnomalySubtask + .. currentmodule:: mpas_analysis.ocean.plot_climatology_map_subtask .. autosummary:: @@ -66,6 +75,20 @@ Ocean tasks PlotClimatologyMapSubtask PlotClimatologyMapSubtask.set_plot_info +.. currentmodule:: mpas_analysis.ocean.plot_depth_integrated_time_series_subtask + +.. autosummary:: + :toctree: generated/ + + PlotDepthIntegratedTimeSeriesSubtask + +.. currentmodule:: mpas_analysis.ocean.plot_hovmoller_subtask + +.. autosummary:: + :toctree: generated/ + + PlotHovmollerSubtask + Sea ice tasks ------------- @@ -235,7 +258,6 @@ Plotting plotting.plot_1D plotting.plot_vertical_section plotting.setup_colormap - plotting.plot_size_y_axis plotting.plot_xtick_format diff --git a/mpas_analysis/config.default b/mpas_analysis/config.default index 0a39892b4..e751ebc1c 100644 --- a/mpas_analysis/config.default +++ b/mpas_analysis/config.default @@ -295,83 +295,91 @@ baseDirectory = /dir/to/seaice/reference # directory where ocean reference simulation results are stored baseDirectory = /dir/to/seaice/reference -[timeSeriesOHC] +[timeSeriesOHCAnomaly] ## options related to plotting time series of ocean heat content (OHC) +## anomalies from year 1 -# plot S, T, and OHC fields themselves, in addition to their anomalies? -plotOriginalFields = False -# compare to observations? -compareWithObservations = False -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] -# Number of points over which to compute moving average for time series (e.g., -#for monthly output, movingAveragePoints=12 corresponds to a 12-month moving -# average window) -movingAveragePointsTimeSeries = 12 -# Number of points over which to compute moving average for Hovmoller plots -movingAveragePointsHovmoller = 12 - -# colormap for OHC -colormapNameOHC = RdYlBu_r -# colormap indices for contour color -colormapIndicesOHC = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] -# colorbar levels/values for contour boundaries -colorbarLevelsOHC = [0, 5, 10, 15, 20, 25, 30, 35, 40] -# contour line levels -contourLevelsOHC = np.arange(0, 41, 4) +# list of regions to plot from the region list in [regions] below +regions = ['global'] -# colormap for OHC anomaly -colormapNameOHCAnomaly = RdBu_r +# approximate depths (m) separating plots of the upper, middle and lower ocean +depths = [700, 2000] + +# preprocessed file prefix, with format OHC..year*.nc +preprocessedFilePrefix = OHC + +# prefix on preprocessed field name, with format ohc_ for suffixes +# 'tot', '700m', '2000m', 'btm' +preprocessedFieldPrefix = ohc + +# Number of points over which to compute moving average(e.g., for monthly +# output, movingAveragePoints=12 corresponds to a 12-month moving average +# window) +movingAveragePoints = 12 + +[hovmollerOHCAnomaly] +## options related to time vs. depth Hovmoller plots of ocean heat content +## (OHC) anomalies from year 1 + +# Note: regions and moving average points are the same as for the time series +# plot + +# colormap +colormapName = RdBu_r # colormap indices for contour color -colormapIndicesOHCAnomaly = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] +colormapIndices = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] # colorbar levels/values for contour boundaries -colorbarLevelsOHCAnomaly = [-2.4, -0.8, -0.4, -0.2, 0, 0.2, 0.4, 0.8, 2.4] +colorbarLevels = [-2.4, -0.8, -0.4, -0.2, 0, 0.2, 0.4, 0.8, 2.4] # contour line levels -contourLevelsOHCAnomaly = np.arange(-2.5, 2.6, 0.5) +contourLevels = np.arange(-2.5, 2.6, 0.5) -# colormap for temperature -colormapNameTemperature = RdYlBu_r -# color indices into colormapName for filled contours -colormapIndicesTemperature = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] -# colormap levels/values for contour boundaries -colorbarLevelsTemperature = [0, 0.6, 1.2, 2, 4, 6, 8, 10, 20] -# contour line levels -contourLevelsTemperature = np.arange(1, 21, 2) +[hovmollerTemperatureAnomaly] +## options related to plotting time series of temperature vs. depth -# colormap for temperature anomaly -colormapNameTemperatureAnomaly = RdBu_r -# color indices into colormapName for filled contours -colormapIndicesTemperatureAnomaly = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] -# colormap levels/values for contour boundaries -colorbarLevelsTemperatureAnomaly = [-1, -0.5, -0.2, -0.05, 0, 0.05, 0.2, 0.5, 1] -# contour line levels -contourLevelsTemperatureAnomaly = np.arange(-1.0, 1.26, 0.25) +# list of regions to plot from the region list in [regions] below +regions = ['global'] + +# Number of points over which to compute moving average(e.g., for monthly +# output, movingAveragePoints=12 corresponds to a 12-month moving average +# window) +movingAveragePoints = 12 -# colormap for salinity -colormapNameSalinity = RdYlBu_r +# colormap +colormapName = RdBu_r # color indices into colormapName for filled contours -colormapIndicesSalinity = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] +colormapIndices = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] # colormap levels/values for contour boundaries -colorbarLevelsSalinity = [34.5, 34.6, 34.65, 34.7, 34.72, 34.74, 34.76, 34.78, 34.8] +colorbarLevels = [-1, -0.5, -0.2, -0.05, 0, 0.05, 0.2, 0.5, 1] # contour line levels -contourLevelsSalinity = np.arange(34.51, 34.82, 0.02) +contourLevels = np.arange(-1.0, 1.26, 0.25) + +[hovmollerSalinityAnomaly] +## options related to plotting time series of salinity vs. depth + +# list of regions to plot from the region list in [regions] below +regions = ['global'] + +# Number of points over which to compute moving average(e.g., for monthly +# output, movingAveragePoints=12 corresponds to a 12-month moving average +# window) +movingAveragePoints = 12 -# colormap for salinity anomaly -colormapNameSalinityAnomaly = RdBu_r +# colormap +colormapName = RdBu_r # color indices into colormapName for filled contours -colormapIndicesSalinityAnomaly = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] +colormapIndices = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] # colormap levels/values for contour boundaries -colorbarLevelsSalinityAnomaly = [-0.1, -0.02, -0.003, -0.001, 0, 0.001, 0.003, 0.02, 0.1] +colorbarLevels = [-0.1, -0.02, -0.003, -0.001, 0, 0.001, 0.003, 0.02, 0.1] # contour line levels -contourLevelsSalinityAnomaly = np.arange(-0.1, 0.11, 0.02) +contourLevels = np.arange(-0.1, 0.11, 0.02) [timeSeriesSST] ## options related to plotting time series of sea surface temperature (SST) # compare to observations? compareWithObservations = True -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] +# list of regions to plot from the region list in [regions] below +regions = ['global'] # Number of points over which to compute moving average (e.g., for monthly # output, movingAveragePoints=12 corresponds to a 12-month moving average # window) @@ -380,10 +388,9 @@ movingAveragePoints = 12 [indexNino34] ## options related to plotting time series of the El Nino 3.4 index -# Specified region for the Nino Index, 5 = Nino34, 3 = Nino3, 4 = Nino4 -# The indexNino34 routine only accepts one value at a time, -# regionIndicesToPlot should be an integer -regionIndicesToPlot = 5 +# Specified region for the Nino Index,'nino3', 'nino4', or 'nino3.4' +# The indexNino34 routine only accepts one value at a time +region = nino3.4 # Data source to read for comparison. There are three options # 1 - ERS SSTv4 -- Updated version of previous -- 1854 - 2016 @@ -486,8 +493,6 @@ movingAveragePointsClimatological = 1 # compare to observations? compareWithObservations = True -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] # Number of points over which to compute moving average (e.g., for monthly # output, movingAveragePoints=12 corresponds to a 12-month moving average # window) diff --git a/mpas_analysis/ocean/__init__.py b/mpas_analysis/ocean/__init__.py index 06c4b722b..a9e4f9c89 100644 --- a/mpas_analysis/ocean/__init__.py +++ b/mpas_analysis/ocean/__init__.py @@ -5,7 +5,10 @@ from .climatology_map_sose import ClimatologyMapSoseTemperature, \ ClimatologyMapSoseSalinity -from .time_series_ohc import TimeSeriesOHC +from .time_series_temperature_anomaly import TimeSeriesTemperatureAnomaly +from .time_series_salinity_anomaly import TimeSeriesSalinityAnomaly +from .time_series_ohc_anomaly import TimeSeriesOHCAnomaly + from .time_series_sst import TimeSeriesSST from .index_nino34 import IndexNino34 from .streamfunction_moc import StreamfunctionMOC diff --git a/mpas_analysis/ocean/compute_anomaly_subtask.py b/mpas_analysis/ocean/compute_anomaly_subtask.py new file mode 100644 index 000000000..befefbf8c --- /dev/null +++ b/mpas_analysis/ocean/compute_anomaly_subtask.py @@ -0,0 +1,171 @@ +# -*- coding: utf-8 -*- + +from __future__ import absolute_import, division, print_function, \ + unicode_literals + +import os + +from ..shared import AnalysisTask + +from ..shared.io import write_netcdf + +from ..shared.timekeeping.utility import get_simulation_start_time + +from ..shared.io.utility import build_config_full_path + +from ..shared.time_series import compute_moving_avg_anomaly_from_start + + +class ComputeAnomalySubtask(AnalysisTask): + """ + A subtask for computing anomalies of moving averages and writing them out. + + Attributes + ---------- + + mpasTimeSeriesTask : ``MpasTimeSeriesTask`` + The task that extracts the time series from MPAS monthly output + + outFileName : str + The file name (usually without full path) where the resulting + data set should be written + + variableList : list of str + Variables to be included in the data set + + movingAveragePoints : int + The number of points (months) used in the moving average used to + smooth the data set + + alter_dataset : function + A function that takes an ``xarray.Dataset`` and returns an + ``xarray.Dataset`` for manipulating the data set (e.g. adding a new + variable computed from others). This operation is performed before + computing moving averages and anomalies, so that these operations are + also performed on any new variables added to the data set. + + Authors + ------- + Xylar Asay-Davis + """ + + def __init__(self, parentTask, mpasTimeSeriesTask, outFileName, + variableList, movingAveragePoints, + subtaskName='computeAnomaly', alter_dataset=None): # {{{ + """ + Construct the analysis task. + + Parameters + ---------- + parentTask : ``AnalysisTask`` + The parent task of which this is a subtask + + mpasTimeSeriesTask : ``MpasTimeSeriesTask`` + The task that extracts the time series from MPAS monthly output + + outFileName : str + The file name (usually without full path) where the resulting + data set should be written + + variableList : list of str + Variables to be included in the data set + + movingAveragePoints : int + The number of points (months) used in the moving average used to + smooth the data set + + subtaskName : str, optional + The name of the subtask + + alter_dataset : function + A function that takes an ``xarray.Dataset`` and returns an + ``xarray.Dataset`` for manipulating the data set (e.g. adding a new + variable computed from others). This operation is performed before + computing moving averages and anomalies, so that these operations + are also performed on any new variables added to the data set. + + Authors + ------- + Xylar Asay-Davis + """ + # first, call the constructor from the base class (AnalysisTask) + super(ComputeAnomalySubtask, self).__init__( + config=parentTask.config, + taskName=parentTask.taskName, + componentName='ocean', + tags=parentTask.tags, + subtaskName=subtaskName) + + self.mpasTimeSeriesTask = mpasTimeSeriesTask + + self.run_after(mpasTimeSeriesTask) + + self.outFileName = outFileName + self.variableList = variableList + self.movingAveragePoints = movingAveragePoints + + self.alter_dataset = alter_dataset + + # }}} + + def setup_and_check(self): # {{{ + """ + Perform steps to set up the analysis and check for errors in the setup. + + Authors + ------- + Xylar Asay-Davis + """ + # first, call setup_and_check from the base class (AnalysisTask), + # which will perform some common setup, including storing: + # self.runDirectory , self.historyDirectory, self.plotsDirectory, + # self.namelist, self.runStreams, self.historyStreams, + # self.calendar + super(ComputeAnomalySubtask, self).setup_and_check() + + self.mpasTimeSeriesTask.add_variables(variableList=self.variableList) + + self.inputFile = self.mpasTimeSeriesTask.outputFile + + # }}} + + def run_task(self): # {{{ + """ + Performs analysis of ocean heat content (OHC) from time-series output. + + Authors + ------- + Xylar Asay-Davis, Milena Veneziani, Greg Streletz + """ + + self.logger.info("\nComputing anomalies...") + + config = self.config + startDate = config.get('timeSeries', 'startDate') + endDate = config.get('timeSeries', 'endDate') + + simulationStartTime = get_simulation_start_time(self.runStreams) + + ds = compute_moving_avg_anomaly_from_start( + timeSeriesFileName=self.inputFile, + variableList=self.variableList, + simulationStartTime=simulationStartTime, + startDate=startDate, + endDate=endDate, + calendar=self.calendar, + movingAveragePoints=self.movingAveragePoints, + alter_dataset=self.alter_dataset) + + outFileName = self.outFileName + if not os.path.isabs(outFileName): + baseDirectory = build_config_full_path( + config, 'output', 'timeSeriesSubdirectory') + + outFileName = '{}/{}'.format(baseDirectory, + outFileName) + + write_netcdf(ds, outFileName) # }}} + + # }}} + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/index_nino34.py b/mpas_analysis/ocean/index_nino34.py index 366c2bfee..5b8fec209 100644 --- a/mpas_analysis/ocean/index_nino34.py +++ b/mpas_analysis/ocean/index_nino34.py @@ -145,7 +145,9 @@ def run_task(self): # {{{ # regionIndex should correspond to NINO34 in surface weighted Average # AM - regionIndex = config.getint('indexNino34', 'regionIndicesToPlot') + regions = config.getExpression('regions', 'regions') + regionToPlot = config.get('indexNino34', 'region') + regionIndex = regions.index(regionToPlot) # Load data: ds = open_mpas_dataset(fileName=self.inputFile, diff --git a/mpas_analysis/ocean/plot_depth_integrated_time_series_subtask.py b/mpas_analysis/ocean/plot_depth_integrated_time_series_subtask.py new file mode 100644 index 000000000..3ff1dc013 --- /dev/null +++ b/mpas_analysis/ocean/plot_depth_integrated_time_series_subtask.py @@ -0,0 +1,400 @@ +# -*- coding: utf-8 -*- + +from __future__ import absolute_import, division, print_function, \ + unicode_literals + +import xarray as xr +import os + +from ..shared import AnalysisTask + +from ..shared.plot.plotting import timeseries_analysis_plot + +from ..shared.generalized_reader import open_multifile_dataset +from ..shared.io import open_mpas_dataset + +from ..shared.timekeeping.utility import date_to_days, days_to_datetime + +from ..shared.io.utility import build_config_full_path + +from ..shared.html import write_image_xml + +from ..shared.time_series import compute_moving_avg + + +class PlotDepthIntegratedTimeSeriesSubtask(AnalysisTask): + """ + Plots a time series, summed or averaged over various depth ranges + + Attributes + ---------- + + regionName : str + The name of the region to plot + + inFileName : str + The file containing the time-depth data set to plot + + outFileLabel : str + The prefix on each plot and associated XML file + + fieldNameInTitle : str + The name of the field being plotted, as used in the plot title + + mpasFieldName : str + The name of the variable in the MPAS timeSeriesStatsMonthly output + + yAxisLabel : str + the y-axis label of the plotted field (including units) + + sectionName : str + A section in the config file where the colormap and contour values + are defined + + thumbnailSuffix : str + The text to be displayed under the thumbnail image, to which the + region name will be prepended + + imageCaption : str + The caption when mousing over the plot or displaying it full + screen + + galleryGroup : str + The name of the group of galleries in which this plot belongs + + groupSubtitle : str + The subtitle of the group in which this plot belongs (or blank + if none) + + groupLink : str + A short name (with no spaces) for the link to the gallery group + + galleryName : str + The name of the gallery in which this plot belongs + + + Authors + ------- + Xylar Asay-Davis, Milena Veneziani, Greg Streletz + """ + + def __init__(self, parentTask, regionName, inFileName, outFileLabel, + fieldNameInTitle, mpasFieldName, yAxisLabel, sectionName, + thumbnailSuffix, imageCaption, galleryGroup, groupSubtitle, + groupLink, galleryName, subtaskName=None): # {{{ + """ + Construct the analysis task. + + Parameters + ---------- + parentTask : ``AnalysisTask`` + The parent task of which this is a subtask + + regionName : str + The name of the region to plot + + inFileName : str + The file containing the time-depth data set to plot + + outFileLabel : str + The prefix on each plot and associated XML file + + fieldNameInTitle : str + The name of the field being plotted, as used in the plot title + + mpasFieldName : str + The name of the variable in the MPAS timeSeriesStatsMonthly output + + yAxisLabel : str + the y-axis label of the plotted field + + sectionName : str + a section in the config file where the colormap and contour values + are defined + + thumbnailSuffix : str + The text to be displayed under the thumbnail image, to which the + region name will be prepended + + imageCaption : str + the caption when mousing over the plot or displaying it full + screen + + galleryGroup : str + the name of the group of galleries in which this plot belongs + + groupSubtitle : str + the subtitle of the group in which this plot belongs (or blank + if none) + + groupLink : str + a short name (with no spaces) for the link to the gallery group + + galleryName : str + the name of the gallery in which this plot belongs + + subtaskName : str + The name of the subtask (``plotTimeSeries`` by default) + + Authors + ------- + Xylar Asay-Davis + """ + + if subtaskName is None: + suffix = regionName[0].upper() + regionName[1:] + subtaskName = 'plotDepthIntegratedTimeSeries{}'.format(suffix) + + # first, call the constructor from the base class (AnalysisTask) + super(PlotDepthIntegratedTimeSeriesSubtask, self).__init__( + config=parentTask.config, + taskName=parentTask.taskName, + componentName='ocean', + tags=parentTask.tags, + subtaskName=subtaskName) + + self.regionName = regionName + self.inFileName = inFileName + self.outFileLabel = outFileLabel + self.fieldNameInTitle = fieldNameInTitle + self.mpasFieldName = mpasFieldName + self.yAxisLabel = yAxisLabel + self.sectionName = sectionName + + # xml/html related variables + self.thumbnailSuffix = thumbnailSuffix + self.imageCaption = imageCaption + self.galleryGroup = galleryGroup + self.groupSubtitle = groupSubtitle + self.groupLink = groupLink + self.galleryName = galleryName + + # }}} + + def setup_and_check(self): # {{{ + """ + Perform steps to set up the analysis and check for errors in the setup. + + Authors + ------- + Xylar Asay-Davis, Greg Streletz + """ + # first, call setup_and_check from the base class (AnalysisTask), + # which will perform some common setup, including storing: + # self.runDirectory , self.historyDirectory, self.plotsDirectory, + # self.namelist, self.runStreams, self.historyStreams, + # self.calendar + super(PlotDepthIntegratedTimeSeriesSubtask, self).setup_and_check() + + config = self.config + + if not os.path.isabs(self.inFileName): + baseDirectory = build_config_full_path( + config, 'output', 'timeSeriesSubdirectory') + + self.inFileName = '{}/{}'.format(baseDirectory, + self.inFileName) + + mainRunName = self.config.get('runs', 'mainRunName') + + self.filePrefix = '{}_{}_{}'.format(self.outFileLabel, + self.regionName, + mainRunName) + self.xmlFileNames = ['{}/{}.xml'.format( + self.plotsDirectory, self.filePrefix)] + + return # }}} + + def run_task(self): # {{{ + """ + Compute vertical agregates of the data and plot the time series + + Authors + ------- + Xylar Asay-Davis, Milena Veneziani, Greg Streletz + """ + + self.logger.info("\nPlotting depth-integrated time series of " + "{}...".format(self.fieldNameInTitle)) + + config = self.config + calendar = self.calendar + + mainRunName = config.get('runs', 'mainRunName') + + plotTitles = config.getExpression('regions', 'plotTitles') + allRegionNames = config.getExpression('regions', 'regions') + regionIndex = allRegionNames.index(self.regionName) + regionNameInTitle = plotTitles[regionIndex] + + startDate = config.get('timeSeries', 'startDate') + endDate = config.get('timeSeries', 'endDate') + + # Load data + self.logger.info(' Load ocean data...') + ds = open_mpas_dataset(fileName=self.inFileName, + calendar=calendar, + variableList=[self.mpasFieldName], + timeVariableNames=None, + startDate=startDate, + endDate=endDate) + ds = ds.isel(nOceanRegionsTmp=regionIndex) + + # Note: restart file, not a mesh file because we need refBottomDepth, + # not in a mesh file + try: + restartFile = self.runStreams.readpath('restart')[0] + except ValueError: + raise IOError('No MPAS-O restart file found: need at least one ' + 'restart file for OHC calculation') + + # Define/read in general variables + self.logger.info(' Read in depth...') + with xr.open_dataset(restartFile) as dsRestart: + # reference depth [m] + depths = dsRestart.refBottomDepth.values + + # add depths as a coordinate to the data set + ds.coords['depth'] = (('nVertLevels',), depths) + + divisionDepths = config.getExpression(self.sectionName, 'depths') + + # for each depth interval to plot, determine the top and bottom depth + topDepths = [0, 0] + divisionDepths + bottomDepths = [depths[-1]] + divisionDepths + [depths[-1]] + + legends = [] + for top, bottom in zip(topDepths, bottomDepths): + if bottom == depths[-1]: + legends.append('{}m-bottom'.format(top)) + else: + legends.append('{}m-{}m'.format(top, bottom)) + + # more possible symbols than we typically use + lines = ['-', '-', '--', '+', 'o', '^', 'v'] + widths = [5, 3, 3, 3, 3, 3, 3] + points = [None, None, None, 300, 300, 300, 300] + + color = 'k' + + xLabel = 'Time [years]' + yLabel = self.yAxisLabel + + title = '{}, {} \n {} (black)'.format(self.fieldNameInTitle, + regionNameInTitle, + mainRunName) + + figureName = '{}/{}.png'.format(self.plotsDirectory, self.filePrefix) + + timeSeries = [] + lineStyles = [] + lineWidths = [] + maxPoints = [] + legendText = [] + + for rangeIndex in range(len(topDepths)): + top = topDepths[rangeIndex] + bottom = bottomDepths[rangeIndex] + field = ds[self.mpasFieldName].where(ds.depth > top) + field = field.where(ds.depth <= bottom) + timeSeries.append(field.sum('nVertLevels')) + + lineStyles.append('{}{}'.format(color, lines[rangeIndex])) + lineWidths.append(widths[rangeIndex]) + maxPoints.append(points[rangeIndex]) + legendText.append(legends[rangeIndex]) + + preprocessedReferenceRunName = config.get( + 'runs', 'preprocessedReferenceRunName') + if preprocessedReferenceRunName != 'None': + preprocessedInputDirectory = config.get( + 'oceanPreprocessedReference', 'baseDirectory') + + self.logger.info(' Load in OHC from preprocessed reference ' + 'run...') + preprocessedFilePrefix = config.get(self.sectionName, + 'preprocessedFilePrefix') + inFilesPreprocessed = '{}/{}.{}.year*.nc'.format( + preprocessedInputDirectory, preprocessedFilePrefix, + preprocessedReferenceRunName) + dsPreprocessed = open_multifile_dataset( + fileNames=inFilesPreprocessed, + calendar=calendar, + config=config, + timeVariableName='xtime') + + yearStart = days_to_datetime(ds.Time.min(), calendar=calendar).year + yearEnd = days_to_datetime(ds.Time.max(), calendar=calendar).year + timeStart = date_to_days(year=yearStart, month=1, day=1, + calendar=calendar) + timeEnd = date_to_days(year=yearEnd, month=12, day=31, + calendar=calendar) + + yearEndPreprocessed = days_to_datetime(dsPreprocessed.Time.max(), + calendar=calendar).year + if yearStart <= yearEndPreprocessed: + dsPreprocessed = dsPreprocessed.sel(Time=slice(timeStart, + timeEnd)) + else: + self.logger.warning('Warning: Preprocessed time series ends ' + 'before the timeSeries startYear and will ' + 'not be plotted.') + preprocessedReferenceRunName = 'None' + + if preprocessedReferenceRunName != 'None': + color = 'r' + title = '{} \n {} (red)'.format(title, + preprocessedReferenceRunName) + + preprocessedFieldPrefix = config.get(self.sectionName, + 'preprocessedFieldPrefix') + + movingAveragePoints = config.getint(self.sectionName, + 'movingAveragePoints') + + suffixes = ['tot'] + ['{}m'.format(depth) for depth in + divisionDepths] + ['btm'] + + # these preprocessed data are OHC *anomalies* + for rangeIndex in range(len(suffixes)): + variableName = '{}_{}'.format(preprocessedFieldPrefix, + suffixes[rangeIndex]) + if variableName in list(dsPreprocessed.data_vars.keys()): + field = dsPreprocessed[variableName] + field = compute_moving_avg(field, movingAveragePoints) + timeSeries.append(field) + else: + self.logger.warning('Warning: Preprocessed variable {} ' + 'not found. Skipping.'.format( + variableName)) + timeSeries.extend(None) + + lineStyles.append('{}{}'.format(color, lines[rangeIndex])) + lineWidths.append(widths[rangeIndex]) + maxPoints.append(points[rangeIndex]) + legendText.append(None) + + timeseries_analysis_plot(config=config, dsvalues=timeSeries, N=None, + title=title, xlabel=xLabel, ylabel=yLabel, + fileout=figureName, lineStyles=lineStyles, + lineWidths=lineWidths, maxPoints=maxPoints, + legendText=legendText, calendar=calendar) + + write_image_xml( + config=config, + filePrefix=self.filePrefix, + componentName='Ocean', + componentSubdirectory='ocean', + galleryGroup=self.galleryGroup, + groupLink=self.groupLink, + galleryName=self.galleryName, + thumbnailDescription='{} {}'.format(self.regionName, + self.thumbnailSuffix), + imageDescription=self.imageCaption, + imageCaption=self.imageCaption) + + # }}} + + # }}} + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/plot_hovmoller_subtask.py b/mpas_analysis/ocean/plot_hovmoller_subtask.py new file mode 100644 index 000000000..2c802c8d7 --- /dev/null +++ b/mpas_analysis/ocean/plot_hovmoller_subtask.py @@ -0,0 +1,292 @@ +# -*- coding: utf-8 -*- + +from __future__ import absolute_import, division, print_function, \ + unicode_literals + +import xarray as xr +import os + +from ..shared import AnalysisTask + +from ..shared.plot.plotting import plot_vertical_section, setup_colormap + +from ..shared.io import open_mpas_dataset + +from ..shared.io.utility import build_config_full_path + +from ..shared.html import write_image_xml + + +class PlotHovmollerSubtask(AnalysisTask): + """ + Plots a time series vs. depth + + Attributes + ---------- + + regionName : str + The name of the region to plot + + inFileName : str + The file containing the time-depth data set to plot + + outFileLabel : str + The prefix on each plot and associated XML file + + fieldNameInTitle : str + The name of the field being plotted, as used in the plot title + + mpasFieldName : str + The name of the variable in the MPAS timeSeriesStatsMonthly output + + unitsLabel : str + The units of the plotted field, to be displayed on color bars + + sectionName : str + A section in the config file where the colormap and contour values + are defined + + thumbnailSuffix : str + The text to be displayed under the thumbnail image, to which the + region name will be prepended + + imageCaption : str + The caption when mousing over the plot or displaying it full + screen + + galleryGroup : str + The name of the group of galleries in which this plot belongs + + groupSubtitle : str + The subtitle of the group in which this plot belongs (or blank + if none) + + groupLink : str + A short name (with no spaces) for the link to the gallery group + + galleryName : str + The name of the gallery in which this plot belongs + + + Authors + ------- + Xylar Asay-Davis, Milena Veneziani, Greg Streletz + """ + + def __init__(self, parentTask, regionName, inFileName, outFileLabel, + fieldNameInTitle, mpasFieldName, unitsLabel, sectionName, + thumbnailSuffix, imageCaption, galleryGroup, groupSubtitle, + groupLink, galleryName, subtaskName=None): # {{{ + """ + Construct the analysis task. + + Parameters + ---------- + parentTask : ``AnalysisTask`` + The parent task of which this is a subtask + + regionName : str + The name of the region to plot + + inFileName : str + The file containing the time-depth data set to plot + + outFileLabel : str + The prefix on each plot and associated XML file + + fieldNameInTitle : str + The name of the field being plotted, as used in the plot title + + mpasFieldName : str + The name of the variable in the MPAS timeSeriesStatsMonthly output + + unitsLabel : str + the units of the plotted field, to be displayed on color bars + + sectionName : str + a section in the config file where the colormap and contour values + are defined + + thumbnailSuffix : str + The text to be displayed under the thumbnail image, to which the + region name will be prepended + + imageCaption : str + the caption when mousing over the plot or displaying it full + screen + + galleryGroup : str + the name of the group of galleries in which this plot belongs + + groupSubtitle : str + the subtitle of the group in which this plot belongs (or blank + if none) + + groupLink : str + a short name (with no spaces) for the link to the gallery group + + galleryName : str + the name of the gallery in which this plot belongs + + subtaskName : str + The name of the subtask (``plotHovmoller`` by default) + + Authors + ------- + Xylar Asay-Davis + """ + + if subtaskName is None: + suffix = regionName[0].upper() + regionName[1:] + subtaskName = 'plotHovmoller{}'.format(suffix) + + # first, call the constructor from the base class (AnalysisTask) + super(PlotHovmollerSubtask, self).__init__( + config=parentTask.config, + taskName=parentTask.taskName, + componentName='ocean', + tags=parentTask.tags, + subtaskName=subtaskName) + + self.regionName = regionName + self.inFileName = inFileName + self.outFileLabel = outFileLabel + self.fieldNameInTitle = fieldNameInTitle + self.mpasFieldName = mpasFieldName + self.unitsLabel = unitsLabel + self.sectionName = sectionName + + # xml/html related variables + self.thumbnailSuffix = thumbnailSuffix + self.imageCaption = imageCaption + self.galleryGroup = galleryGroup + self.groupSubtitle = groupSubtitle + self.groupLink = groupLink + self.galleryName = galleryName + + # }}} + + def setup_and_check(self): # {{{ + """ + Perform steps to set up the analysis and check for errors in the setup. + + Authors + ------- + Xylar Asay-Davis, Greg Streletz + """ + # first, call setup_and_check from the base class (AnalysisTask), + # which will perform some common setup, including storing: + # self.runDirectory , self.historyDirectory, self.plotsDirectory, + # self.namelist, self.runStreams, self.historyStreams, + # self.calendar + super(PlotHovmollerSubtask, self).setup_and_check() + + config = self.config + + if not os.path.isabs(self.inFileName): + baseDirectory = build_config_full_path( + config, 'output', 'timeSeriesSubdirectory') + + self.inFileName = '{}/{}'.format(baseDirectory, + self.inFileName) + + mainRunName = self.config.get('runs', 'mainRunName') + + self.filePrefix = '{}_{}_{}'.format(self.outFileLabel, + self.regionName, + mainRunName) + self.xmlFileNames = ['{}/{}.xml'.format( + self.plotsDirectory, self.filePrefix)] + + return # }}} + + def run_task(self): # {{{ + """ + Make the Hovmoller plot from the time series. + + Authors + ------- + Xylar Asay-Davis, Milena Veneziani, Greg Streletz + """ + + self.logger.info("\nPlotting {} trends vs. depth...".format( + self.fieldNameInTitle)) + + config = self.config + + mainRunName = config.get('runs', 'mainRunName') + + plotTitles = config.getExpression('regions', 'plotTitles') + allRegionNames = config.getExpression('regions', 'regions') + regionIndex = allRegionNames.index(self.regionName) + regionNameInTitle = plotTitles[regionIndex] + + startDate = self.config.get('timeSeries', 'startDate') + endDate = self.config.get('timeSeries', 'endDate') + + # Load data + self.logger.info(' Load ocean data...') + ds = open_mpas_dataset(fileName=self.inFileName, + calendar=self.calendar, + variableList=[self.mpasFieldName], + timeVariableNames=None, + startDate=startDate, + endDate=endDate) + ds = ds.isel(nOceanRegionsTmp=regionIndex) + + # Note: restart file, not a mesh file because we need refBottomDepth, + # not in a mesh file + try: + restartFile = self.runStreams.readpath('restart')[0] + except ValueError: + raise IOError('No MPAS-O restart file found: need at least one ' + 'restart file for OHC calculation') + + # Define/read in general variables + self.logger.info(' Read in depth...') + with xr.open_dataset(restartFile) as dsRestart: + # reference depth [m] + depth = dsRestart.refBottomDepth.values + + Time = ds.Time.values + field = ds[self.mpasFieldName].values.transpose() + + xLabel = 'Time [years]' + yLabel = 'Depth [m]' + + title = '{}, {} \n {}'.format(self.fieldNameInTitle, regionNameInTitle, + mainRunName) + + figureName = '{}/{}.png'.format(self.plotsDirectory, self.filePrefix) + + (colormapName, colorbarLevels) = setup_colormap(config, + self.sectionName) + + contourLevels = config.getExpression(self.sectionName, + 'contourLevels', + usenumpyfunc=True) + + plot_vertical_section(config, Time, depth, field, + colormapName, colorbarLevels, contourLevels, + self.unitsLabel, title, xLabel, yLabel, + figureName, linewidths=1, xArrayIsTime=True, + calendar=self.calendar) + + write_image_xml( + config=config, + filePrefix=self.filePrefix, + componentName='Ocean', + componentSubdirectory='ocean', + galleryGroup=self.galleryGroup, + groupLink=self.groupLink, + galleryName=self.galleryName, + thumbnailDescription='{} {}'.format(self.regionName, + self.thumbnailSuffix), + imageDescription=self.imageCaption, + imageCaption=self.imageCaption) + + # }}} + + # }}} + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/time_series_ohc.py b/mpas_analysis/ocean/time_series_ohc.py deleted file mode 100644 index b701d6e48..000000000 --- a/mpas_analysis/ocean/time_series_ohc.py +++ /dev/null @@ -1,690 +0,0 @@ -# -*- coding: utf-8 -*- - -from __future__ import absolute_import, division, print_function, \ - unicode_literals - -import numpy as np -import netCDF4 - -from ..shared import AnalysisTask - -from ..shared.plot.plotting import timeseries_analysis_plot, \ - plot_vertical_section, setup_colormap - -from ..shared.generalized_reader import open_multifile_dataset -from ..shared.io import open_mpas_dataset - -from ..shared.timekeeping.utility import get_simulation_start_time, \ - date_to_days, days_to_datetime, string_to_datetime - -from ..shared.io.utility import build_config_full_path, make_directories, \ - check_path_exists -from ..shared.html import write_image_xml - - -class TimeSeriesOHC(AnalysisTask): - """ - Performs analysis of ocean heat content (OHC) from time-series output. - - Attributes - ---------- - - mpasTimeSeriesTask : ``MpasTimeSeriesTask`` - The task that extracts the time series from MPAS monthly output - - Authors - ------- - Xylar Asay-Davis, Milena Veneziani, Greg Streletz - """ - - def __init__(self, config, mpasTimeSeriesTask): # {{{ - """ - Construct the analysis task. - - Parameters - ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options - - mpasTimeSeriesTask : ``MpasTimeSeriesTask`` - The task that extracts the time series from MPAS monthly output - - Authors - ------- - Xylar Asay-Davis - """ - # first, call the constructor from the base class (AnalysisTask) - super(TimeSeriesOHC, self).__init__( - config=config, - taskName='timeSeriesOHC', - componentName='ocean', - tags=['timeSeries', 'ohc']) - - self.mpasTimeSeriesTask = mpasTimeSeriesTask - - self.run_after(mpasTimeSeriesTask) - - # }}} - - def setup_and_check(self): # {{{ - """ - Perform steps to set up the analysis and check for errors in the setup. - - Raises - ------ - OSError - If files are not present - - Authors - ------- - Xylar Asay-Davis, Greg Streletz - """ - # first, call setup_and_check from the base class (AnalysisTask), - # which will perform some common setup, including storing: - # self.runDirectory , self.historyDirectory, self.plotsDirectory, - # self.namelist, self.runStreams, self.historyStreams, - # self.calendar - super(TimeSeriesOHC, self).setup_and_check() - - config = self.config - - self.startDate = self.config.get('timeSeries', 'startDate') - self.endDate = self.config.get('timeSeries', 'endDate') - - self.variables = {} - for suffix in ['avgLayerTemperature', 'avgLayerSalinity', - 'sumLayerMaskValue', 'avgLayerArea', - 'avgLayerThickness']: - self.variables[suffix] = \ - 'timeMonthly_avg_avgValueWithinOceanLayerRegion_' + suffix - - self.mpasTimeSeriesTask.add_variables( - variableList=self.variables.values()) - - self.inputFile = self.mpasTimeSeriesTask.outputFile - - if config.get('runs', 'preprocessedReferenceRunName') != 'None': - check_path_exists(config.get('oceanPreprocessedReference', - 'baseDirectory')) - - self.startDate = self.config.get('timeSeries', 'startDate') - self.endDate = self.config.get('timeSeries', 'endDate') - - mainRunName = config.get('runs', 'mainRunName') - regions = config.getExpression('regions', 'regions') - regionIndicesToPlot = config.getExpression('timeSeriesOHC', - 'regionIndicesToPlot') - - self.xmlFileNames = [] - self.filePrefixes = {} - - regions = [regions[index] for index in regionIndicesToPlot] - - configSectionName = 'timeSeriesOHC' - plotOriginalFields = config.getboolean(configSectionName, - 'plotOriginalFields') - - if plotOriginalFields: - plotTypes = ['TZ', 'TAnomalyZ', 'SZ', 'SAnomalyZ', 'OHCZ', - 'OHCAnomalyZ', 'OHC', 'OHCAnomaly'] - else: - plotTypes = ['TAnomalyZ', 'SAnomalyZ', 'OHCAnomalyZ', 'OHCAnomaly'] - - for region in regions: - for plotType in plotTypes: - filePrefix = '{}_{}_{}'.format(plotType, region, mainRunName) - self.xmlFileNames.append('{}/{}.xml'.format( - self.plotsDirectory, filePrefix)) - self.filePrefixes[plotType, region] = filePrefix - - return # }}} - - def run_task(self): # {{{ - """ - Performs analysis of ocean heat content (OHC) from time-series output. - - Authors - ------- - Xylar Asay-Davis, Milena Veneziani, Greg Streletz - """ - - self.logger.info("\nPlotting OHC time series and T, S, and OHC " - "vertical trends...") - - simulationStartTime = get_simulation_start_time(self.runStreams) - config = self.config - calendar = self.calendar - variables = self.variables - - # read parameters from config file - mainRunName = config.get('runs', 'mainRunName') - preprocessedReferenceRunName = \ - config.get('runs', 'preprocessedReferenceRunName') - preprocessedInputDirectory = config.get('oceanPreprocessedReference', - 'baseDirectory') - - configSectionName = 'timeSeriesOHC' - - plotOriginalFields = config.getboolean(configSectionName, - 'plotOriginalFields') - - compareWithObservations = config.getboolean(configSectionName, - 'compareWithObservations') - - movingAveragePointsTimeSeries = config.getint( - configSectionName, 'movingAveragePointsTimeSeries') - - movingAveragePointsHovmoller = config.getint( - configSectionName, 'movingAveragePointsHovmoller') - - regions = config.getExpression('regions', 'regions') - plotTitles = config.getExpression('regions', 'plotTitles') - regionIndicesToPlot = config.getExpression(configSectionName, - 'regionIndicesToPlot') - - outputDirectory = build_config_full_path(config, 'output', - 'timeseriesSubdirectory') - - make_directories(outputDirectory) - - regionNames = config.getExpression('regions', 'regions') - - # Note: input file, not a mesh file because we need dycore specific - # fields such as refBottomDepth and namelist fields such as - # config_density0, as well as simulationStartTime, that are not - # guaranteed to be in the mesh file. - try: - restartFile = self.runStreams.readpath('restart')[0] - except ValueError: - raise IOError('No MPAS-O restart file found: need at least one ' - 'restart file for OHC calculation') - - # Define/read in general variables - self.logger.info(' Read in depth and compute specific depth ' - 'indexes...') - ncFile = netCDF4.Dataset(restartFile, mode='r') - # reference depth [m] - depth = ncFile.variables['refBottomDepth'][:] - ncFile.close() - - k700m = np.where(depth > 700.)[0][0] - 1 - k2000m = np.where(depth > 2000.)[0][0] - 1 - - kbtm = len(depth)-1 - - # Load data - self.logger.info(' Load ocean data...') - ds = open_mpas_dataset(fileName=self.inputFile, - calendar=calendar, - variableList=self.variables.values(), - startDate=self.startDate, - endDate=self.endDate) - # rename the variables to shorter names for convenience - renameDict = dict((v, k) for k, v in variables.items()) - ds.rename(renameDict, inplace=True) - - timeStart = string_to_datetime(self.startDate) - timeEnd = string_to_datetime(self.endDate) - - # Select year-1 data and average it (for later computing anomalies) - timeStartFirstYear = string_to_datetime(simulationStartTime) - if timeStartFirstYear < timeStart: - startDateFirstYear = simulationStartTime - firstYear = int(startDateFirstYear[0:4]) - endDateFirstYear = '{:04d}-12-31_23:59:59'.format(firstYear) - dsFirstYear = open_mpas_dataset( - fileName=self.inputFile, - calendar=calendar, - variableList=self.variables.values(), - startDate=startDateFirstYear, - endDate=endDateFirstYear) - - dsFirstYear.rename(renameDict, inplace=True) - - firstYearAvgLayerTemperature = dsFirstYear['avgLayerTemperature'] - firstYearavgLayerSalinity = dsFirstYear['avgLayerSalinity'] - else: - firstYearAvgLayerTemperature = ds['avgLayerTemperature'] - firstYearavgLayerSalinity = ds['avgLayerSalinity'] - firstYear = timeStart.year - - timeStartFirstYear = date_to_days(year=firstYear, month=1, day=1, - calendar=calendar) - timeEndFirstYear = date_to_days(year=firstYear, month=12, day=31, - hour=23, minute=59, second=59, - calendar=calendar) - - firstYearAvgLayerTemperature = firstYearAvgLayerTemperature.sel( - Time=slice(timeStartFirstYear, timeEndFirstYear)) - firstYearAvgLayerTemperature = \ - firstYearAvgLayerTemperature.mean('Time') - - firstYearavgLayerSalinity = firstYearavgLayerSalinity.sel( - Time=slice(timeStartFirstYear, timeEndFirstYear)) - firstYearavgLayerSalinity = \ - firstYearavgLayerSalinity.mean('Time') - - self.logger.info(' Compute temperature and salinity anomalies...') - - ds['avgLayerTemperatureAnomaly'] = \ - (ds['avgLayerTemperature'] - firstYearAvgLayerTemperature) - - ds['avgLayerSalinityAnomaly'] = \ - (ds['avgLayerSalinity'] - firstYearavgLayerSalinity) - - yearStart = days_to_datetime(ds.Time.min(), calendar=calendar).year - yearEnd = days_to_datetime(ds.Time.max(), calendar=calendar).year - timeStart = date_to_days(year=yearStart, month=1, day=1, - calendar=calendar) - timeEnd = date_to_days(year=yearEnd, month=12, day=31, - calendar=calendar) - - if preprocessedReferenceRunName != 'None': - self.logger.info(' Load in OHC from preprocessed reference ' - 'run...') - inFilesPreprocessed = '{}/OHC.{}.year*.nc'.format( - preprocessedInputDirectory, preprocessedReferenceRunName) - dsPreprocessed = open_multifile_dataset( - fileNames=inFilesPreprocessed, - calendar=calendar, - config=config, - timeVariableName='xtime') - yearEndPreprocessed = days_to_datetime(dsPreprocessed.Time.max(), - calendar=calendar).year - if yearStart <= yearEndPreprocessed: - dsPreprocessedTimeSlice = \ - dsPreprocessed.sel(Time=slice(timeStart, timeEnd)) - else: - self.logger.warning('Preprocessed time series ends before the ' - 'timeSeries startYear and will not be ' - 'plotted.') - preprocessedReferenceRunName = 'None' - - # Add the OHC to the data set - ds = self._compute_ohc(ds, regionNames) - - unitsScalefactor = 1e-22 - - self.logger.info(' Compute OHC and make plots...') - for regionIndex in regionIndicesToPlot: - region = regions[regionIndex] - - # Plot temperature, salinity and OHC anomalies (and full fields) - # trends with depth - - # First T vs depth/time - x = ds.Time.values - y = depth - z = ds.avgLayerTemperatureAnomaly.isel( - nOceanRegionsTmp=regionIndex) - z = z.transpose() - - colorbarLabel = '[$^\circ$C]' - xLabel = 'Time [years]' - yLabel = 'Depth [m]' - - title = 'Temperature Anomaly, {} \n {}'.format(plotTitles[regionIndex], mainRunName) - - figureName = '{}/{}.png'.format(self.plotsDirectory, self.filePrefixes['TAnomalyZ', region]) - - (colormapName, colorbarLevels) = setup_colormap(config, - configSectionName, - suffix='TemperatureAnomaly') - - contourLevels = \ - config.getExpression(configSectionName, - 'contourLevels{}'.format('TemperatureAnomaly'), - usenumpyfunc=True) - - plot_vertical_section(config, x, y, z, - colormapName, colorbarLevels, contourLevels, - colorbarLabel, title, xLabel, yLabel, - figureName, linewidths=1, xArrayIsTime=True, - N=movingAveragePointsHovmoller, calendar=calendar) - - caption = 'Trend of {} Temperature Anomaly vs depth from Year ' \ - '{}'.format(region, simulationStartTime[0:4]) - write_image_xml( - config=config, - filePrefix=self.filePrefixes['TAnomalyZ', region], - componentName='Ocean', - componentSubdirectory='ocean', - galleryGroup='Trends vs Depth', - groupLink='trendsvsdepth', - thumbnailDescription=u'{} ΔT'.format(region), - imageDescription=caption, - imageCaption=caption) - - if plotOriginalFields: - z = ds['avgLayerTemperature'].isel(nOceanRegionsTmp=regionIndex) - z = z.transpose() - - title = 'Temperature, {} \n {}'.format(plotTitles[regionIndex], mainRunName) - - figureName = '{}/{}.png'.format(self.plotsDirectory, self.filePrefixes['TZ', region]) - - (colormap, colorbarLevels) = setup_colormap(config, - configSectionName, - suffix='Temperature') - - contourLevels = config.getExpression(configSectionName, - 'contourLevels{}'.format('Temperature'), - usenumpyfunc=True) - - plot_vertical_section(config, x, y, z, - colormap, colorbarLevels, contourLevels, - colorbarLabel, title, xLabel, yLabel, - figureName, linewidths=1, xArrayIsTime=True, - N=movingAveragePointsHovmoller, calendar=calendar) - - caption = 'Trend of {} Temperature vs depth'.format(region) - write_image_xml( - config=config, - filePrefix=self.filePrefixes['TZ', region], - componentName='Ocean', - componentSubdirectory='ocean', - galleryGroup='Trends vs Depth', - groupLink='trendsvsdepth', - thumbnailDescription=u'{} T'.format(region), - imageDescription=caption, - imageCaption=caption) - - - # Second S vs depth/time - z = ds.avgLayerSalinityAnomaly.isel(nOceanRegionsTmp=regionIndex) - z = z.transpose() - - title = 'Salinity Anomaly, {} \n {}'.format(plotTitles[regionIndex], mainRunName) - - colorbarLabel = '[PSU]' - - figureName = '{}/{}.png'.format(self.plotsDirectory, self.filePrefixes['SAnomalyZ', region]) - - (colormapName, colorbarLevels) = setup_colormap(config, - configSectionName, - suffix='SalinityAnomaly') - - contourLevels = \ - config.getExpression(configSectionName, - 'contourLevels{}'.format('SalinityAnomaly'), - usenumpyfunc=True) - - plot_vertical_section(config, x, y, z, - colormapName, colorbarLevels, contourLevels, - colorbarLabel, title, xLabel, yLabel, - figureName, linewidths=1, xArrayIsTime=True, - N=movingAveragePointsHovmoller, calendar=calendar) - - caption = 'Trend of {} Salinity Anomaly vs depth from Year ' \ - '{}'.format(region, simulationStartTime[0:4]) - write_image_xml( - config=config, - filePrefix=self.filePrefixes['SAnomalyZ', region], - componentName='Ocean', - componentSubdirectory='ocean', - galleryGroup='Trends vs Depth', - groupLink='trendsvsdepth', - thumbnailDescription=u'{} ΔS'.format(region), - imageDescription=caption, - imageCaption=caption) - - if plotOriginalFields: - z = ds['avgLayerSalinity'].isel(nOceanRegionsTmp=regionIndex) - z = z.transpose() - - title = 'Salinity, {} \n {}'.format(plotTitles[regionIndex], mainRunName) - - figureName = '{}/{}.png'.format(self.plotsDirectory, self.filePrefixes['SZ', region]) - - (colormapName, colorbarLevels) = setup_colormap(config, - configSectionName, - suffix='Salinity') - - contourLevels = config.getExpression(configSectionName, - 'contourLevels{}'.format('Salinity'), - usenumpyfunc=True) - - plot_vertical_section(config, x, y, z, - colormapName, colorbarLevels, contourLevels, - colorbarLabel, title, xLabel, yLabel, - figureName, linewidths=1, xArrayIsTime=True, - N=movingAveragePointsHovmoller, calendar=calendar) - - caption = 'Trend of {} Salinity vs depth'.format(region) - write_image_xml( - config=config, - filePrefix=self.filePrefixes['SZ', region], - componentName='Ocean', - componentSubdirectory='ocean', - galleryGroup='Trends vs Depth', - groupLink='trendsvsdepth', - thumbnailDescription=u'{} S'.format(region), - imageDescription=caption, - imageCaption=caption) - - - # Third OHC vs depth/time - ohcAnomaly = ds.ohcAnomaly.isel(nOceanRegionsTmp=regionIndex) - ohcAnomalyScaled = unitsScalefactor*ohcAnomaly - z = ohcAnomalyScaled.transpose() - - title = 'OHC Anomaly, {} \n {}'.format(plotTitles[regionIndex], mainRunName) - - colorbarLabel = '[x$10^{22}$ J]' - - figureName = '{}/{}.png'.format(self.plotsDirectory, self.filePrefixes['OHCAnomalyZ', region]) - - (colormap, colorbarLevels) = setup_colormap(config, - configSectionName, - suffix='OHCAnomaly') - - contourLevels = \ - config.getExpression(configSectionName, - 'contourLevels{}'.format('OHCAnomaly'), - usenumpyfunc=True) - - plot_vertical_section(config, x, y, z, - colormap, colorbarLevels, contourLevels, - colorbarLabel, title, xLabel, yLabel, - figureName, linewidths=1, xArrayIsTime=True, - N=movingAveragePointsHovmoller, calendar=calendar) - - caption = 'Trend of {} OHC Anomaly vs depth from Year ' \ - '{}'.format(region, simulationStartTime[0:4]) - write_image_xml( - config=config, - filePrefix=self.filePrefixes['OHCAnomalyZ', region], - componentName='Ocean', - componentSubdirectory='ocean', - galleryGroup='Trends vs Depth', - groupLink='trendsvsdepth', - thumbnailDescription=u'{} ΔOHC'.format(region), - imageDescription=caption, - imageCaption=caption) - - if plotOriginalFields: - ohc = ds.ohc.isel(nOceanRegionsTmp=regionIndex) - ohcScaled = unitsScalefactor*ohc - z = ohcScaled.transpose() - - title = 'OHC, {} \n {}'.format(plotTitles[regionIndex], mainRunName) - - figureName = '{}/{}.png'.format(self.plotsDirectory, self.filePrefixes['OHCZ', region]) - - (colormap, colorbarLevels) = setup_colormap(config, - configSectionName, - suffix='OHC') - - contourLevels = config.getExpression(configSectionName, - 'contourLevels{}'.format('OHC'), - usenumpyfunc=True) - - plot_vertical_section(config, x, y, z, - colormap, colorbarLevels, contourLevels, - colorbarLabel, title, xLabel, yLabel, - figureName, linewidths=1, xArrayIsTime=True, - N=movingAveragePointsHovmoller, calendar=calendar) - - caption = 'Trend of {} OHC vs depth'.format(region) - write_image_xml( - config=config, - filePrefix=self.filePrefixes['OHCZ', region], - componentName='Ocean', - componentSubdirectory='ocean', - galleryGroup='Trends vs Depth', - groupLink='trendsvsdepth', - thumbnailDescription=u'{} OHC'.format(region), - imageDescription=caption, - imageCaption=caption) - - - # Now plot OHC timeseries - - # OHC over 0-bottom depth range: - ohcAnomalyTotal = unitsScalefactor*ohcAnomaly.sum('nVertLevels') - - # OHC over 0-700m depth range: - ohcAnomaly700m = unitsScalefactor*ohcAnomaly[:, 0:k700m].sum('nVertLevels') - - # OHC over 700m-2000m depth range: - ohcAnomaly2000m = \ - unitsScalefactor*ohcAnomaly[:, k700m+1:k2000m].sum('nVertLevels') - - # OHC over 2000m-bottom depth range: - ohcAnomalyBottom = unitsScalefactor*ohcAnomaly[:, k2000m+1:kbtm].sum('nVertLevels') - - xLabel = 'Time [years]' - yLabel = '$\Delta$OHC [x$10^{22}$ J]' - - title = 'OHC anomaly, {} \n {}'.format(plotTitles[regionIndex], - mainRunName) - - figureName = '{}/{}.png'.format(self.plotsDirectory, - self.filePrefixes['OHCAnomaly', region]) - - if preprocessedReferenceRunName != 'None': - # these preprocessed data are OHC *anomalies* - ohcPreprocessedTotal = dsPreprocessedTimeSlice.ohc_tot - ohcPreprocessed700m = dsPreprocessedTimeSlice.ohc_700m - ohcPreprocessed2000m = dsPreprocessedTimeSlice.ohc_2000m - ohcPreprocessedBottom = dsPreprocessedTimeSlice.ohc_btm - title = '{} (black lines) \n {} (red lines)'.format( - title, preprocessedReferenceRunName) - timeseries_analysis_plot(config, [ohcAnomalyTotal, - ohcAnomaly700m, - ohcAnomaly2000m, - ohcAnomalyBottom, - ohcPreprocessedTotal, - ohcPreprocessed700m, - ohcPreprocessed2000m, - ohcPreprocessedBottom], - movingAveragePointsTimeSeries, title, - xLabel, yLabel, figureName, - lineStyles=['k-', 'k-', 'k--', 'k+', - 'r-', 'r-', 'r--', 'r+'], - lineWidths=[5, 3, 3, 3, - 5, 3, 3, 3], - maxPoints=[None, None, None, 300, - None, None, None, 300], - legendText=['0-bottom', '0-700m', - '700-2000m', '2000m-bottom', - None, None, None, None], - calendar=calendar) - - if (not compareWithObservations and - preprocessedReferenceRunName == 'None'): - timeseries_analysis_plot(config, [ohcAnomalyTotal, - ohcAnomaly700m, - ohcAnomaly2000m, - ohcAnomalyBottom], - movingAveragePointsTimeSeries, title, - xLabel, yLabel, figureName, - lineStyles=['k-', 'k-', 'k--', 'k+'], - lineWidths=[5, 3, 3, 3], - maxPoints=[None, None, None, 300], - legendText=['0-bottom', '0-700m', - '700-2000m', '2000m-bottom'], - calendar=calendar) - - caption = 'Running Mean of the Anomaly in {} Ocean Heat Content ' \ - 'from Year {}'.format(region, simulationStartTime[0:4]) - write_image_xml( - config=config, - filePrefix=self.filePrefixes['OHCAnomaly', region], - componentName='Ocean', - componentSubdirectory='ocean', - galleryGroup='Time Series', - groupLink='timeseries', - thumbnailDescription=u'{} ΔOHC'.format(region), - imageDescription=caption, - imageCaption=caption) - - if plotOriginalFields: - ohcTotal = unitsScalefactor*ohc.sum('nVertLevels') - ohc700m = unitsScalefactor*ohc[:, 0:k700m].sum('nVertLevels') - ohc2000m = \ - unitsScalefactor*ohc[:, k700m+1:k2000m].sum('nVertLevels') - ohcBottom = unitsScalefactor*ohc[:, k2000m+1:kbtm].sum('nVertLevels') - - title = 'OHC, {} \n {}'.format(plotTitles[regionIndex], - mainRunName) - - figureName = '{}/{}.png'.format(self.plotsDirectory, - self.filePrefixes['OHC', region]) - - timeseries_analysis_plot(config, [ohcTotal, - ohc700m, - ohc2000m, - ohcBottom], - movingAveragePointsTimeSeries, title, - xLabel, yLabel, figureName, - lineStyles=['k-', 'k-', 'k--', 'k+'], - lineWidths=[5, 3, 3, 3], - legendText=['0-bottom', '0-700m', - '700-2000m', '2000m-bottom'], - calendar=calendar) - - caption = 'Running Mean of {} Ocean Heat Content'.format(region) - write_image_xml( - config=config, - filePrefix=self.filePrefixes['OHC', region], - componentName='Ocean', - componentSubdirectory='ocean', - galleryGroup='Time Series', - groupLink='timeseries', - thumbnailDescription=u'{} OHC'.format(region), - imageDescription=caption, - imageCaption=caption) - - # }}} - - def _compute_ohc(self, ds, regionNames): # {{{ - ''' - Compute the OHC time series. - ''' - - # specific heat [J/(kg*degC)] - cp = self.namelist.getfloat('config_specific_heat_sea_water') - # [kg/m3] - rho = self.namelist.getfloat('config_density0') - - ds['ohc'] = rho*cp*ds['sumLayerMaskValue'] * \ - ds['avgLayerArea'] * ds['avgLayerThickness'] * \ - ds['avgLayerTemperature'] - ds.ohc.attrs['units'] = 'J' - ds.ohc.attrs['description'] = 'Ocean heat content in each region' - - ds['ohcAnomaly'] = rho*cp*ds['sumLayerMaskValue'] * \ - ds['avgLayerArea'] * ds['avgLayerThickness'] * \ - ds.avgLayerTemperatureAnomaly - ds.ohcAnomaly.attrs['units'] = 'J' - ds.ohcAnomaly.attrs['description'] = \ - 'Ocean heat content anomaly in each region' - - ds['regionNames'] = ('nOceanRegionsTmp', regionNames) - - return ds # }}} - - # }}} - -# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/time_series_ohc_anomaly.py b/mpas_analysis/ocean/time_series_ohc_anomaly.py new file mode 100644 index 000000000..cab550a5e --- /dev/null +++ b/mpas_analysis/ocean/time_series_ohc_anomaly.py @@ -0,0 +1,142 @@ +# -*- coding: utf-8 -*- + +from __future__ import absolute_import, division, print_function, \ + unicode_literals + +from ..shared import AnalysisTask + +from .compute_anomaly_subtask import ComputeAnomalySubtask +from .plot_hovmoller_subtask import PlotHovmollerSubtask +from .plot_depth_integrated_time_series_subtask import \ + PlotDepthIntegratedTimeSeriesSubtask + + +class TimeSeriesOHCAnomaly(AnalysisTask): + """ + Performs analysis of ocean heat content (OHC) from time-series output. + + Authors + ------- + Xylar Asay-Davis, Milena Veneziani, Greg Streletz + """ + + def __init__(self, config, mpasTimeSeriesTask): # {{{ + """ + Construct the analysis task. + + Parameters + ---------- + config : instance of MpasAnalysisConfigParser + Contains configuration options + + mpasTimeSeriesTask : ``MpasTimeSeriesTask`` + The task that extracts the time series from MPAS monthly output + + Authors + ------- + Xylar Asay-Davis + """ + # first, call the constructor from the base class (AnalysisTask) + super(TimeSeriesOHCAnomaly, self).__init__( + config=config, + taskName='timeSeriesOHCAnomaly', + componentName='ocean', + tags=['timeSeries', 'ohc']) + + sectionName = 'timeSeriesOHCAnomaly' + regionNames = config.getExpression(sectionName, 'regions') + movingAveragePoints = config.getint(sectionName, 'movingAveragePoints') + + self.variableDict = {} + for suffix in ['avgLayerTemperature', 'sumLayerMaskValue', + 'avgLayerArea', 'avgLayerThickness']: + key = 'timeMonthly_avg_avgValueWithinOceanLayerRegion_' + suffix + self.variableDict[key] = suffix + + mpasFieldName = 'ohc' + + timeSeriesFileName = 'regionAveragedOHCAnomaly.nc' + + anomalyTask = ComputeAnomalySubtask( + parentTask=self, + mpasTimeSeriesTask=mpasTimeSeriesTask, + outFileName=timeSeriesFileName, + variableList=list(self.variableDict.keys()), + movingAveragePoints=movingAveragePoints, + alter_dataset=self._compute_ohc) + self.add_subtask(anomalyTask) + + for regionName in regionNames: + caption = 'Trend of {} OHC Anomaly vs depth'.format( + regionName) + plotTask = PlotHovmollerSubtask( + parentTask=self, + regionName=regionName, + inFileName=timeSeriesFileName, + outFileLabel='ohcAnomalyZ', + fieldNameInTitle='OHC Anomaly', + mpasFieldName=mpasFieldName, + unitsLabel=r'[$\times 10^{22}$ J]', + sectionName='hovmollerOHCAnomaly', + thumbnailSuffix=u'ΔOHC', + imageCaption=caption, + galleryGroup='Trends vs Depth', + groupSubtitle=None, + groupLink='trendsvsdepth', + galleryName=None) + + plotTask.run_after(anomalyTask) + self.add_subtask(plotTask) + + caption = 'Running Mean of the Anomaly in {} Ocean Heat ' \ + 'Content'.format(regionName) + plotTask = PlotDepthIntegratedTimeSeriesSubtask( + parentTask=self, + regionName=regionName, + inFileName=timeSeriesFileName, + outFileLabel='ohcAnomaly', + fieldNameInTitle='OHC Anomaly', + mpasFieldName=mpasFieldName, + yAxisLabel=r'$\Delta$OHC [$\times 10^{22}$ J]', + sectionName='timeSeriesOHCAnomaly', + thumbnailSuffix=u'ΔOHC', + imageCaption=caption, + galleryGroup='Time Series', + groupSubtitle=None, + groupLink='timeseries', + galleryName=None) + + plotTask.run_after(anomalyTask) + self.add_subtask(plotTask) + + # }}} + + def _compute_ohc(self, ds): # {{{ + ''' + Compute the OHC time series. + ''' + + # regionNames = self.config.getExpression('regions', 'regions') + # ds['regionNames'] = ('nOceanRegionsTmp', regionNames) + + # for convenience, rename the variables to simpler, shorter names + ds = ds.rename(self.variableDict) + + # specific heat [J/(kg*degC)] + cp = self.namelist.getfloat('config_specific_heat_sea_water') + # [kg/m3] + rho = self.namelist.getfloat('config_density0') + + unitsScalefactor = 1e-22 + + ds['ohc'] = unitsScalefactor*rho*cp*ds['sumLayerMaskValue'] * \ + ds['avgLayerArea'] * ds['avgLayerThickness'] * \ + ds['avgLayerTemperature'] + ds.ohc.attrs['units'] = '$10^{22}$ J' + ds.ohc.attrs['description'] = 'Ocean heat content in each region' + + return ds # }}} + + # }}} + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/time_series_salinity_anomaly.py b/mpas_analysis/ocean/time_series_salinity_anomaly.py new file mode 100644 index 000000000..dba021329 --- /dev/null +++ b/mpas_analysis/ocean/time_series_salinity_anomaly.py @@ -0,0 +1,88 @@ +# -*- coding: utf-8 -*- + +from __future__ import absolute_import, division, print_function, \ + unicode_literals + +from ..shared import AnalysisTask + +from .compute_anomaly_subtask import ComputeAnomalySubtask +from .plot_hovmoller_subtask import PlotHovmollerSubtask + + +class TimeSeriesSalinityAnomaly(AnalysisTask): + """ + Performs analysis of time series of salinity anomalies from the first + simulation year as a function of depth. + + Authors + ------- + Xylar Asay-Davis + """ + + def __init__(self, config, mpasTimeSeriesTask): # {{{ + """ + Construct the analysis task. + + Parameters + ---------- + config : instance of MpasAnalysisConfigParser + Contains configuration options + + mpasTimeSeriesTask : ``MpasTimeSeriesTask`` + The task that extracts the time series from MPAS monthly output + + Authors + ------- + Xylar Asay-Davis + """ + # first, call the constructor from the base class (AnalysisTask) + super(TimeSeriesSalinityAnomaly, self).__init__( + config=config, + taskName='timeSeriesSalinityAnomaly', + componentName='ocean', + tags=['timeSeries', 'salinity']) + + sectionName = 'hovmollerSalinityAnomaly' + regionNames = config.getExpression(sectionName, 'regions') + movingAveragePoints = config.getint(sectionName, 'movingAveragePoints') + + mpasFieldName = 'timeMonthly_avg_avgValueWithinOceanLayerRegion_' \ + 'avgLayerSalinity' + + timeSeriesFileName = 'regionAveragedSalinityAnomaly.nc' + + anomalyTask = ComputeAnomalySubtask( + parentTask=self, + mpasTimeSeriesTask=mpasTimeSeriesTask, + outFileName=timeSeriesFileName, + variableList=[mpasFieldName], + movingAveragePoints=movingAveragePoints) + self.add_subtask(anomalyTask) + + for regionName in regionNames: + caption = 'Trend of {} Salinity Anomaly vs depth'.format( + regionName) + plotTask = PlotHovmollerSubtask( + parentTask=self, + regionName=regionName, + inFileName=timeSeriesFileName, + outFileLabel='SAnomalyZ', + fieldNameInTitle='Salinity Anomaly', + mpasFieldName=mpasFieldName, + unitsLabel='[PSU]', + sectionName=sectionName, + thumbnailSuffix=u'ΔS', + imageCaption=caption, + galleryGroup='Trends vs Depth', + groupSubtitle=None, + groupLink='trendsvsdepth', + galleryName=None) + + plotTask.run_after(anomalyTask) + self.add_subtask(plotTask) + + # }}} + + # }}} + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/time_series_sst.py b/mpas_analysis/ocean/time_series_sst.py index 1d671140f..bcdeda694 100644 --- a/mpas_analysis/ocean/time_series_sst.py +++ b/mpas_analysis/ocean/time_series_sst.py @@ -96,15 +96,11 @@ def setup_and_check(self): # {{{ self.inputFile = self.mpasTimeSeriesTask.outputFile mainRunName = config.get('runs', 'mainRunName') - regions = config.getExpression('regions', 'regions') - regionIndicesToPlot = config.getExpression('timeSeriesSST', - 'regionIndicesToPlot') + regions = config.getExpression('timeSeriesSST', 'regions') self.xmlFileNames = [] self.filePrefixes = {} - regions = [regions[index] for index in regionIndicesToPlot] - for region in regions: filePrefix = 'sst_{}_{}'.format(region, mainRunName) self.xmlFileNames.append('{}/{}.xml'.format(self.plotsDirectory, @@ -141,17 +137,16 @@ def run_task(self): # {{{ regions = config.getExpression('regions', 'regions') plotTitles = config.getExpression('regions', 'plotTitles') - regionIndicesToPlot = config.getExpression('timeSeriesSST', - 'regionIndicesToPlot') + regionsToPlot = config.getExpression('timeSeriesSST', 'regions') + + regionIndicesToPlot = [regions.index(region) for region in + regionsToPlot] outputDirectory = build_config_full_path(config, 'output', 'timeseriesSubdirectory') make_directories(outputDirectory) - regionNames = config.getExpression('regions', 'regions') - regionNames = [regionNames[index] for index in regionIndicesToPlot] - dsSST = open_mpas_dataset(fileName=self.inputFile, calendar=calendar, variableList=self.variableList, diff --git a/mpas_analysis/ocean/time_series_temperature_anomaly.py b/mpas_analysis/ocean/time_series_temperature_anomaly.py new file mode 100644 index 000000000..4897795f3 --- /dev/null +++ b/mpas_analysis/ocean/time_series_temperature_anomaly.py @@ -0,0 +1,88 @@ +# -*- coding: utf-8 -*- + +from __future__ import absolute_import, division, print_function, \ + unicode_literals + +from ..shared import AnalysisTask + +from .compute_anomaly_subtask import ComputeAnomalySubtask +from .plot_hovmoller_subtask import PlotHovmollerSubtask + + +class TimeSeriesTemperatureAnomaly(AnalysisTask): + """ + Performs analysis of time series of temperature anomalies from the first + simulation year as a function of depth. + + Authors + ------- + Xylar Asay-Davis + """ + + def __init__(self, config, mpasTimeSeriesTask): # {{{ + """ + Construct the analysis task. + + Parameters + ---------- + config : instance of MpasAnalysisConfigParser + Contains configuration options + + mpasTimeSeriesTask : ``MpasTimeSeriesTask`` + The task that extracts the time series from MPAS monthly output + + Authors + ------- + Xylar Asay-Davis + """ + # first, call the constructor from the base class (AnalysisTask) + super(TimeSeriesTemperatureAnomaly, self).__init__( + config=config, + taskName='timeSeriesTemperatureAnomaly', + componentName='ocean', + tags=['timeSeries', 'temperature']) + + sectionName = 'hovmollerTemperatureAnomaly' + regionNames = config.getExpression(sectionName, 'regions') + movingAveragePoints = config.getint(sectionName, 'movingAveragePoints') + + mpasFieldName = 'timeMonthly_avg_avgValueWithinOceanLayerRegion_' \ + 'avgLayerTemperature' + + timeSeriesFileName = 'regionAveragedTemperatureAnomaly.nc' + + anomalyTask = ComputeAnomalySubtask( + parentTask=self, + mpasTimeSeriesTask=mpasTimeSeriesTask, + outFileName=timeSeriesFileName, + variableList=[mpasFieldName], + movingAveragePoints=movingAveragePoints) + self.add_subtask(anomalyTask) + + for regionName in regionNames: + caption = 'Trend of {} Temperature Anomaly vs depth'.format( + regionName) + plotTask = PlotHovmollerSubtask( + parentTask=self, + regionName=regionName, + inFileName=timeSeriesFileName, + outFileLabel='TAnomalyZ', + fieldNameInTitle='Temperature Anomaly', + mpasFieldName=mpasFieldName, + unitsLabel='[$^\circ$C]', + sectionName=sectionName, + thumbnailSuffix=u'ΔT', + imageCaption=caption, + galleryGroup='Trends vs Depth', + groupSubtitle=None, + groupLink='trendsvsdepth', + galleryName=None) + + plotTask.run_after(anomalyTask) + self.add_subtask(plotTask) + + # }}} + + # }}} + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/shared/io/__init__.py b/mpas_analysis/shared/io/__init__.py index b03a21757..daa194070 100644 --- a/mpas_analysis/shared/io/__init__.py +++ b/mpas_analysis/shared/io/__init__.py @@ -1,4 +1,4 @@ from .namelist_streams_interface import NameList, StreamsFile from .utility import paths from .write_netcdf import write_netcdf -from .mpas_reader import open_mpas_dataset +from .mpas_reader import open_mpas_dataset, subset_variables diff --git a/mpas_analysis/shared/io/mpas_reader.py b/mpas_analysis/shared/io/mpas_reader.py index 4e03a2589..720d9f2fd 100644 --- a/mpas_analysis/shared/io/mpas_reader.py +++ b/mpas_analysis/shared/io/mpas_reader.py @@ -34,7 +34,9 @@ def open_mpas_dataset(fileName, calendar, timeVariableNames : str or list of 2 str, optional The name of the time variable (typically ``'xtime'`` - or ``['xtime_startMonthly', 'xtime_endMonthly']``) + or ``['xtime_startMonthly', 'xtime_endMonthly']``), or ``None`` if + time does not need to be parsed (and is already in the ``Time`` + variable) variableList : list of strings, optional If present, a list of variables to be included in the data set @@ -63,7 +65,8 @@ def open_mpas_dataset(fileName, calendar, ds = xarray.open_dataset(fileName, decode_cf=True, decode_times=False, lock=False) - ds = _parse_dataset_time(ds, timeVariableNames, calendar) + if timeVariableNames is not None: + ds = _parse_dataset_time(ds, timeVariableNames, calendar) if startDate is not None and endDate is not None: if isinstance(startDate, six.string_types): diff --git a/mpas_analysis/shared/plot/plotting.py b/mpas_analysis/shared/plot/plotting.py index 33d5258b9..91dbeea8d 100644 --- a/mpas_analysis/shared/plot/plotting.py +++ b/mpas_analysis/shared/plot/plotting.py @@ -21,7 +21,7 @@ from matplotlib.ticker import FuncFormatter, FixedLocator import numpy as np from functools import partial -from mpl_toolkits.axes_grid1 import make_axes_locatable, ImageGrid +from mpl_toolkits.axes_grid1 import make_axes_locatable from ..timekeeping.utility import days_to_datetime, date_to_days @@ -115,8 +115,12 @@ def timeseries_analysis_plot(config, dsvalues, N, title, xlabel, ylabel, dsvalue = dsvalues[dsIndex] if dsvalue is None: continue - mean = pd.Series.rolling(dsvalue.to_pandas(), N, center=True).mean() - mean = xr.DataArray.from_series(mean) + if N == 1 or N is None: + mean = dsvalue + else: + mean = pd.Series.rolling(dsvalue.to_pandas(), N, + center=True).mean() + mean = xr.DataArray.from_series(mean) minDays.append(mean.Time.min()) maxDays.append(mean.Time.max()) @@ -126,7 +130,7 @@ def timeseries_analysis_plot(config, dsvalues, N, title, xlabel, ylabel, stride = int(round(nTime/float(maxPoints[dsIndex]))) mean = mean.isel(Time=slice(0, None, stride)) - plt.plot(mean['Time'], mean, + plt.plot(mean['Time'].values, mean.values, lineStyles[dsIndex], linewidth=lineWidths[dsIndex], label=legendText[dsIndex]) diff --git a/mpas_analysis/shared/time_series/__init__.py b/mpas_analysis/shared/time_series/__init__.py index 61f54bd5c..7106aab94 100644 --- a/mpas_analysis/shared/time_series/__init__.py +++ b/mpas_analysis/shared/time_series/__init__.py @@ -1,2 +1,6 @@ from .time_series import cache_time_series from .mpas_time_series_task import MpasTimeSeriesTask + +from .anomaly import compute_moving_avg_anomaly_from_start +from .moving_average import compute_moving_avg + diff --git a/mpas_analysis/shared/time_series/anomaly.py b/mpas_analysis/shared/time_series/anomaly.py new file mode 100644 index 000000000..46a2d4ae8 --- /dev/null +++ b/mpas_analysis/shared/time_series/anomaly.py @@ -0,0 +1,86 @@ +# -*- coding: utf-8 -*- + +from __future__ import absolute_import, division, print_function, \ + unicode_literals + +from ..io import open_mpas_dataset +from .moving_average import compute_moving_avg + + +def compute_moving_avg_anomaly_from_start(timeSeriesFileName, variableList, + simulationStartTime, startDate, + endDate, calendar, + movingAveragePoints=12, + alter_dataset=None): # {{{ + + ''' + Compute the rolling mean of the anomaly of a quantity from the beginning + of the simulation (such that the rolling mean starts at zero by definition) + + Parameters + ---------- + timeSeriesFileName : str + a file produced by ``MpasTimeSeriesTask`` containing variables, the + anomaly and rolling mean of which is to be computed + + variableList : list of str + variable names to include in the resulting data set + + simulationStartTime : str + the start date of the simulation + + startDate, endDate : str + the start and end dates of the time series + + calendar : {'gregorian', 'gregoraian_noleap'} + The calendar used in the MPAS run + + movingAveragePoints : int, optional + The number of points (months) over which to perform the rolling average + of the data set + + alter_dataset : function + A function for manipulating the data set (e.g. computing new + variables), taking an ``xarray.Dataset`` as input argument and + returning an ``xarray.Dataset`` + + Returns + ------- + ds : ``xarray.Dataset`` + The anomaly of the rolling time mean from the start of the simulation + + Authors + ------- + Xylar Asay-Davis + ''' + + ds = open_mpas_dataset(fileName=timeSeriesFileName, + calendar=calendar, + variableList=variableList, + startDate=startDate, + endDate=endDate) + + if alter_dataset is not None: + ds = alter_dataset(ds) + + dsStart = open_mpas_dataset( + fileName=timeSeriesFileName, + calendar=calendar, + variableList=variableList, + startDate=simulationStartTime) + + if alter_dataset is not None: + dsStart = alter_dataset(dsStart) + + dsStart = dsStart.isel(Time=slice(0, movingAveragePoints)).mean('Time') + + for variable in ds.data_vars: + ds[variable] = ds[variable] - dsStart[variable] + + ds = compute_moving_avg(ds) + + return ds + + # }}} + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/shared/time_series/moving_average.py b/mpas_analysis/shared/time_series/moving_average.py new file mode 100644 index 000000000..576c891ba --- /dev/null +++ b/mpas_analysis/shared/time_series/moving_average.py @@ -0,0 +1,38 @@ +# -*- coding: utf-8 -*- + +from __future__ import absolute_import, division, print_function, \ + unicode_literals + + +def compute_moving_avg(ds, movingAveragePoints=12): # {{{ + + ''' + Compute the rolling mean of a data set + + Parameters + ---------- + ds : ``xarray.Dataset`` + a dataset to be averaged + + movingAveragePoints : int, optional + The number of points (months) over which to perform the rolling average + of the data set + + Returns + ------- + ds : ``xarray.Dataset`` + The anomaly of the rolling time mean from the start of the simulation + + Authors + ------- + Xylar Asay-Davis + ''' + + ds = ds.rolling(Time=movingAveragePoints, + center=True).mean().dropna('Time') + + return ds + + # }}} + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/run_mpas_analysis b/run_mpas_analysis index f1b13b590..d3a4fac02 100755 --- a/run_mpas_analysis +++ b/run_mpas_analysis @@ -86,7 +86,13 @@ def build_analysis_list(config): # {{{ analyses.append(oceanTimeSeriesTask) analyses.append(ocean.TimeSeriesAntarcticMelt(config, oceanTimeSeriesTask)) - analyses.append(ocean.TimeSeriesOHC(config, oceanTimeSeriesTask)) + analyses.append(ocean.TimeSeriesTemperatureAnomaly(config, + oceanTimeSeriesTask)) + analyses.append(ocean.TimeSeriesSalinityAnomaly(config, + oceanTimeSeriesTask)) + analyses.append(ocean.TimeSeriesOHCAnomaly(config, + oceanTimeSeriesTask)) + analyses.append(ocean.TimeSeriesSST(config, oceanTimeSeriesTask)) analyses.append(ocean.MeridionalHeatTransport(config, oceanClimatolgyTask)) analyses.append(ocean.StreamfunctionMOC(config, oceanClimatolgyTask))