From 42c5611058c65fea1c41a6649fd26a9aa700c03f Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 9 May 2017 22:50:56 +0200 Subject: [PATCH] Change climatology to classes Change functionality of the climatology module to be handled by 3 new classes: Climatology - the base class MpasClimatology - for computing, caching and/or remapping MPAS climatologies ObservationClimatology - for computing and/or remapping observational climatologies These classes know more about tasks than the previous standalone functions, and can therefore automatically handle much more of the process of setting up climatologies. This makes the tasks themselves shorter and cleaner, and hopefully easier to understand. Functionality has been added to AnalysisTask to cache the times only from a given stream. This is useful for later determining which files need to be opened when computing a given climatology or time series and which have already been cached. Climatology caching in the MpasClimatology class has been updated to use the time cache. Note: the time cache is a dictionary, and thus has been stored in a python ".pickle" file, since NetCDF doesn't support python dictionaries. Since many tasks depend on the results from the cached times, a task has been added to cache these times before calling other tasks. Each analysis task now has a list of prerequisite tasks that must run before it can be run. In addition to supporting caching times, in the future this functionality could be used to build tasks out of multiple subtasks (e.g. one to compute the climatology and one to do the plotting). This would be particularly useful if multiple tasks depend on the same climatology (which is not currently the case). All analysis tasks that use the climatology module have been updated to use the new classes. All tests related to climatologies and the AnalysisTask base class have been updated. The functions "read" and "create" in MeshDescriptors have been made static functions that return an instance of the class, saving the need to call the empty constructor before calling these methods. mpas_xarray has been updated so that a coordinate can be passed to the variable list. This is useful if the only information needed from a data set is the Time coordinate (as is the case for time caching). Cleanup in config.default: * rename regridded to remapped * remove overwrite flags (we assume they are false and the better way to handle this is to delete the cache files) --- config.default | 44 +- mpas_analysis/ocean/climatology_map.py | 303 ++- mpas_analysis/ocean/index_nino34.py | 21 +- .../ocean/meridional_heat_transport.py | 110 +- mpas_analysis/ocean/streamfunction_moc.py | 109 +- mpas_analysis/ocean/time_series_ohc.py | 2 +- mpas_analysis/ocean/time_series_sst.py | 2 +- mpas_analysis/sea_ice/climatology_map.py | 411 ++-- .../sea_ice/sea_ice_analysis_task.py | 3 +- mpas_analysis/sea_ice/time_series.py | 2 +- mpas_analysis/shared/analysis_task.py | 410 +++- .../shared/cache_dataset_times_task.py | 153 ++ mpas_analysis/shared/climatology/__init__.py | 7 +- .../shared/climatology/climatology.py | 1826 ++++++++++------- .../shared/generalized_reader/__init__.py | 1 + .../generalized_reader/generalized_reader.py | 2 +- mpas_analysis/shared/grid/grid.py | 189 +- .../shared/interpolation/remapper.py | 47 - mpas_analysis/shared/io/__init__.py | 3 +- mpas_analysis/shared/io/utility.py | 16 +- .../shared/mpas_xarray/mpas_xarray.py | 21 +- mpas_analysis/shared/timekeeping/utility.py | 104 +- .../ocean_maps.py | 4 +- mpas_analysis/test/test_analysis_task | 1 + mpas_analysis/test/test_analysis_task.py | 121 +- mpas_analysis/test/test_climatology.py | 340 ++- mpas_analysis/test/test_climatology/mpas-o_in | 1083 ++++++++++ .../test/test_climatology/streams.ocean | 16 + mpas_analysis/test/test_interpolate.py | 12 +- run_analysis.py | 422 ++-- 30 files changed, 3826 insertions(+), 1959 deletions(-) create mode 100644 mpas_analysis/shared/cache_dataset_times_task.py create mode 120000 mpas_analysis/test/test_analysis_task create mode 100644 mpas_analysis/test/test_climatology/mpas-o_in create mode 100644 mpas_analysis/test/test_climatology/streams.ocean diff --git a/config.default b/config.default index 5005b627d..08c539801 100644 --- a/config.default +++ b/config.default @@ -41,7 +41,7 @@ parallelTaskCount = 1 # Prefix on the commnd line before a parallel task (e.g. 'srun -n 1 python') # Default is no prefix (run_analysis.py is executed directly) -commandPrefix = +commandPrefix = [input] ## options related to reading in the results to be analyzed @@ -99,9 +99,10 @@ scratchSubdirectory = scratch plotsSubdirectory = plots logsSubdirectory = logs mpasClimatologySubdirectory = clim/mpas -mpasRegriddedClimSubdirectory = clim/mpas/regridded +mpasRemappedClimSubdirectory = clim/mpas/remapped mappingSubdirectory = mapping timeSeriesSubdirectory = timeseries +timeCacheSubdirectory = timecache # a list of analyses to generate. Valid names are: # 'timeSeriesOHC', 'timeSeriesSST', 'climatologyMapSST', @@ -153,9 +154,6 @@ comparisonLonResolution = 0.5 # the interpolation method # mpasMappingFile = /path/to/mapping/file -# overwrite files when building climatologies? -overwriteMpasClimatology = False - # interpolation order for model and observation results. Likely values are # 'bilinear', 'neareststod' (nearest neighbor) or 'conserve' mpasInterpolationMethod = bilinear @@ -230,14 +228,11 @@ interpolationMethod = bilinear # The directories where observation climatologies will be stored if they need # to be computed. If a relative path is supplied, it is relative to the output # base directory. If an absolute path is supplied, this should point to -# cached climatology files on the desired comparison grid, in which case -# overwriteObsClimatology should be False. If cached regridded files are -# supplied, there is no need to provide cached files before regridding. +# cached climatology files on the desired comparison grid. If cached remapped +# files are supplied, there is no need to provide cached files before +# remapping. climatologySubdirectory = clim/obs -regriddedClimSubdirectory = clim/obs/regridded - -# overwrite files when building climatologies? -overwriteObsClimatology = False +remappedClimSubdirectory = clim/obs/remapped [oceanReference] ## options related to ocean reference run with which the results will be @@ -252,11 +247,11 @@ baseDirectory = /dir/to/ocean/reference # directory where ocean reference simulation results are stored baseDirectory = /dir/to/ocean/reference - + [seaIceObservations] ## options related to sea ice observations with which the results will be ## compared - + # directory where sea ice observations are stored baseDirectory = /dir/to/seaice/observations areaNH = IceArea_timeseries/iceAreaNH_climo.nc @@ -290,14 +285,11 @@ interpolationMethod = bilinear # The directories where observation climatologies will be stored if they need # to be computed. If a relative path is supplied, it is relative to the output # base directory. If an absolute path is supplied, this should point to -# cached climatology files on the desired comparison grid, in which case -# overwriteObsClimatology should be False. If cached regridded files are -# supplied, there is no need to provide cached files before regridding. +# cached climatology files on the desired comparison grid. If cached remapped +# files are supplied, there is no need to provide cached files before +# remapping. climatologySubdirectory = clim/obs -regriddedClimSubdirectory = clim/obs/regridded - -# overwrite files when building climatologies? -overwriteObsClimatology = False +remappedClimSubdirectory = clim/obs/remapped [seaIceReference] ## options related to sea ice reference run with which the results will be @@ -399,7 +391,7 @@ regionMaskFiles = /path/to/MOCregional/mapping/file # is handled automatically. If the MOC calculation encounters memory problems, # consider setting maxChunkSize to a number significantly lower than nEdges # in your MPAS mesh so that the calculation will be divided into smaller -# pieces. +# pieces. # Note, need to use maxChunkSize for the 18to6 # maxChunkSize = 1000 @@ -450,7 +442,7 @@ titleFontSize = 18 polarPlot = False [climatologyMapSST] -## options related to plotting horizontally regridded climatologies of +## options related to plotting horizontally remapped climatologies of ## sea surface temperature (SST) against reference model results and ## observations @@ -473,7 +465,7 @@ colorbarLevelsDifference = [-5, -3, -2, -1, 0, 1, 2, 3, 5] comparisonTimes = ['JFM', 'JAS', 'ANN'] [climatologyMapSSS] -## options related to plotting horizontally regridded climatologies of +## options related to plotting horizontally remapped climatologies of ## sea surface salinity (SSS) against reference model results and observations # colormap for model/observations @@ -495,7 +487,7 @@ colorbarLevelsDifference = [-3, -2, -1, -0.5, 0, 0.5, 1, 2, 3] comparisonTimes = ['JFM', 'JAS', 'ANN'] [climatologyMapMLD] -## options related to plotting horizontally regridded climatologies of +## options related to plotting horizontally remapped climatologies of ## mixed layer depth (MLD) against reference model results and observations # colormap for model/observations @@ -517,7 +509,7 @@ colorbarLevelsDifference = [-150, -80, -30, -10, 0, 10, 30, 80, 150] comparisonTimes = ['JFM', 'JAS', 'ANN'] [climatologyMapSeaIceConcThick] -## options related to plotting horizontally regridded climatologies of +## options related to plotting horizontally remapped climatologies of ## sea ice concentration and thickness against reference model results and ## observations diff --git a/mpas_analysis/ocean/climatology_map.py b/mpas_analysis/ocean/climatology_map.py index a70d768fd..1f00ab572 100644 --- a/mpas_analysis/ocean/climatology_map.py +++ b/mpas_analysis/ocean/climatology_map.py @@ -11,8 +11,6 @@ import xarray as xr import datetime import numpy as np -import os -import warnings from ..shared.analysis_task import AnalysisTask @@ -20,25 +18,18 @@ setup_colormap from ..shared.constants import constants -from ..shared.io.utility import build_config_full_path +from ..shared.io import build_config_full_path from ..shared.generalized_reader.generalized_reader \ import open_multifile_dataset from ..shared.timekeeping.utility import get_simulation_start_time -from ..shared.climatology import get_lat_lon_comparison_descriptor, \ - get_remapper, get_mpas_climatology_file_names, \ - get_observation_climatology_file_names, \ - compute_climatology, cache_climatologies, update_start_end_year, \ - remap_and_write_climatology - -from ..shared.grid import MpasMeshDescriptor, LatLonGridDescriptor +from ..shared.grid import LatLonGridDescriptor +from ..shared.climatology import MpasClimatology, ObservationClimatology from ..shared.mpas_xarray import mpas_xarray -from ..shared.interpolation import Remapper - class ClimatologyMapOcean(AnalysisTask): # {{{ """ @@ -49,6 +40,46 @@ class ClimatologyMapOcean(AnalysisTask): # {{{ Luke Van Roekel, Xylar Asay-Davis, Milena Veneziani """ + def __init__(self, config, taskName, componentName, tags=[]): # {{{ + """ + Construct the analysis task. + + Parameters + ---------- + config : instance of MpasAnalysisConfigParser + Contains configuration options + + taskName : str + The name of the task, typically the same as the class name except + starting with lowercase (e.g. 'myTask' for class 'MyTask') + + componentName : {'ocean', 'seaIce'} + The name of the component (same as the folder where the task + resides) + + tags : list of str, optional + Tags used to describe the task (e.g. 'timeSeries', 'climatology', + horizontalMap', 'index', 'transect'). These are used to determine + which tasks are generated (e.g. 'all_transect' or 'no_climatology' + in the 'generate' flags) + + Authors + ------- + Xylar Asay-Davis + """ + # call the constructor from the base class (AnalysisTask) + super(ClimatologyMapOcean, self).__init__( + config, taskName, componentName, tags, + prerequisiteTasks=['cacheOceanTimeSeriesStatsTimes']) + + # by default, we don't override "useNcremap" from the config file. + # However, some child classes will set this explicitly, since results + # are better for remapping observations with the masking provided in + # MPAS-Analysis directly + self.useNcremapObs = None + + # }}} + def setup_and_check(self): # {{{ """ Perform steps to set up the analysis and check for errors in the setup. @@ -68,6 +99,8 @@ def setup_and_check(self): # {{{ analysisOptionName='config_am_timeseriesstatsmonthly_enable', raiseException=True) + self.simulationStartTime = get_simulation_start_time(self.runStreams) + # }}} def run(self): # {{{ @@ -89,34 +122,8 @@ def run(self): # {{{ # get local versions of member variables for convenience config = self.config - calendar = self.calendar fieldName = self.fieldName - simulationStartTime = get_simulation_start_time(self.runStreams) - - # get a list of timeSeriesStats output files from the streams file, - # reading only those that are between the start and end dates - startDate = config.get('climatology', 'startDate') - endDate = config.get('climatology', 'endDate') - streamName = \ - self.historyStreams.find_stream(self.streamMap['timeSeriesStats']) - inputFiles = self.historyStreams.readpath(streamName, - startDate=startDate, - endDate=endDate, - calendar=calendar) - print '\n Reading files:\n' \ - ' {} through\n {}'.format( - os.path.basename(inputFiles[0]), - os.path.basename(inputFiles[-1])) - - mainRunName = config.get('runs', 'mainRunName') - - overwriteMpasClimatology = config.getWithDefault( - 'climatology', 'overwriteMpasClimatology', False) - - overwriteObsClimatology = config.getWithDefault( - 'oceanObservations', 'overwriteObsClimatology', False) - try: restartFileName = self.runStreams.readpath('restart')[0] except ValueError: @@ -125,27 +132,6 @@ def run(self): # {{{ outputTimes = config.getExpression(self.taskName, 'comparisonTimes') - comparisonDescriptor = get_lat_lon_comparison_descriptor(config) - - varList = [fieldName] - - ds = open_multifile_dataset(fileNames=inputFiles, - calendar=calendar, - config=config, - simulationStartTime=simulationStartTime, - timeVariableName='Time', - variableList=varList, - iselValues=self.iselValues, - variableMap=self.variableMap, - startDate=startDate, - endDate=endDate) - - changed, startYear, endYear = update_start_end_year(ds, config, - calendar) - - mpasDescriptor = MpasMeshDescriptor( - restartFileName, meshName=config.get('input', 'mpasMeshName')) - parallel = self.config.getint('execute', 'parallelTaskCount') > 1 if parallel: # avoid writing the same mapping file from multiple processes @@ -153,122 +139,78 @@ def run(self): # {{{ else: mappingFilePrefix = 'map' - mpasRemapper = get_remapper( - config=config, sourceDescriptor=mpasDescriptor, - comparisonDescriptor=comparisonDescriptor, - mappingFileSection='climatology', - mappingFileOption='mpasMappingFile', - mappingFilePrefix=mappingFilePrefix, - method=config.get('climatology', 'mpasInterpolationMethod')) - - obsDescriptor = LatLonGridDescriptor() - obsDescriptor.read(fileName=self.obsFileName, latVarName='lat', - lonVarName='lon') - - origObsRemapper = Remapper(comparisonDescriptor, obsDescriptor) - (colormapResult, colorbarLevelsResult) = setup_colormap( config, self.taskName, suffix='Result') (colormapDifference, colorbarLevelsDifference) = setup_colormap( config, self.taskName, suffix='Difference') - dsObs = None - obsRemapperBuilt = False + # we don't have any way to know the observation lat/lon and grid + # without this, so we have to read the observational data set even + # if it's already been remapped. + dsObs = self._build_observational_dataset() + # create a descriptor of the observation grid using the lat/lon + # coordinates + obsDescriptor = LatLonGridDescriptor.read(ds=dsObs) # Interpolate and compute biases for months in outputTimes: monthValues = constants.monthDictionary[months] - (climatologyFileName, climatologyPrefix, regriddedFileName) = \ - get_mpas_climatology_file_names( - config=config, - fieldName=fieldName, - monthNames=months, - mpasMeshName=mpasDescriptor.meshName, - comparisonGridName=comparisonDescriptor.meshName) - - if (overwriteMpasClimatology or - not os.path.exists(regriddedFileName)): - seasonalClimatology = cache_climatologies( - ds, monthValues, config, climatologyPrefix, calendar, + mpasClimatology = MpasClimatology( + task=self, + fieldName=fieldName, + monthNames=months, + streamName='timeSeriesStats', + meshFileName=restartFileName, + comparisonGrid='latlon', + mappingFileSection='climatology', + mappingFileOption='mpasMappingFile', + mappingFilePrefix=mappingFilePrefix, + method=config.get('climatology', 'mpasInterpolationMethod')) + + if mpasClimatology.remappedDataSet is None: + # the remapped climatology hasn't been cached yet + mpasClimatology.cache( + openDataSetFunc=self._open_mpas_dataset_part, printProgress=True) + mpasClimatology.remap_and_write() - if seasonalClimatology is None: - # apparently, there was no data available to create the - # climatology - warnings.warn('no data to create {} climatology for ' - '{}'.format(fieldName, months)) - continue - - remappedClimatology = remap_and_write_climatology( - config, seasonalClimatology, climatologyFileName, - regriddedFileName, mpasRemapper) - - else: - - remappedClimatology = xr.open_dataset(regriddedFileName) - - modelOutput = remappedClimatology[fieldName].values - lon = remappedClimatology['lon'].values - lat = remappedClimatology['lat'].values + modelOutput = \ + mpasClimatology.remappedDataSet[self.mpasFieldName].values + lon = mpasClimatology.remappedDataSet['lon'].values + lat = mpasClimatology.remappedDataSet['lat'].values lonTarg, latTarg = np.meshgrid(lon, lat) - # now the observations - (climatologyFileName, regriddedFileName) = \ - get_observation_climatology_file_names( - config=config, fieldName=fieldName, monthNames=months, - componentName='ocean', remapper=origObsRemapper) - - if (overwriteObsClimatology or - not os.path.exists(regriddedFileName)): - - if dsObs is None: - # load the observations the first time - dsObs = self._build_observational_dataset() - - seasonalClimatology = compute_climatology( - dsObs, monthValues, maskVaries=True) - - if not obsRemapperBuilt: - seasonalClimatology.load() - seasonalClimatology.close() - seasonalClimatology.to_netcdf(climatologyFileName) - # make the remapper for the climatology - obsDescriptor = LatLonGridDescriptor() - obsDescriptor.read(fileName=climatologyFileName, - latVarName='lat', - lonVarName='lon') - - obsRemapper = get_remapper( - config=config, sourceDescriptor=obsDescriptor, - comparisonDescriptor=comparisonDescriptor, - mappingFileSection='oceanObservations', - mappingFileOption='{}ClimatologyMappingFile'.format( - fieldName), - mappingFilePrefix='map_obs_{}'.format(fieldName), - method=config.get('oceanObservations', - 'interpolationMethod')) - - obsRemapperBuilt = True - - if obsRemapper is None: - # no need to remap because the observations are on the - # comparison grid already - remappedClimatology = seasonalClimatology - else: - remappedClimatology = \ - remap_and_write_climatology( - config, seasonalClimatology, climatologyFileName, - regriddedFileName, obsRemapper) - - else: - - remappedClimatology = xr.open_dataset(regriddedFileName) - observations = remappedClimatology[self.obsFieldName].values + obsClimatology = \ + ObservationClimatology( + task=self, + fieldName=self.obsFieldName, + monthNames=months, + obsGridDescriptor=obsDescriptor, + comparisonGrid='latlon', + mappingFileSection='oceanObservations', + mappingFileOption='{}ClimatologyMappingFile'.format( + fieldName), + mappingFilePrefix='map_obs_{}'.format(fieldName), + method=config.get('oceanObservations', + 'interpolationMethod')) + + if obsClimatology.remappedDataSet is None: + # the remapped climatology hasn't been cached yet + obsClimatology.compute(ds=dsObs, monthValues=monthValues, + maskVaries=True) + obsClimatology.remap_and_write(useNcremap=self.useNcremapObs) + + observations = \ + obsClimatology.remappedDataSet[self.obsFieldName].values bias = modelOutput - observations + startYear = mpasClimatology.startYear + endYear = mpasClimatology.endYear + + mainRunName = config.get('runs', 'mainRunName') outFileName = '{}/{}_{}_{}_years{:04d}-{:04d}.png'.format( self.plotsDirectory, self.outFileLabel, mainRunName, months, startYear, endYear) @@ -293,6 +235,40 @@ def run(self): # {{{ # }}} + def _open_mpas_dataset_part(self, inputFileNames, startDate, + endDate): # {{{ + """ + Open part of a data set between the given start and end date, used + to cache a climatology of the data set. + + Parameters + ---------- + inputFileNames : list of str + File names in the multifile data set to open + + startDate, endDate : float + start and end date to which to crop the Time dimension (given in + days since 0001-01-01) + + Authors + ------- + Xylar Asay-Davis + """ + varList = [self.mpasFieldName] + + ds = open_multifile_dataset( + fileNames=inputFileNames, + calendar=self.calendar, + config=self.config, + simulationStartTime=self.simulationStartTime, + timeVariableName='Time', + variableList=varList, + iselValues=self.iselValues, + variableMap=self.variableMap, + startDate=startDate, + endDate=endDate) + return ds # }}} + # }}} @@ -343,7 +319,7 @@ def setup_and_check(self): # {{{ # self.runDirectory , self.historyDirectory, self.plotsDirectory, # self.namelist, self.runStreams, self.historyStreams, # self.calendar, self.namelistMap, self.streamMap, self.variableMap - super(ClimatologyMapOcean, self).setup_and_check() + super(ClimatologyMapSST, self).setup_and_check() observationsDirectory = build_config_full_path( self.config, 'oceanObservations', @@ -354,6 +330,7 @@ def setup_and_check(self): # {{{ observationsDirectory) self.iselValues = {'nVertLevels': 0} + self.mpasFieldName = 'temperature' self.obsFieldName = 'SST' @@ -452,7 +429,7 @@ def setup_and_check(self): # {{{ # self.runDirectory , self.historyDirectory, self.plotsDirectory, # self.namelist, self.runStreams, self.historyStreams, # self.calendar, self.namelistMap, self.streamMap, self.variableMap - super(ClimatologyMapOcean, self).setup_and_check() + super(ClimatologyMapSSS, self).setup_and_check() observationsDirectory = build_config_full_path( self.config, 'oceanObservations', @@ -463,6 +440,7 @@ def setup_and_check(self): # {{{ observationsDirectory) self.iselValues = {'nVertLevels': 0} + self.mpasFieldName = 'salinity' self.obsFieldName = 'SSS' @@ -519,6 +497,10 @@ def __init__(self, config): # {{{ Xylar Asay-Davis """ + # We want to use the "online" remapper because it does a better job of + # masking missing values for MLD. + self.useNcremapObs = False + self.fieldName = 'mld' self.fieldNameInTitle = 'MLD' @@ -544,7 +526,7 @@ def setup_and_check(self): # {{{ # self.runDirectory , self.historyDirectory, self.plotsDirectory, # self.namelist, self.runStreams, self.historyStreams, # self.calendar, self.namelistMap, self.streamMap, self.variableMap - super(ClimatologyMapOcean, self).setup_and_check() + super(ClimatologyMapMLD, self).setup_and_check() observationsDirectory = build_config_full_path( self.config, 'oceanObservations', @@ -555,6 +537,7 @@ def setup_and_check(self): # {{{ observationsDirectory) self.iselValues = None + self.mpasFieldName = 'mld' self.obsFieldName = 'mld_dt_mean' diff --git a/mpas_analysis/ocean/index_nino34.py b/mpas_analysis/ocean/index_nino34.py index f7b6ea74c..6b48c198d 100644 --- a/mpas_analysis/ocean/index_nino34.py +++ b/mpas_analysis/ocean/index_nino34.py @@ -5,13 +5,14 @@ from scipy import signal, stats import os -from ..shared.climatology import climatology +from ..shared.climatology import Climatology from ..shared.constants import constants -from ..shared.io.utility import build_config_full_path +from ..shared.io import build_config_full_path from ..shared.generalized_reader.generalized_reader \ import open_multifile_dataset -from ..shared.timekeeping.utility import get_simulation_start_time +from ..shared.timekeeping.utility import get_simulation_start_time, \ + add_years_months_days_in_month from ..shared.plot.plotting import nino34_timeseries_plot, nino34_spectra_plot @@ -139,7 +140,7 @@ def run(self): # {{{ print ' Compute NINO3.4 index...' regionSST = ds.avgSurfaceTemperature.isel(nOceanRegions=regionIndex) - nino34 = self._compute_nino34_index(regionSST, calendar) + nino34 = self._compute_nino34_index(regionSST) # Compute the observational index over the entire time range # nino34Obs = compute_nino34_index(dsObs.sst, calendar) @@ -183,7 +184,7 @@ def run(self): # {{{ obsTitle, figureName, linewidths=2) # }}} - def _compute_nino34_index(self, regionSST, calendar): # {{{ + def _compute_nino34_index(self, regionSST): # {{{ """ Computes nino34 index time series. It follow the standard nino34 algorithm, i.e., @@ -217,15 +218,13 @@ def _compute_nino34_index(self, regionSST, calendar): # {{{ raise ValueError('regionSST should be an xarray DataArray') # add 'month' data array so we can group by month below. - regionSST = climatology.add_years_months_days_in_month(regionSST, - calendar) + regionSST = add_years_months_days_in_month(regionSST, self.calendar) # Compute monthly average and anomaly of climatology of SST - monthlyClimatology = \ - climatology.compute_monthly_climatology(regionSST, - maskVaries=False) + monthlyClimatology = Climatology(task=self) + monthlyClimatology.compute_monthly(regionSST, maskVaries=False) - anomaly = regionSST.groupby('month') - monthlyClimatology + anomaly = regionSST.groupby('month') - monthlyClimatology.dataSet # Remove the long term trend from the anomalies detrendedAnomal = signal.detrend(anomaly.values) diff --git a/mpas_analysis/ocean/meridional_heat_transport.py b/mpas_analysis/ocean/meridional_heat_transport.py index 160c13e35..0a799f522 100644 --- a/mpas_analysis/ocean/meridional_heat_transport.py +++ b/mpas_analysis/ocean/meridional_heat_transport.py @@ -8,14 +8,14 @@ from ..shared.plot.plotting import plot_vertical_section,\ setup_colormap, plot_1D -from ..shared.io.utility import build_config_full_path, make_directories +from ..shared.io import build_config_full_path, make_directories from ..shared.generalized_reader.generalized_reader \ import open_multifile_dataset from ..shared.timekeeping.utility import get_simulation_start_time -from ..shared.climatology.climatology import cache_climatologies +from ..shared.climatology import MpasClimatology from ..shared.analysis_task import AnalysisTask @@ -45,10 +45,11 @@ def __init__(self, config): # {{{ ''' # first, call the constructor from the base class (AnalysisTask) super(MeridionalHeatTransport, self).__init__( - config=config, - taskName='meridionalHeatTransport', - componentName='ocean', - tags=['climatology']) + config=config, + taskName='meridionalHeatTransport', + componentName='ocean', + tags=['climatology'], + prerequisiteTasks=['cacheOceanTimeSeriesStatsTimes']) # }}} def setup_and_check(self): # {{{ @@ -80,23 +81,8 @@ def setup_and_check(self): # {{{ analysisOptionName='config_am_meridionalheattransport_enable', raiseException=True) - # Get a list of timeSeriesStats output files from the streams file, - # reading only those that are between the start and end dates - # First a list necessary for theMHT climatology - streamName = self.historyStreams.find_stream( - self.streamMap['timeSeriesStats']) - self.startDate = config.get('climatology', 'startDate') - self.endDate = config.get('climatology', 'endDate') - self.inputFiles = \ - self.historyStreams.readpath(streamName, - startDate=self.startDate, - endDate=self.endDate, - calendar=self.calendar) self.simulationStartTime = get_simulation_start_time(self.runStreams) - self.startYear = config.getint('climatology', 'startYear') - self.endYear = config.getint('climatology', 'endYear') - self.sectionName = 'meridionalHeatTransport' # Read in obs file information @@ -168,39 +154,16 @@ def run(self): # {{{ # Then we will need to add another section for regions with a loop # over number of regions. ###################################################################### - variableList = ['avgMeridionalHeatTransportLat', - 'avgMeridionalHeatTransportLatZ'] print '\n Compute and plot global meridional heat transport' - outputDirectory = build_config_full_path(config, 'output', - 'mpasClimatologySubdirectory') - - print '\n List of files for climatologies:\n' \ - ' {} through\n {}'.format( - os.path.basename(self.inputFiles[0]), - os.path.basename(self.inputFiles[-1])) - - make_directories(outputDirectory) - - print ' Load data...' - ds = open_multifile_dataset( - fileNames=self.inputFiles, - calendar=self.calendar, - config=config, - simulationStartTime=self.simulationStartTime, - timeVariableName='Time', - variableList=variableList, - variableMap=self.variableMap, - startDate=self.startDate, - endDate=self.endDate) + annualClimatology = MpasClimatology(task=self, + fieldName='mht', + monthNames='ANN', + streamName='timeSeriesStats') - # Compute annual climatology - cachePrefix = '{}/meridionalHeatTransport'.format(outputDirectory) - annualClimatology = cache_climatologies(ds, monthDictionary['ANN'], - config, cachePrefix, - self.calendar, - printProgress=True) + annualClimatology.cache(openDataSetFunc=self._open_mht_part, + printProgress=True) # **** Plot MHT **** # Define plotting variables @@ -212,14 +175,15 @@ def run(self): # {{{ print ' Plot global MHT...' # Plot 1D MHT (zonally averaged, depth integrated) x = binBoundaryMerHeatTrans - y = annualClimatology.avgMeridionalHeatTransportLat + y = annualClimatology.dataSet.avgMeridionalHeatTransportLat xLabel = 'latitude [deg]' yLabel = 'meridional heat transport [PW]' title = 'Global MHT (ANN, years {:04d}-{:04d})\n {}'.format( - self.startYear, self.endYear, mainRunName) + annualClimatology.startYear, annualClimatology.endYear, + mainRunName) figureName = '{}/mht_{}_years{:04d}-{:04d}.png'.format( self.plotsDirectory, mainRunName, - self.startYear, self.endYear) + annualClimatology.startYear, annualClimatology.endYear) if self.observationsFile is not None: # Load in observations dsObs = xr.open_dataset(self.observationsFile) @@ -251,7 +215,7 @@ def run(self): # {{{ # normalize 2D MHT by layer thickness MHTLatZ = \ - annualClimatology.avgMeridionalHeatTransportLatZ.values.T[:, :] + annualClimatology.dataSet.avgMeridionalHeatTransportLatZ.values.T for k in range(nVertLevels): MHTLatZ[k, :] = MHTLatZ[k, :]/refLayerThickness[k] @@ -261,10 +225,11 @@ def run(self): # {{{ xLabel = 'latitude [deg]' yLabel = 'depth [m]' title = 'Global MHT (ANN, years {:04d}-{:04d})\n {}'.format( - self.startYear, self.endYear, mainRunName) + annualClimatology.startYear, annualClimatology.endYear, + mainRunName) figureName = '{}/mhtZ_{}_years{:04d}-{:04d}.png'.format( self.plotsDirectory, mainRunName, - self.startYear, self.endYear) + annualClimatology.startYear, annualClimatology.endYear) colorbarLabel = '[PW/m]' contourLevels = config.getExpression(self.sectionName, 'contourLevelsGlobal', @@ -280,6 +245,39 @@ def run(self): # {{{ invertYAxis=False) # }}} + def _open_mht_part(self, inputFileNames, startDate, endDate): # {{{ + """ + Open part of the MHT data set between the given start and end date, + used to cache a climatology of the data set. + + Parameters + ---------- + inputFileNames : list of str + File names in the multifile data set to open + + startDate, endDate : float + start and end date to which to crop the Time dimension (given in + days since 0001-01-01) + + Authors + ------- + Xylar Asay-Davis + """ + variableList = ['avgMeridionalHeatTransportLat', + 'avgMeridionalHeatTransportLatZ'] + + ds = open_multifile_dataset( + fileNames=inputFileNames, + calendar=self.calendar, + config=self.config, + simulationStartTime=self.simulationStartTime, + timeVariableName='Time', + variableList=variableList, + variableMap=self.variableMap, + startDate=startDate, + endDate=endDate) + + return ds # }}} # }}} # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/streamfunction_moc.py b/mpas_analysis/ocean/streamfunction_moc.py index a04cf5a53..61020eb88 100644 --- a/mpas_analysis/ocean/streamfunction_moc.py +++ b/mpas_analysis/ocean/streamfunction_moc.py @@ -9,7 +9,7 @@ from ..shared.plot.plotting import plot_vertical_section,\ timeseries_analysis_plot, setup_colormap -from ..shared.io.utility import build_config_full_path, make_directories +from ..shared.io import build_config_full_path, make_directories from ..shared.generalized_reader.generalized_reader \ import open_multifile_dataset @@ -17,8 +17,7 @@ from ..shared.timekeeping.utility import get_simulation_start_time, \ days_to_datetime -from ..shared.climatology.climatology import update_start_end_year, \ - cache_climatologies +from ..shared.climatology import MpasClimatology from ..shared.analysis_task import AnalysisTask @@ -55,10 +54,11 @@ def __init__(self, config): # {{{ ''' # first, call the constructor from the base class (AnalysisTask) super(StreamfunctionMOC, self).__init__( - config=config, - taskName='streamfunctionMOC', - componentName='ocean', - tags=['streamfunction', 'moc', 'climatology', 'timeSeries']) + config=config, + taskName='streamfunctionMOC', + componentName='ocean', + tags=['streamfunction', 'moc', 'climatology', 'timeSeries'], + prerequisiteTasks=['cacheOceanTimeSeriesStatsTimes']) # }}} @@ -98,18 +98,8 @@ def setup_and_check(self): # {{{ # First a list necessary for the streamfunctionMOC climatology streamName = self.historyStreams.find_stream( self.streamMap['timeSeriesStats']) - self.startDateClimo = config.get('climatology', 'startDate') - self.endDateClimo = config.get('climatology', 'endDate') - self.inputFilesClimo = \ - self.historyStreams.readpath(streamName, - startDate=self.startDateClimo, - endDate=self.endDateClimo, - calendar=self.calendar) self.simulationStartTime = get_simulation_start_time(self.runStreams) - self.startYearClimo = config.getint('climatology', 'startYear') - self.endYearClimo = config.getint('climatology', 'endYear') - # Then a list necessary for the streamfunctionMOC Atlantic timeseries self.startDateTseries = config.get('timeSeries', 'startDate') self.endDateTseries = config.get('timeSeries', 'endDate') @@ -141,11 +131,6 @@ def run(self): # {{{ print "\n Plotting streamfunction of Meridional Overturning " \ "Circulation (MOC)..." - print '\n List of files for climatologies:\n' \ - ' {} through\n {}'.format( - os.path.basename(self.inputFilesClimo[0]), - os.path.basename(self.inputFilesClimo[-1])) - print '\n List of files for time series:\n' \ ' {} through\n {}'.format( os.path.basename(self.inputFilesTseries[0]), @@ -163,8 +148,10 @@ def run(self): # {{{ # sectionName, dictClimo, # dictTseries) else: - self._cache_velocity_climatologies() - self._compute_moc_climo_postprocess() + velocityClimatology = self._cache_velocity_climatologies() + self.startYearClimo = velocityClimatology.startYear + self.endYearClimo = velocityClimatology.endYear + self._compute_moc_climo_postprocess(velocityClimatology) dsMOCTimeSeries = self._compute_moc_time_series_postprocess() # **** Plot MOC **** @@ -244,42 +231,56 @@ def _load_mesh(self): # {{{ def _cache_velocity_climatologies(self): # {{{ '''compute yearly velocity climatologies and cache them''' - variableList = ['avgNormalVelocity', - 'avgVertVelocityTop'] + velocityClimatology = MpasClimatology( + task=self, + fieldName='meanVelocity', + monthNames='ANN', + streamName='timeSeriesStats') - config = self.config + # compute and cache the velocity climatology + velocityClimatology.cache(openDataSetFunc=self._open_velcoity_part, + printProgress=True) - outputDirectory = build_config_full_path(config, 'output', - 'mpasClimatologySubdirectory') + return velocityClimatology # }}} - make_directories(outputDirectory) + def _open_velcoity_part(self, inputFileNames, startDate, endDate): # {{{ + """ + Open part of the monthly mean velocity data set between the given + start and end date, used to cache a climatology of the data set. - chunking = config.getExpression(self.sectionName, 'maxChunkSize') + Parameters + ---------- + inputFileNames : list of str + File names in the multifile data set to open + + startDate, endDate : float + start and end date to which to crop the Time dimension (given in + days since 0001-01-01) + + Authors + ------- + Xylar Asay-Davis + """ + + variableList = ['avgNormalVelocity', + 'avgVertVelocityTop'] + + chunking = self.config.getExpression(self.sectionName, 'maxChunkSize') ds = open_multifile_dataset( - fileNames=self.inputFilesClimo, + fileNames=inputFileNames, calendar=self.calendar, - config=config, + config=self.config, simulationStartTime=self.simulationStartTime, timeVariableName='Time', variableList=variableList, variableMap=self.variableMap, - startDate=self.startDateClimo, - endDate=self.endDateClimo, + startDate=startDate, + endDate=endDate, chunking=chunking) - # update the start and end year in config based on the real extend of - # ds - update_start_end_year(ds, config, self.calendar) - - cachePrefix = '{}/meanVelocity'.format(outputDirectory) - - # compute and cache the velocity climatology - cache_climatologies(ds, monthDictionary['ANN'], - config, cachePrefix, self.calendar, - printProgress=True) - # }}} + return ds # }}} - def _compute_moc_climo_postprocess(self): # {{{ + def _compute_moc_climo_postprocess(self, velocityClimatology): # {{{ '''compute mean MOC streamfunction as a post-process''' @@ -340,19 +341,7 @@ def _compute_moc_climo_postprocess(self): # {{{ outputDirectory, self.startYearClimo, self.endYearClimo) if not os.path.exists(outputFileClimo): - print ' Load data...' - - cachePrefix = '{}/meanVelocity'.format(outputDirectory) - - if self.startYearClimo == self.endYearClimo: - yearString = '{:04d}'.format(self.startYearClimo) - velClimoFile = '{}_year{}.nc'.format(cachePrefix, yearString) - else: - yearString = '{:04d}-{:04d}'.format(self.startYearClimo, - self.endYearClimo) - velClimoFile = '{}_years{}.nc'.format(cachePrefix, yearString) - - annualClimatology = xr.open_dataset(velClimoFile) + annualClimatology = velocityClimatology.dataSet # Convert to numpy arrays # (can result in a memory error for large array size) diff --git a/mpas_analysis/ocean/time_series_ohc.py b/mpas_analysis/ocean/time_series_ohc.py index ba5c1b883..4c7685539 100644 --- a/mpas_analysis/ocean/time_series_ohc.py +++ b/mpas_analysis/ocean/time_series_ohc.py @@ -14,7 +14,7 @@ from ..shared.time_series import time_series -from ..shared.io.utility import build_config_full_path, make_directories, \ +from ..shared.io import build_config_full_path, make_directories, \ check_path_exists diff --git a/mpas_analysis/ocean/time_series_sst.py b/mpas_analysis/ocean/time_series_sst.py index d1e6e240d..c52fe758d 100644 --- a/mpas_analysis/ocean/time_series_sst.py +++ b/mpas_analysis/ocean/time_series_sst.py @@ -12,7 +12,7 @@ from ..shared.time_series import time_series -from ..shared.io.utility import build_config_full_path, make_directories, \ +from ..shared.io import build_config_full_path, make_directories, \ check_path_exists diff --git a/mpas_analysis/sea_ice/climatology_map.py b/mpas_analysis/sea_ice/climatology_map.py index 91c75702d..e443f96c4 100644 --- a/mpas_analysis/sea_ice/climatology_map.py +++ b/mpas_analysis/sea_ice/climatology_map.py @@ -9,17 +9,13 @@ from ..shared.constants import constants -from ..shared.climatology import get_lat_lon_comparison_descriptor, \ - get_remapper, get_mpas_climatology_file_names, \ - get_observation_climatology_file_names, \ - cache_climatologies, update_start_end_year, \ - remap_and_write_climatology -from ..shared.grid import MpasMeshDescriptor, LatLonGridDescriptor +from ..shared.climatology import MpasClimatology, ObservationClimatology +from ..shared.grid import LatLonGridDescriptor from ..shared.plot.plotting import plot_polar_comparison, \ setup_colormap -from ..shared.io.utility import build_config_full_path +from ..shared.io import build_config_full_path from ..shared.generalized_reader.generalized_reader \ import open_multifile_dataset @@ -27,7 +23,7 @@ from .sea_ice_analysis_task import SeaIceAnalysisTask -class ClimatologyMapSeaIce(SeaIceAnalysisTask): +class ClimatologyMapSeaIce(SeaIceAnalysisTask): # {{{ """ General comparison of 2-d model fields against data. Currently only supports sea ice concentration and sea ice thickness @@ -55,7 +51,8 @@ def __init__(self, config): # {{{ config=config, taskName='climatologyMapSeaIceConcThick', componentName='seaIce', - tags=['climatology', 'horizontalMap']) + tags=['climatology', 'horizontalMap'], + prerequisiteTasks=['cacheSeaIceTimeSeriesStatsTimes']) # }}} @@ -90,57 +87,12 @@ def run(self): # {{{ print "\nPlotting 2-d maps of sea-ice concentration and thickness " \ "climatologies..." - # get a list of timeSeriesStatsMonthly output files from the streams - # file, reading only those that are between the start and end dates - startDate = self.config.get('climatology', 'startDate') - endDate = self.config.get('climatology', 'endDate') - streamName = self.historyStreams.find_stream( - self.streamMap['timeSeriesStats']) - fileNames = self.historyStreams.readpath(streamName, - startDate=startDate, - endDate=endDate, - calendar=self.calendar) - print '\n Reading files:\n' \ - ' {} through\n {}'.format( - os.path.basename(fileNames[0]), - os.path.basename(fileNames[-1])) - # Load data - print ' Load sea-ice data...' - self.ds = open_multifile_dataset( - fileNames=fileNames, calendar=self.calendar, config=self.config, - simulationStartTime=self.simulationStartTime, - timeVariableName='Time', - variableList=['iceAreaCell', 'iceVolumeCell'], - variableMap=self.variableMap, startDate=startDate, - endDate=endDate) - - # Compute climatologies (first motnhly and then seasonally) - print ' Compute seasonal climatologies...' - - changed, startYear, endYear = update_start_end_year(self.ds, - self.config, - self.calendar) - - mpasDescriptor = MpasMeshDescriptor( - self.restartFileName, - meshName=self.config.get('input', 'mpasMeshName')) - - comparisonDescriptor = get_lat_lon_comparison_descriptor(self.config) - parallel = self.config.getint('execute', 'parallelTaskCount') > 1 if parallel: # avoid writing the same mapping file from multiple processes - mappingFilePrefix = 'map_{}'.format(self.taskName) + self.mappingFilePrefix = 'map_{}'.format(self.taskName) else: - mappingFilePrefix = 'map' - - self.mpasRemapper = get_remapper( - config=self.config, sourceDescriptor=mpasDescriptor, - comparisonDescriptor=comparisonDescriptor, - mappingFileSection='climatology', - mappingFileOption='mpasMappingFile', - mappingFilePrefix=mappingFilePrefix, - method=self.config.get('climatology', 'mpasInterpolationMethod')) + self.mappingFilePrefix = 'map' self._compute_and_plot_concentration() @@ -159,17 +111,8 @@ def _compute_and_plot_concentration(self): print ' Make ice concentration plots...' config = self.config - calendar = self.calendar - ds = self.ds mainRunName = config.get('runs', 'mainRunName') - startYear = config.getint('climatology', 'startYear') - endYear = config.getint('climatology', 'endYear') - overwriteMpasClimatology = config.getWithDefault( - 'climatology', 'overwriteMpasClimatology', False) - - overwriteObsClimatology = config.getWithDefault( - 'seaIceObservations', 'overwriteObsClimatology', False) subtitle = 'Ice concentration' @@ -178,94 +121,40 @@ def _compute_and_plot_concentration(self): 'DJF': ('SH', 'Winter'), 'JJA': ('SH', 'Summer')} - obsFileNames = {} - regriddedObsFileNames = {} - obsRemappers = {} - - comparisonDescriptor = self.mpasRemapper.destinationDescriptor - - buildObsClimatologies = overwriteObsClimatology for months in hemisphereSeasons: hemisphere, season = hemisphereSeasons[months] + # used in _open_sea_ice_part + self.fieldName = 'iceAreaCell' climFieldName = 'iceConcentration' - for obsName in ['NASATeam', 'Bootstrap']: - key = (months, obsName) - obsFileName = build_config_full_path( - config, 'seaIceObservations', - 'concentration{}{}_{}'.format(obsName, hemisphere, months)) - obsFieldName = '{}_{}_{}'.format(climFieldName, hemisphere, - obsName) - - obsDescriptor = LatLonGridDescriptor() - obsDescriptor.read(fileName=obsFileName, latVarName='t_lat', - lonVarName='t_lon') - obsRemapper = get_remapper( - config=config, sourceDescriptor=obsDescriptor, - comparisonDescriptor=comparisonDescriptor, - mappingFileSection='seaIceObservations', - mappingFileOption='seaIceClimatologyMappingFile', - mappingFilePrefix='map_obs_seaIce', - method=config.get('seaIceObservations', - 'interpolationMethod')) - obsRemappers[key] = obsRemapper - - if not os.path.isfile(obsFileName): - raise OSError('Obs file {} not found.'.format( - obsFileName)) - (climatologyFileName, regriddedFileName) = \ - get_observation_climatology_file_names( - config=config, fieldName=obsFieldName, - monthNames=months, componentName=self.componentName, - remapper=obsRemapper) - - obsFileNames[key] = obsFileName - regriddedObsFileNames[key] = regriddedFileName - - if not os.path.exists(regriddedFileName): - buildObsClimatologies = True - - for months in hemisphereSeasons: - hemisphere, season = hemisphereSeasons[months] - monthValues = constants.monthDictionary[months] - field = 'iceAreaCell' - climFieldName = 'iceConcentration' - - # interpolate the model results - mpasMeshName = self.mpasRemapper.sourceDescriptor.meshName - comparisonGridName = \ - self.mpasRemapper.destinationDescriptor.meshName - (climatologyFileName, climatologyPrefix, regriddedFileName) = \ - get_mpas_climatology_file_names( - config=config, - fieldName=climFieldName, - monthNames=months, - mpasMeshName=mpasMeshName, - comparisonGridName=comparisonGridName) - - if (overwriteMpasClimatology or - not os.path.exists(regriddedFileName)): - seasonalClimatology = cache_climatologies( - ds, monthValues, config, climatologyPrefix, calendar, - printProgress=True) - if seasonalClimatology is None: + mpasClimatology = MpasClimatology( + task=self, + fieldName=climFieldName, + monthNames=months, + streamName='timeSeriesStats', + meshFileName=self.restartFileName, + comparisonGrid='latlon', + mappingFileSection='climatology', + mappingFileOption='mpasMappingFile', + mappingFilePrefix=self.mappingFilePrefix, + method=config.get('climatology', 'mpasInterpolationMethod')) + + if mpasClimatology.remappedDataSet is None: + mpasClimatology.cache(openDataSetFunc=self._open_sea_ice_part, + printProgress=True) + if mpasClimatology.dataSet is None: # apparently, there was no data available to create the # climatology warnings.warn('no data to create sea ice concentration ' 'climatology for {}'.format(months)) continue - remappedClimatology = remap_and_write_climatology( - config, seasonalClimatology, climatologyFileName, - regriddedFileName, self.mpasRemapper) + mpasClimatology.remap_and_write() - else: - - remappedClimatology = xr.open_dataset(regriddedFileName) - - iceConcentration = remappedClimatology[field].values - lon = remappedClimatology['lon'].values - lat = remappedClimatology['lat'].values + iceConcentration = \ + mpasClimatology.remappedDataSet[self.fieldName].values + lon = mpasClimatology.remappedDataSet['lon'].values + lat = mpasClimatology.remappedDataSet['lat'].values lonTarg, latTarg = np.meshgrid(lon, lat) @@ -292,29 +181,47 @@ def _compute_and_plot_concentration(self): # ice concentrations from NASATeam (or Bootstrap) algorithm for obsName in ['NASATeam', 'Bootstrap']: - obsFieldName = 'AICE' - - key = (months, obsName) - regriddedFileName = regriddedObsFileNames[key] - - if buildObsClimatologies: - obsFileName = obsFileNames[key] + obsFileName = build_config_full_path( + config, 'seaIceObservations', + 'concentration{}{}_{}'.format(obsName, hemisphere, months)) + if not os.path.isfile(obsFileName): + raise OSError('Obs file {} not found.'.format( + obsFileName)) - seasonalClimatology = xr.open_dataset(obsFileName) + obsFieldName = '{}_{}_{}'.format(climFieldName, hemisphere, + obsName) + obsDescriptor = LatLonGridDescriptor.read(fileName=obsFileName, + latVarName='t_lat', + lonVarName='t_lon') + + obsClimatology = \ + ObservationClimatology( + task=self, + fieldName=obsFieldName, + monthNames=months, + obsGridDescriptor=obsDescriptor, + comparisonGrid='latlon', + mappingFileSection='seaIceObservations', + mappingFileOption='seaIceClimatologyMappingFile', + mappingFilePrefix='map_obs_seaIce', + method=config.get('seaIceObservations', + 'interpolationMethod')) - remappedClimatology = remap_and_write_climatology( - config, seasonalClimatology, climatologyFileName, - regriddedFileName, obsRemappers[key]) + if obsClimatology.remappedDataSet is None: + obsClimatology.dataSet = xr.open_dataset(obsFileName) + obsClimatology.remap_and_write() - obsIceConcentration = remappedClimatology[obsFieldName].values + obsIceConcentration = \ + obsClimatology.remappedDataSet['AICE'].values difference = iceConcentration - obsIceConcentration title = '{} ({}, years {:04d}-{:04d})'.format( - subtitle, months, startYear, endYear) + subtitle, months, mpasClimatology.startYear, + mpasClimatology.endYear) fileout = '{}/iceconc{}{}_{}_{}_years{:04d}-{:04d}.png'.format( self.plotsDirectory, obsName, hemisphere, mainRunName, - months, startYear, endYear) + months, mpasClimatology.startYear, mpasClimatology.endYear) plot_polar_comparison( config, lonTarg, @@ -349,113 +256,50 @@ def _compute_and_plot_thickness(self): print ' Make ice thickness plots...' config = self.config - calendar = self.calendar - ds = self.ds subtitle = 'Ice thickness' plotsDirectory = build_config_full_path(config, 'output', 'plotsSubdirectory') mainRunName = config.get('runs', 'mainRunName') - startYear = config.getint('climatology', 'startYear') - endYear = config.getint('climatology', 'endYear') - overwriteMpasClimatology = config.getWithDefault( - 'climatology', 'overwriteMpasClimatology', False) - - overwriteObsClimatology = config.getWithDefault( - 'seaIceObservations', 'overwriteObsClimatology', False) - - obsFileNames = {} - regriddedObsFileNames = {} - obsRemappers = {} - - comparisonDescriptor = self.mpasRemapper.destinationDescriptor - - # build a list of regridded observations files - buildObsClimatologies = overwriteObsClimatology - for months in ['FM', 'ON']: - climFieldName = 'iceThickness' - for hemisphere in ['NH', 'SH']: - key = (months, hemisphere) - obsFileName = build_config_full_path( - config, 'seaIceObservations', - 'thickness{}_{}'.format(hemisphere, months)) - if not os.path.isfile(obsFileName): - raise OSError('Obs file {} not found.'.format( - obsFileName)) - - obsFieldName = '{}_{}'.format(climFieldName, hemisphere) - obsDescriptor = LatLonGridDescriptor() - obsDescriptor.read(fileName=obsFileName, latVarName='t_lat', - lonVarName='t_lon') - obsRemapper = get_remapper( - config=config, sourceDescriptor=obsDescriptor, - comparisonDescriptor=comparisonDescriptor, - mappingFileSection='seaIceObservations', - mappingFileOption='seaIceClimatologyMappingFile', - mappingFilePrefix='map_obs_seaIce', - method=config.get('seaIceObservations', - 'interpolationMethod')) - obsRemappers[key] = obsRemapper - - (climatologyFileName, regriddedFileName) = \ - get_observation_climatology_file_names( - config=config, fieldName=obsFieldName, - monthNames=months, componentName=self.componentName, - remapper=obsRemapper) - - obsFileNames[key] = obsFileName - regriddedObsFileNames[key] = regriddedFileName - - if not os.path.exists(regriddedFileName): - buildObsClimatologies = True for months in ['FM', 'ON']: - monthValues = constants.monthDictionary[months] - field = 'iceVolumeCell' + self.fieldName = 'iceVolumeCell' climFieldName = 'iceThickness' - # interpolate the model results - mpasMeshName = self.mpasRemapper.sourceDescriptor.meshName - comparisonGridName = \ - self.mpasRemapper.destinationDescriptor.meshName - (climatologyFileName, climatologyPrefix, regriddedFileName) = \ - get_mpas_climatology_file_names( - config=config, - fieldName=climFieldName, - monthNames=months, - mpasMeshName=mpasMeshName, - comparisonGridName=comparisonGridName) - - if (overwriteMpasClimatology or - not os.path.exists(climatologyFileName)): - seasonalClimatology = cache_climatologies( - ds, monthValues, config, climatologyPrefix, calendar, - printProgress=True) - if seasonalClimatology is None: + mpasClimatology = MpasClimatology( + task=self, + fieldName=climFieldName, + monthNames=months, + streamName='timeSeriesStats', + meshFileName=self.restartFileName, + comparisonGrid='latlon', + mappingFileSection='climatology', + mappingFileOption='mpasMappingFile', + mappingFilePrefix=self.mappingFilePrefix, + method=config.get('climatology', 'mpasInterpolationMethod')) + + if mpasClimatology.remappedDataSet is None: + mpasClimatology.cache(openDataSetFunc=self._open_sea_ice_part, + printProgress=True) + if mpasClimatology.dataSet is None: # apparently, there was no data available to create the # climatology - warnings.warn('no data to create sea ice thickness ' + warnings.warn('no data to create sea ice concentration ' 'climatology for {}'.format(months)) continue - remappedClimatology = remap_and_write_climatology( - config, seasonalClimatology, climatologyFileName, - regriddedFileName, self.mpasRemapper) - - else: - - remappedClimatology = xr.open_dataset(regriddedFileName) + mpasClimatology.remap_and_write() - iceThickness = remappedClimatology[field].values + iceThickness = \ + mpasClimatology.remappedDataSet[self.fieldName].values iceThickness = ma.masked_values(iceThickness, 0) - lon = remappedClimatology['lon'].values - lat = remappedClimatology['lat'].values + lon = mpasClimatology.remappedDataSet['lon'].values + lat = mpasClimatology.remappedDataSet['lat'].values lonTarg, latTarg = np.meshgrid(lon, lat) for hemisphere in ['NH', 'SH']: - obsFieldName = 'HI' (colormapResult, colorbarLevelsResult) = setup_colormap( config, @@ -473,20 +317,38 @@ def _compute_and_plot_thickness(self): 'climatologyMapSeaIceConcThick', 'minimumLatitude{}'.format(hemisphere)) - # now the observations - key = (months, hemisphere) - regriddedFileName = regriddedObsFileNames[key] + obsFileName = build_config_full_path( + config, 'seaIceObservations', + 'thickness{}_{}'.format(hemisphere, months)) + if not os.path.isfile(obsFileName): + raise OSError('Obs file {} not found.'.format( + obsFileName)) - if buildObsClimatologies: - obsFileName = obsFileNames[key] + obsFieldName = '{}_{}'.format(climFieldName, hemisphere) - seasonalClimatology = xr.open_dataset(obsFileName) + obsDescriptor = LatLonGridDescriptor.read(fileName=obsFileName, + latVarName='t_lat', + lonVarName='t_lon') + + obsClimatology = \ + ObservationClimatology( + task=self, + fieldName=obsFieldName, + monthNames=months, + obsGridDescriptor=obsDescriptor, + comparisonGrid='latlon', + mappingFileSection='seaIceObservations', + mappingFileOption='seaIceClimatologyMappingFile', + mappingFilePrefix='map_obs_seaIce', + method=config.get('seaIceObservations', + 'interpolationMethod')) - remappedClimatology = remap_and_write_climatology( - config, seasonalClimatology, climatologyFileName, - regriddedFileName, obsRemappers[key]) + if obsClimatology.remappedDataSet is None: + obsClimatology.dataSet = xr.open_dataset(obsFileName) + obsClimatology.remap_and_write() - obsIceThickness = remappedClimatology[obsFieldName].values + obsIceThickness = \ + obsClimatology.remappedDataSet['HI'].values # Mask thickness fields obsIceThickness = ma.masked_values(obsIceThickness, 0) @@ -499,12 +361,12 @@ def _compute_and_plot_thickness(self): difference = iceThickness - obsIceThickness - title = '{} ({}, years {:04d}-{:04d})'.format(subtitle, months, - startYear, - endYear) + title = '{} ({}, years {:04d}-{:04d})'.format( + subtitle, months, mpasClimatology.startYear, + mpasClimatology.endYear) fileout = '{}/icethick{}_{}_{}_years{:04d}-{:04d}.png'.format( - plotsDirectory, hemisphere, mainRunName, months, startYear, - endYear) + plotsDirectory, hemisphere, mainRunName, months, + mpasClimatology.startYear, mpasClimatology.endYear) plot_polar_comparison( config, lonTarg, @@ -528,4 +390,39 @@ def _compute_and_plot_thickness(self): # }}} + def _open_sea_ice_part(self, inputFileNames, startDate, endDate): # {{{ + """ + Open part of a sea-ice data set between the given start and end date, + used to cache a climatology of the data set. + + Parameters + ---------- + inputFileNames : list of str + File names in the multifile data set to open + + startDate, endDate : float + start and end date to which to crop the Time dimension (given in + days since 0001-01-01) + + Authors + ------- + Xylar Asay-Davis + """ + + variableList = [self.fieldName] + + ds = open_multifile_dataset( + fileNames=inputFileNames, + calendar=self.calendar, + config=self.config, + simulationStartTime=self.simulationStartTime, + timeVariableName='Time', + variableList=variableList, + variableMap=self.variableMap, + startDate=startDate, + endDate=endDate) + + return ds # }}} + + # }}} # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python \ No newline at end of file diff --git a/mpas_analysis/sea_ice/sea_ice_analysis_task.py b/mpas_analysis/sea_ice/sea_ice_analysis_task.py index 53897bb6b..7ceeaeb48 100644 --- a/mpas_analysis/sea_ice/sea_ice_analysis_task.py +++ b/mpas_analysis/sea_ice/sea_ice_analysis_task.py @@ -1,7 +1,6 @@ from ..shared.analysis_task import AnalysisTask -from ..shared.io import StreamsFile -from ..shared.io.utility import build_config_full_path +from ..shared.io import StreamsFile, build_config_full_path from ..shared.timekeeping.utility import get_simulation_start_time diff --git a/mpas_analysis/sea_ice/time_series.py b/mpas_analysis/sea_ice/time_series.py index 41b5c2379..b50002e5c 100644 --- a/mpas_analysis/sea_ice/time_series.py +++ b/mpas_analysis/sea_ice/time_series.py @@ -6,7 +6,7 @@ from ..shared.plot.plotting import timeseries_analysis_plot, \ timeseries_analysis_plot_polar -from ..shared.io.utility import build_config_full_path, check_path_exists, \ +from ..shared.io import build_config_full_path, check_path_exists, \ make_directories from ..shared.timekeeping.utility import date_to_days, days_to_datetime, \ diff --git a/mpas_analysis/shared/analysis_task.py b/mpas_analysis/shared/analysis_task.py index 3596a34ee..71fe3d014 100644 --- a/mpas_analysis/shared/analysis_task.py +++ b/mpas_analysis/shared/analysis_task.py @@ -4,13 +4,22 @@ Authors ------- Xylar Asay-Davis - ''' import warnings +import pickle +import os +from collections import OrderedDict +import numpy + +from .constants import constants +from .generalized_reader import open_multifile_dataset + +from .timekeeping.utility import days_to_datetime, string_to_days_since_date, \ + add_years_months_days_in_month -from .io import NameList, StreamsFile -from .io.utility import build_config_full_path, make_directories +from .io import NameList, StreamsFile, build_config_full_path, \ + make_directories from .variable_namelist_stream_maps.ocean_maps import oceanNamelistMap, \ oceanStreamMap, oceanVariableMap @@ -23,12 +32,74 @@ class AnalysisTask(object): # {{{ ''' The base class for analysis tasks. + Attributes + ---------- + config : instance of MpasAnalysisConfigParser + Contains configuration options + + taskName : str + The name of the task, typically the same as the class name except + starting with lowercase (e.g. 'myTask' for class 'MyTask') + + componentName : {'ocean', 'seaIce'} + The name of the component (same as the folder where the task + resides) + + tags : list of str + Tags used to describe the task (e.g. 'timeSeries', 'climatology', + horizontalMap', 'index', 'transect'). These are used to determine + which tasks are generated (e.g. 'all_transect' or 'no_climatology' + in the 'generate' flags) + + prerequisiteTasks : list of str, optional + Names of tasks that must complete before this task can run. + Typically, this will include one or more tasks of the form + ``cacheTimes``, e.g. + ``cacheOceanTimeSeriesStatsTimes`` + + runDirectory : str + the base input directory for namelists, streams files and restart files + + historyDirectory : str + the base input directory for history files + + plotsDirectory : str + the directory for writing plots (which is also created if it doesn't + exist) + + namelist : ``NameList`` object + the namelist reader + + runStreams : ``StreamsFile`` object + the streams file reader for streams in the run directory (e.g. restart + files) + + historyStreams : ``StreamsFile`` object + the streams file reader for streams in the history directory (most + streams other than restart files) + + calendar : {'gregorian', 'gregorian_noleap'} + the name of the calendar + + namelistMap : dict + A map between names of namelist options used by MPAS-Analysis and + those in various MPAS versions + + streamMap : dict + a map between names of streams used by MPAS-Analysis and those in + various MPAS versions + + variableMap : dict + a map between names of variables within streams used by MPAS-Analysis + and those in various MPAS versions + Authors ------- Xylar Asay-Davis ''' - def __init__(self, config, taskName, componentName, tags=[]): # {{{ + def __init__(self, config, taskName, componentName, tags=[], + prerequisiteTasks=None): # {{{ ''' Construct the analysis task. @@ -55,6 +126,12 @@ def __init__(self, config, taskName, componentName, tags=[]): # {{{ which tasks are generated (e.g. 'all_transect' or 'no_climatology' in the 'generate' flags) + prerequisiteTasks : list of str, optional + Names of tasks that must complete before this task can run. + Typically, this will include one or more tasks of the form + ``cacheTimes``, e.g. + ``cacheOceanTimeSeriesStatsTimes`` + Authors ------- Xylar Asay-Davis @@ -62,33 +139,14 @@ def __init__(self, config, taskName, componentName, tags=[]): # {{{ self.config = config self.taskName = taskName self.componentName = componentName - self.tags = tags # }}} + self.tags = tags + self.prerequisiteTasks = prerequisiteTasks # }}} def setup_and_check(self): # {{{ ''' Perform steps to set up the analysis (e.g. reading namelists and streams files). - After this call, the following member variables are set: - self.runDirectory : the base input directory for namelists, streams - files and restart files - self.historyDirectory : the base input directory for history files - self.plotsDirectory : the directory for writing plots (which is - also created if it doesn't exist) - self.namelist : the namelist reader - self.runStreams : the streams file reader for streams in the run - directory (e.g. restart files) - self.historyStreams : the streams file reader for streams in the - history directory (most streams other than restart files) - self.calendar : the name of the calendar ('gregorian' or - 'gregoraian_noleap') - self.namelistMap : a map between names of namelist options used by - MPAS-Analysis and those in various MPAS versions - self.streamMap : a map between names of streams used by - MPAS-Analysis and those in various MPAS versions - self.variableMap : a map between names of variables within streams - used by MPAS-Analysis and those in various MPAS versions - Individual tasks (children classes of this base class) should first call this method to perform basic setup, then, check whether the configuration is correct for a given analysis and perform additional, @@ -96,6 +154,51 @@ def setup_and_check(self): # {{{ necessary observations and other data files are found, then, determine the list of files to be read when the analysis is run. + If the task includes ``climatology``, ``timeSeries`` or ``index`` tags, + ``startDate`` and ``endDate`` config options are computed from + ``startYear`` and ``endYear``config options. + + After this call, the following attributes are set. + + Attributes + ---------- + runDirectory : str + the base input directory for namelists, streams files and restart + files + + historyDirectory : str + the base input directory for history files + + plotsDirectory : str + the directory for writing plots (which is also created if it + doesn't exist) + + namelist : ``NameList`` object + the namelist reader + + runStreams : ``StreamsFile`` object + the streams file reader for streams in the run directory (e.g. + restart files) + + historyStreams : ``StreamsFile`` object + the streams file reader for streams in the history directory (most + streams other than restart files) + + calendar : {'gregorian', 'gregorian_noleap'} + the name of the calendar + + namelistMap : dict + A map between names of namelist options used by MPAS-Analysis and + those in various MPAS versions + + streamMap : dict + a map between names of streams used by MPAS-Analysis and those in + various MPAS versions + + variableMap : dict + a map between names of variables within streams used by + MPAS-Analysis and those in various MPAS versions + Authors ------- Xylar Asay-Davis @@ -311,6 +414,263 @@ def set_start_end_date(self, section): # {{{ self.config.getint(section, 'endYear')) self.config.set(section, 'endDate', endDate) # }}} + def update_start_end_date(self, section, streamName): # {{{ + ''' + Update the start and end dates (and years) based on the times found + in the given stream. Cache the times if they are not already cached. + + Parameters + ---------- + section : str + The name of a section in the config file containing ``startYear`` + and ``endYear`` options. ``section`` is typically one of + ``climatology``, ``timeSeries`` or ``index`` + + streamName : str + The name of a stream from which to read (and cache) the times + + Returns + ------- + changed : bool + Whether the start and end dates were updated. + + Authors + ------- + Xylar Asay-Davis + ''' + + startDate = self.config.get(section, 'startDate') + endDate = self.config.get(section, 'endDate') + startDate = string_to_days_since_date(dateString=startDate, + calendar=self.calendar) + endDate = string_to_days_since_date(dateString=endDate, + calendar=self.calendar) + + inFileNames = self.get_input_file_names( + streamName, startAndEndDateSection=section) + + fullTimeCache = self.cache_multifile_dataset_times( + inFileNames, streamName, timeVariableName='Time') + + # find only those cached times between starDate and endDate + times = [] + for fileName in fullTimeCache: + localTimes = fullTimeCache[fileName]['times'] + mask = numpy.logical_and(localTimes >= startDate, + localTimes < endDate) + if numpy.count_nonzero(mask) == 0: + continue + + times.extend(list(localTimes[mask])) + + requestedStartYear = self.config.getint('climatology', 'startYear') + requestedEndYear = self.config.getint('climatology', 'endYear') + + startYear = days_to_datetime(numpy.amin(times), + calendar=self.calendar).year + endYear = days_to_datetime(numpy.amax(times), + calendar=self.calendar).year + changed = False + if startYear != requestedStartYear or endYear != requestedEndYear: + message = "{} start and/or end year different from " \ + "requested\n" \ + "requested: {:04d}-{:04d}\n" \ + "actual: {:04d}-{:04d}\n".format(section, + requestedStartYear, + requestedEndYear, + startYear, + endYear) + warnings.warn(message) + self.config.set(section, 'startYear', str(startYear)) + self.config.set(section, 'endYear', str(endYear)) + + startDate = '{:04d}-01-01_00:00:00'.format(startYear) + self.config.set(section, 'startDate', startDate) + endDate = '{:04d}-12-31_23:59:59'.format(endYear) + self.config.set(section, 'endDate', endDate) # }}} + + changed = True + + return changed # }}} + + def get_input_file_names(self, streamName, + startDate=None, endDate=None, + startAndEndDateSection=None): # {{{ + ''' + Get a list of input files corresponding to the given stream and + optionally bounded by the start and end dates found in the given + section of the config file. + + Parameters + ---------- + streamName : str + The name of a stream to check. If ``self.streamMap`` is defined, + the streamName will be mapped to the corresponding name in the + streams file + + startDate, endDate : float, optional + start and end date to use in determining which files to include in + the list + + startAndEndDateSection : str, optional + If ``startDate`` and ``endDate`` arguments are not supplied, the + name of a section in the config file containing ``startDate`` and + ``endDate`` options to use instead. ``startAndEndDateSection`` is + typically one of ``climatology``, ``timeSeries`` or ``index``. + + Raises + ------ + RuntimeError + If no files are found in the desired date range. + + Authors + ------- + Xylar Asay-Davis + ''' + + if startDate is None and endDate is None and \ + startAndEndDateSection is not None: + startDate = self.config.get(startAndEndDateSection, 'startDate') + endDate = self.config.get(startAndEndDateSection, 'endDate') + + if self.streamMap is not None: + streamName = \ + self.historyStreams.find_stream(self.streamMap[streamName]) + inputFileNames = self.historyStreams.readpath(streamName, + startDate=startDate, + endDate=endDate, + calendar=self.calendar) + + if len(inputFileNames) == 0: + raise RuntimeError('No input files could be found in stream {} ' + 'between {} and {}'.format(streamName, + startDate, endDate)) + return inputFileNames # }}} + + def cache_multifile_dataset_times(self, inFileNames, streamName, + timeVariableName='Time'): # {{{ + """ + Creates a cache file of the times in each file of a multifile data set. + This is useful when caching climatologies and time series as a + simulation evolves, since files that have already been processed will + not need to be opened to find out which times they contain. + + Parameters + ---------- + inFileNames : list of str + A list of file paths to read + + streamName : str + The name of a stream, used to build the name of the cache file + + timeVariableName : string, optional + The name of the time variable (typically 'Time' if using a + variableMap or 'xtime' if not using a variableMap) + + Author + ------ + Xylar Asay-Davis + """ + + timeCacheDirectory = build_config_full_path( + self.config, 'output', 'timeCacheSubdirectory') + + make_directories(timeCacheDirectory) + + cacheFileName = '{}/{}_{}_times.pickle'.format(timeCacheDirectory, + self.componentName, + streamName) + if os.path.exists(cacheFileName): + with open(cacheFileName, 'rb') as handle: + inTimeCache = pickle.load(handle) + else: + inTimeCache = OrderedDict() + + # add files already in the time cache to the list of files to check + # (and potentially update) + fileNames = list(set(inFileNames + inTimeCache.keys())) + fileNames.sort() + + fileNames = [os.path.abspath(fileName) for fileName in fileNames] + + if hasattr(self, 'simulationStartTime'): + simulationStartTime = self.simulationStartTime + else: + simulationStartTime = None + + filesToRead = [] + for fileName in fileNames: + read = True + if fileName in inTimeCache.keys(): + dateModified = inTimeCache[fileName]['dateModified'] + if dateModified == os.path.getmtime(fileName): + # the file is already in the cache and hasn't been + # changed so we don't need to update the times in the + # cache + read = False + if read: + filesToRead.append(fileName) + + if len(filesToRead) == 0: + outTimeCache = inTimeCache + else: + print '\n Caching times from files:\n' \ + ' {} through\n {}'.format( + os.path.basename(filesToRead[0]), + os.path.basename(filesToRead[-1])) + + outTimeCache = OrderedDict() + for fileName in fileNames: + read = True + if fileName in filesToRead: + # open the files one at a time so we know which time is in + # which file + ds = open_multifile_dataset( + fileNames=[fileName], + calendar=self.calendar, + config=self.config, + simulationStartTime=simulationStartTime, + timeVariableName=timeVariableName, + variableList=['Time'], + variableMap=self.variableMap) + + ds = add_years_months_days_in_month(ds, self.calendar) + + dateModified = os.path.getmtime(fileName) + times = ds.Time.values + datetimes = days_to_datetime(times, calendar=self.calendar) + years = numpy.array([date.year for date in datetimes]) + months = numpy.array([date.month for date in datetimes]) + if 'startTime' in ds.coords and 'endTime' in ds.coords: + daysInMonth = ds.endTime.values - ds.startTime.values + else: + if self.calendar == 'gregorian': + message = 'The MPAS run used the Gregorian ' \ + 'calendar but does not appear to ' \ + 'have\n' \ + 'supplied start and end times. ' \ + 'Climatologies will be computed with\n' \ + 'month durations ignoring leap years.' + warnings.warn(message) + + daysInMonth = numpy.array( + [constants.daysInMonth[month-1] for month + in ds.month.values], float) + del ds + outTimeCache[fileName] = {'times': times, + 'years': years, + 'months': months, + 'daysInMonth': daysInMonth, + 'dateModified': dateModified} + else: + outTimeCache[fileName] = inTimeCache[fileName] + + with open(cacheFileName, 'wb') as handle: + pickle.dump(outTimeCache, handle, + protocol=pickle.HIGHEST_PROTOCOL) + + return outTimeCache # }}} + # }}} diff --git a/mpas_analysis/shared/cache_dataset_times_task.py b/mpas_analysis/shared/cache_dataset_times_task.py new file mode 100644 index 000000000..9c7ffc1d5 --- /dev/null +++ b/mpas_analysis/shared/cache_dataset_times_task.py @@ -0,0 +1,153 @@ +from .analysis_task import AnalysisTask + +from ..shared.io import StreamsFile, build_config_full_path +from ..shared.timekeeping.utility import get_simulation_start_time + + +class CacheDatasetTimesTask(AnalysisTask): # {{{ + ''' + A task for caching the times in a multifile data sets of an MPAS analysis + member for later use. Since analysis member may be used by multiple tasks + (indeed the ``timeSeriesStats`` is currently used by *all* tasks), it is + important that this time information gets processed once before all other + tasks get run in parallel. + + Authors + ------- + Xylar Asay-Davis + ''' + + def __init__(self, config, componentName, streamName, + startAndEndDateSections, namelistOption=None): # {{{ + ''' + Construct an analysis task for caching the times in the multifile + data set in the given component and stream. The name of the task + includes the component and stream name. For example, if + ``component='ocean'`` and ``streamName='timeSeriesStats``, then the + task name is ``cacheOceanTimeSeriesStatsTimes``. + + Parameters + ---------- + config : instance of MpasAnalysisConfigParser + Contains configuration options + + componentName : {'ocean', 'seaIce'} + The name of the component (same as the folder where the task + resides) + + streamName : str + The name of the stream from which the climatology data set will + be read, used to cache times and update their bounds in the + configuration parser. + + startAndEndDateSections : list, {'climatology', 'timeSeries', 'index'} + The name of sections in the config file containing ``startDate`` + and ``endDate`` options as a list. + + namelistOption : str, optional + The name of a namelist option (e.g. + ``config_am_timeseriesstatsmonthly_enable``) that should be set to + true of the required analysis member has been enabled. If this + option is ``None`` (the default), no check is performed + + Authors + ------- + Xylar Asay-Davis + ''' + + upperComponent = componentName[0].upper() + componentName[1:] + upperStream = streamName[0].upper() + streamName[1:] + + taskName = 'cache{}{}Times'.format( + upperComponent, upperStream) + + # first, call the constructor from the base class (AnalysisTask). + super(CacheDatasetTimesTask, self).__init__( + config=config, + taskName=taskName, + componentName=componentName, + tags=startAndEndDateSections) + + self.streamName = streamName + self.namelistOption = namelistOption + self.startAndEndDateSections = startAndEndDateSections + + # }}} + + def setup_and_check(self): # {{{ + ''' + Perform steps to set up the analysis and check for errors in the setup. + + Raises + ------ + ValueError: if startAndEndDateSections is not a list containing + {'climatology', 'timeSeries', 'index'} + + 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, self.namelistMap, self.streamMap, self.variableMap + super(CacheDatasetTimesTask, self).setup_and_check() + + if self.namelistOption is not None: + self.check_analysis_enabled( + analysisOptionName=self.namelistOption, + raiseException=True) + + for section in self.startAndEndDateSections: + if not self.config.has_section(section): + raise ValueError('Config file does not have a section ' + '{}.'.format(section)) + for option in ['startDate', 'endDate']: + if not self.config.has_option(section, option): + raise ValueError('Config section {} does not have ' + 'expected option {}.'.format(section, + option)) + + # }}} + + def run(self): # {{{ + ''' + The main method of the task that performs the analysis task. + + Authors + ------- + Xylar Asay-Davis + ''' + print "" + print "Caching times from the {} component for the {} stream " \ + "...".format(self.componentName, self.streamName) + + try: + self.simulationStartTime = get_simulation_start_time( + self.runStreams) + except IOError as e: + if self.componentName == 'ocean': + raise e + else: + # try the ocean stream instead + runDirectory = build_config_full_path(self.config, 'input', + 'runSubdirectory') + oceanStreamsFileName = build_config_full_path( + self.config, 'input', 'oceanStreamsFileName') + oceanStreams = StreamsFile(oceanStreamsFileName, + streamsdir=runDirectory) + self.simulationStartTime = \ + get_simulation_start_time(oceanStreams) + + for sectionName in self.startAndEndDateSections: + # internally, this will also cache the times for the data set + self.update_start_end_date(section=sectionName, + streamName=self.streamName) + + # }}} + +# }}} + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/shared/climatology/__init__.py b/mpas_analysis/shared/climatology/__init__.py index 2442345da..ddf70a536 100644 --- a/mpas_analysis/shared/climatology/__init__.py +++ b/mpas_analysis/shared/climatology/__init__.py @@ -1,5 +1,2 @@ -from .climatology import get_lat_lon_comparison_descriptor, get_remapper, \ - get_mpas_climatology_file_names, get_observation_climatology_file_names, \ - compute_monthly_climatology, compute_climatology, cache_climatologies, \ - update_start_end_year, add_years_months_days_in_month, \ - remap_and_write_climatology +from .climatology import Climatology, \ + MpasClimatology, ObservationClimatology \ No newline at end of file diff --git a/mpas_analysis/shared/climatology/climatology.py b/mpas_analysis/shared/climatology/climatology.py index c15364722..0bd3346a6 100644 --- a/mpas_analysis/shared/climatology/climatology.py +++ b/mpas_analysis/shared/climatology/climatology.py @@ -4,25 +4,26 @@ Authors ------- Xylar Asay-Davis - -Last Modified -------------- -04/13/2017 """ import xarray as xr import os import numpy import warnings +from collections import OrderedDict +import numbers from ..constants import constants -from ..timekeeping.utility import days_to_datetime +from ..timekeeping.utility import string_to_days_since_date, \ + add_years_months_days_in_month -from ..io.utility import build_config_full_path, make_directories, fingerprint_generator +from ..io import build_config_full_path, make_directories +from ..io.utility import fingerprint_generator from ..interpolation import Remapper -from ..grid import LatLonGridDescriptor, ProjectionGridDescriptor +from ..grid import MpasMeshDescriptor, LatLonGridDescriptor, \ + ProjectionGridDescriptor def get_lat_lon_comparison_descriptor(config): # {{{ @@ -43,10 +44,6 @@ def get_lat_lon_comparison_descriptor(config): # {{{ Authors ------- Xylar Asay-Davis - - Last Modified - ------------- - 04/13/2017 """ climSection = 'climatology' @@ -62,923 +59,1158 @@ def get_lat_lon_comparison_descriptor(config): # {{{ lat = numpy.linspace(constants.latmin, constants.latmax, nLat) lon = numpy.linspace(constants.lonmin, constants.lonmax, nLon) - descriptor = LatLonGridDescriptor() - descriptor.create(lat, lon, units='degrees') + descriptor = LatLonGridDescriptor.create(lat, lon, units='degrees') return descriptor # }}} -def get_remapper(config, sourceDescriptor, comparisonDescriptor, - mappingFileSection, mappingFileOption, mappingFilePrefix, - method): # {{{ +class Climatology(object): # {{{ """ - Given config options and descriptions of the source and comparison grids, - returns a ``Remapper`` object that can be used to remap from source files - or data sets to corresponding data sets on the comparison grid. - - If necessary, creates the mapping file containing weights and indices - needed to perform remapping. + A class for computing climatologies from a monthly mean data set. - Parameters + Attributes ---------- - config : instance of ``MpasAnalysisConfigParser`` - Contains configuration options - - sourceDescriptor : ``MeshDescriptor`` subclass object - A description of the source mesh or grid - - comparisonDescriptor : ``MeshDescriptor`` subclass object - A description of the comparison grid - - mappingFileSection, mappingFileOption : str - Section and option in ``config`` where the name of the mapping file - may be given, or where it will be stored if a new mapping file is - created - - mappingFilePrefix : str - A prefix to be prepended to the mapping file name - - method : {'bilinear', 'neareststod', 'conserve'} - The method of interpolation used. - - Returns - ------- - remapper : ``Remapper`` object - A remapper that can be used to remap files or data sets from the source - grid or mesh to the comparison grid. - - Authors - ------- - Xylar Asay-Davis - - Last Modified - ------------- - 04/13/2017 - """ - - if config.has_option(mappingFileSection, mappingFileOption): - # a mapping file was supplied, so we'll use that name - mappingFileName = config.get(mappingFileSection, mappingFileOption) - else: - if _matches_comparison(sourceDescriptor, comparisonDescriptor): - # no need to remap - mappingFileName = None - - else: - # we need to build the path to the mapping file and an appropriate - # file name - mappingSubdirectory = build_config_full_path(config, 'output', - 'mappingSubdirectory') - - make_directories(mappingSubdirectory) - - mappingFileName = '{}/{}_{}_to_{}_{}.nc'.format( - mappingSubdirectory, mappingFilePrefix, - sourceDescriptor.meshName, comparisonDescriptor.meshName, - method) - - config.set(mappingFileSection, mappingFileOption, - mappingFileName) - - remapper = Remapper(sourceDescriptor, comparisonDescriptor, - mappingFileName) - - remapper.build_mapping_file(method=method) - - return remapper # }}} - - -def get_mpas_climatology_file_names(config, fieldName, monthNames, - mpasMeshName, - comparisonGridName=None): # {{{ - """ - Given config options, the name of a field and a string identifying the - months in a seasonal climatology, returns the full path for MPAS - climatology files before and after regridding. + task : ``AnalysisTask`` object + An analysis task for which the climatology is needed. ``task`` is + used to get config options, define the component name, map streams, + namelists and variables, etc. - Parameters - ---------- config : instance of MpasAnalysisConfigParser Contains configuration options - fieldName : str - Name of the field being mapped, used as a prefix for the climatology - file name. + calendar : {'gregorian', 'gregorian_noleap'} + the name of the calendar monthNames : str - A string identifying the months in a seasonal climatology (e.g. 'JFM') - - mpasMeshName : str - The name of the MPAS mesh - - comparisonGridName : str, optional - The name of the comparison grid (if any) - - Returns - ------- - climatologyFileName : str - The absolute path to a file where the climatology should be stored - before regridding. - - climatologyPrefix : str - The prfix including absolute path for climatology cache files before - regridding. - - regriddedFileName : str - The absolute path to a file where the climatology should be stored - after regridding if ``comparisonGridName`` is supplied - - Authors - ------- - Xylar Asay-Davis - - Last Modified - ------------- - 05/05/2017 - """ - - climSection = 'climatology' - startYear = config.getint(climSection, 'startYear') - endYear = config.getint(climSection, 'endYear') - - climatologyDirectory = build_config_full_path( - config, 'output', 'mpasClimatologySubdirectory') - - make_directories(climatologyDirectory) - - climatologyPrefix = '{}/{}_{}_{}'.format(climatologyDirectory, fieldName, - mpasMeshName, monthNames) - - yearString, fileSuffix = _get_year_string(startYear, endYear) - climatologyFileName = '{}_{}.nc'.format(climatologyPrefix, fileSuffix) - - if comparisonGridName is None: - return (climatologyFileName, climatologyPrefix) - else: - regriddedDirectory = build_config_full_path( - config, 'output', 'mpasRegriddedClimSubdirectory') - - make_directories(regriddedDirectory) - - regriddedFileName = '{}/{}_{}_to_{}_{}_{}.nc'.format( - regriddedDirectory, fieldName, mpasMeshName, - comparisonGridName, monthNames, fileSuffix) - - return (climatologyFileName, climatologyPrefix, - regriddedFileName) - - # }}} - + The months that make up the climatology, used in constructing the + names of cache files. If provided, ``monthNames`` should be one of + the keys of ``monthDictionary`` in the ``constants`` module. -def get_observation_climatology_file_names(config, fieldName, monthNames, - componentName, remapper): # {{{ - """ - Given config options, the name of a field and a string identifying the - months in a seasonal climatology, returns the full path for observation - climatology files before and after regridding. + monthValues : list of int + A list of integer months that make up the season in ``monthNames``, + taken from ``monthDictionary`` in the ``constants`` module. The + entries are sorted in ascending order. - Parameters - ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + sourceDescriptor : ``MeshDescriptor`` object + A descriptor of the source grid or mesh used for remapping. This + attribute should be set by one of the subclasses of ``Climatology`` if + remapping will be performed. - fieldName : str - Name of the field being mapped, used as a prefix for the climatology - file name. - - monthNames : str - A string identifying the months in a seasonal climatology (e.g. 'JFM') + comparisonDescriptor : ``MeshDescriptor`` object + If the name of a comparison grid name was supplied, this is an object + describing that comparison grid (e.g. for remapping). If no comparison + grid name was supplied, this is ``None``. remapper : ``Remapper`` object - A remapper that used to remap files or data sets from the - observation grid to a comparison grid + A remapper between the source and comparison grids. Available only + after calling the ``create_remapper`` method. - Returns - ------- - climatologyFileName : str - The absolute path to a file where the climatology should be stored - before regridding. + dataSet : ``xarray.Dataset`` or ``xarray.DataArray`` object + A data set containing the climatology (before remapping, if + applicable). Available only after calling the ``compute`` or + ``compute_monthly`` method. - regriddedFileName : str - The absolute path to a file where the climatology should be stored - after regridding. + climatologyFileName : str + The name of the file where ``dataSet`` (the climatology data set before + remapping) should be stored. This attribute should be set by one of + the subclasses of ``Climatology`` if remapping will be performed. + + remappedDataSet : ``xarray.Dataset`` or ``xarray.DataArray`` object + A data set containing the remapped climatology. Available only after + calling the ``remap_and_write`` method. + + remappedFileName : str + The name of the file where ``remappedDataSet`` (the climatology data + set after remapping) should be stored. This attribute should be set by + one of the subclasses of ``Climatology`` if remapping will be + performed. + + Examples + -------- + Here is an example of how to use a ``Climatology`` object within a + task. In the example, ``ds`` is an xarray data set with a ``Time`` + dimension. The result of the example code is that the xarray data set + ``climatology.dataSet`` is available for computation and plotting. + + >>> from ..shared.climatology import Climatology + >>> climatology = Climatology( + task=self, + monthNames='ANN') + >>> climatology.compute(ds=ds) + >>> print climatology.dataSet + + Computing a monthly climatology would look like this: + >>> from ..shared.climatology import Climatology + >>> climatology = Climatology(task=self) + >>> climatology.compute_monthly(ds=ds) + >>> print climatology.dataSet Authors ------- Xylar Asay-Davis - - Last Modified - ------------- - 03/03/2017 """ - obsSection = '{}Observations'.format(componentName) - - climatologyDirectory = build_config_full_path( - config=config, section='output', - relativePathOption='climatologySubdirectory', - relativePathSection=obsSection) - - regriddedDirectory = build_config_full_path( - config=config, section='output', - relativePathOption='regriddedClimSubdirectory', - relativePathSection=obsSection) - - obsGridName = remapper.sourceDescriptor.meshName - comparisonGridName = remapper.destinationDescriptor.meshName + def __init__(self, task, monthNames=None, comparisonGrid=None): + """ + Create a new climatology. + + Parameters + ---------- + task : ``AnalysisTask`` object + An analysis task for which the climatology is needed. ``task`` is + used to get config options, define the component name, map streams, + namelists and variables, etc. + + monthNames : str, optional + The months that make up the climatology, used in constructing the + names of cache files. If provided, ``monthNames`` should be one of + the keys of ``monthDictionary`` in the ``constants`` module. + + comparisonGrid : {'latlon'}, optional + The name of the comparison grid to use for remapping (if any). + If ``comparisonGrid=None`` (the default), no remapping will be + performed. + + Raises + ------ + ValueError + If comarisonGrid does not describe a known comparions grid + + Authors + ------- + Xylar Asay-Davis + """ + self.task = task + self.config = task.config + self.calendar = task.calendar + self.monthNames = monthNames + if monthNames is None: + self.monthValues = None + else: + monthValues = constants.monthDictionary[monthNames] + if isinstance(monthValues, numbers.Integral): + self.monthValues = [monthValues] + else: + self.monthValues = sorted(monthValues) + if comparisonGrid == 'latlon': + self.comparisonDescriptor = \ + get_lat_lon_comparison_descriptor(task.config) + elif comparisonGrid is None: + self.comparisonDescriptor = None + else: + raise ValueError('Unknown comaprison grid type {}'.format( + comparisonGrid)) + + def create_remapper(self, mappingFileSection, mappingFileOption, + mappingFilePrefix, method): # {{{ + """ + Creates an attribute ``remapper``, a ``Remapper`` object, that can be + used to remap from source files or data sets to corresponding data sets + on the comparison grid. + + This call requires that the ``sourceDescriptor`` attribute be assigned, + which happens automatically for some of the subclasses of + ``Climatology`` (e.g. ``MpasClimatology`` and + ``ObservationClimatology``). + + If necessary, creates the mapping file containing weights and indices + needed to perform remapping. + + Parameters + ---------- + mappingFileSection, mappingFileOption : str + Section and option in ``config`` where the name of the mapping file + may be given, or where it will be stored if a new mapping file is + created + + mappingFilePrefix : str + A prefix to be prepended to the mapping file name + + method : {'bilinear', 'neareststod', 'conserve'} + The method of interpolation used. + + Returns + ------- + remapper : ``Remapper`` object + A remapper between the source and comparison grids, also stored + as a ``remapper`` attribute. + + Authors + ------- + Xylar Asay-Davis + """ + + config = self.config + + if config.has_option(mappingFileSection, mappingFileOption): + # a mapping file was supplied, so we'll use that name + mappingFileName = config.get(mappingFileSection, mappingFileOption) + else: + if self._matches_comparison(self.sourceDescriptor, + self.comparisonDescriptor): + # no need to remap + mappingFileName = None - climatologyFileName = '{}/{}_{}_{}.nc'.format( - climatologyDirectory, fieldName, obsGridName, monthNames) - regriddedFileName = '{}/{}_{}_to_{}_{}.nc'.format( - regriddedDirectory, fieldName, obsGridName, comparisonGridName, - monthNames) + else: + # we need to build the path to the mapping file and an + # appropriate file name + mappingSubdirectory = \ + build_config_full_path(config, 'output', + 'mappingSubdirectory') + + make_directories(mappingSubdirectory) + + mappingFileName = '{}/{}_{}_to_{}_{}.nc'.format( + mappingSubdirectory, mappingFilePrefix, + self.sourceDescriptor.meshName, + self.comparisonDescriptor.meshName, + method) + + config.set(mappingFileSection, mappingFileOption, + mappingFileName) + + remapper = Remapper(self.sourceDescriptor, self.comparisonDescriptor, + mappingFileName) + + remapper.build_mapping_file(method=method) + + self.remapper = remapper + return remapper # }}} + + def compute(self, ds, monthValues=None, maskVaries=True): # {{{ + """ + Compute a monthly, seasonal or annual climatology data set from a data + set. The mean is weighted but the number of days in each month of + the data set, ignoring values masked out with NaNs. If the month + coordinate is not present, a data array ``month`` will be added based + on ``Time`` and the provided calendar. + + Parameters + ---------- + ds : ``xarray.Dataset`` or ``xarray.DataArray`` object + A data set with a ``Time`` coordinate expressed as days since + 0001-01-01 or ``month`` coordinate + + monthValues : int or array-like of ints, optional + A single month or an array of months to be averaged together. + If this option is not provided, the value of ``monthValues`` passed + to ``__init__`` will be used. + + maskVaries: bool, optional + If the mask (where variables in ``ds`` are ``NaN``) varies with + time. If not, the weighted average does not need make extra effort + to account for the mask. Most MPAS fields will have masks that + don't vary in time, whereas observations may sometimes be present + only at some times and not at others, requiring + ``maskVaries = True``. + + Returns + ------- + dataSet : object of same type as ``ds`` + A data set without the ``'Time'`` coordinate containing the mean + of ds over all months in monthValues, weighted by the number of + days in each month. Also stored as the ``dataSet`` attribute. + + Authors + ------- + Xylar Asay-Davis + """ + + if monthValues is None: + monthValues = self.monthValues + + ds = add_years_months_days_in_month(ds, self.calendar) + + mask = xr.zeros_like(ds.month, bool) + + for month in monthValues: + mask = xr.ufuncs.logical_or(mask, ds.month == month) + + climatologyMonths = ds.where(mask, drop=True) + + self.dataSet = self._compute_masked_mean(climatologyMonths, maskVaries) + + return self.dataSet # }}} + + def compute_monthly(self, ds, maskVaries=True): # {{{ + """ + Compute monthly climatologies from a data set. The mean is + weighted by the number of days in each month of the data set, + ignoring values masked out with NaNs. If the month coordinate is + not present, a data array ``month`` will be added based on + ``Time``. + + Parameters + ---------- + ds : ``xarray.Dataset`` or ``xarray.DataArray`` object + A data set with a ``Time`` coordinate expressed as days since + 0001-01-01 or ``month`` coordinate + + maskVaries: bool, optional + If the mask (where variables in ``ds`` are ``NaN``) varies with + time. If not, the weighted average does not need make extra + effort to account for the mask. Most MPAS fields will have + masks that don't vary in time, whereas observations may + sometimes be present only at some times and not at others, + requiring ``maskVaries = True``. + + Returns + ------- + dataSet : object of same type as ``ds`` + A data set with the same fields at ``ds`` but in which all data + from each month has been averaged together to create a + climatology for each month. Also stored as the ``dataSet`` + attribute + + Authors + ------- + Xylar Asay-Davis + """ + + def compute_one_month_climatology(ds): + monthValues = list(ds.month.values) + return self.compute(ds, monthValues, maskVaries) + + ds = add_years_months_days_in_month(ds, self.calendar) + + self.dataSet = \ + ds.groupby('month').apply(compute_one_month_climatology) + + return self.dataSet # }}} + + def remap_and_write(self, useNcremap=None): # {{{ + """ + Remap the climatology data set produced by a call to ``compute`` or + ``cache``, write the result to an output file, and return the remapped + data set. + + This call requires that ``climatologyFileName`` and + ``remappedFileName`` attributes be assigned, which happens + automatically for some of the subclasses of ``Climatology`` + (e.g. ``MpasClimatology`` and ``ObservationClimatology``). + + Note that the files named by attributes ``climatologyFileName`` and + ``remappedFileName`` will be overwritten if they exist, so if + this behavior is not desired, the calling code should skip this call + if the files exist and simply load the contents of + ``remappedFileName``. + + Parameters + ---------- + useNcremap : bool, optional + If present, overrides ``useNcremap`` from the config file. This + is useful for data sets that cannot be handle by ncremap (or are + masked better with the "online" remapper). + + Returns + ------- + remappedDataSet : ``xarray.DataSet`` or ``xarray.DataArray`` object + A data set containing the remapped climatology. Also stored in + the ``remappedDataSet`` attribute. + + Authors + ------- + Xylar Asay-Davis + """ + if useNcremap is None: + useNcremap = self.config.getboolean('climatology', 'useNcremap') + + if self.remapper.mappingFileName is None: + # no remapping is needed + self.remappedDataSet = self.dataSet + else: + if useNcremap: + if not os.path.exists(self.climatologyFileName): + self.dataSet.to_netcdf(self.climatologyFileName) + self.remapper.remap_file(inFileName=self.climatologyFileName, + outFileName=self.remappedFileName, + overwrite=True) + self.remappedDataSet = xr.open_dataset(self.remappedFileName) + else: + renormalizationThreshold = self.config.getfloat( + 'climatology', 'renormalizationThreshold') + + self.remappedDataSet = self.remapper.remap( + self.dataSet, renormalizationThreshold) + self.remappedDataSet.to_netcdf(self.remappedFileName) + return self.remappedDataSet # }}} + + def _compute_masked_mean(self, ds, maskVaries): # {{{ + ''' + Compute the time average of data set, masked out where the variables + in ds are NaN and, if ``maskVaries == True``, weighting by the number + of days used to compute each monthly mean time in ds. + + Authors + ------- + Xylar Asay-Davis + ''' + def ds_to_weights(ds): + # make an identical data set to ds but replacing all data arrays + # with notnull applied to that data array + weights = ds.copy(deep=True) + if isinstance(ds, xr.core.dataarray.DataArray): + weights = ds.notnull() + elif isinstance(ds, xr.core.dataset.Dataset): + for var in ds.data_vars: + weights[var] = ds[var].notnull() + else: + raise TypeError('ds must be an instance of either ' + 'xarray.Dataset or xarray.DataArray.') - make_directories(climatologyDirectory) + return weights - if not _matches_comparison(remapper.sourceDescriptor, - remapper.destinationDescriptor): - make_directories(regriddedDirectory) + if maskVaries: + dsWeightedSum = (ds * ds.daysInMonth).sum(dim='Time', + keep_attrs=True) - return (climatologyFileName, regriddedFileName) # }}} + weights = ds_to_weights(ds) + weightSum = (weights * ds.daysInMonth).sum(dim='Time') -def compute_monthly_climatology(ds, calendar=None, maskVaries=True): # {{{ - """ - Compute monthly climatologies from a data set. The mean is weighted but - the number of days in each month of the data set, ignoring values masked - out with NaNs. If the month coordinate is not present, a data array - ``month`` will be added based on ``Time`` and the provided calendar. + timeMean = dsWeightedSum / weightSum.where(weightSum > 0.) + else: + days = ds.daysInMonth.sum(dim='Time') + + dsWeightedSum = (ds * ds.daysInMonth).sum(dim='Time', + keep_attrs=True) + + timeMean = dsWeightedSum / days.where(days > 0.) + + return timeMean # }}} + + def _matches_comparison(self, obsDescriptor, comparisonDescriptor): # {{{ + ''' + Determine if the two meshes are the same + + Authors + ------- + Xylar Asay-Davis + ''' + + if isinstance(obsDescriptor, ProjectionGridDescriptor) and \ + isinstance(comparisonDescriptor, ProjectionGridDescriptor): + # pretty hard to determine if projections are the same, so we'll + # rely on the grid names + match = \ + obsDescriptor.meshName == comparisonDescriptor.meshName and \ + len(obsDescriptor.x) == len(comparisonDescriptor.x) and \ + len(obsDescriptor.y) == len(comparisonDescriptor.y) and \ + numpy.all(numpy.isclose(obsDescriptor.x, + comparisonDescriptor.x)) and \ + numpy.all(numpy.isclose(obsDescriptor.y, + comparisonDescriptor.y)) + elif isinstance(obsDescriptor, LatLonGridDescriptor) and \ + isinstance(comparisonDescriptor, LatLonGridDescriptor): + match = \ + ((('degree' in obsDescriptor.units and + 'degree' in comparisonDescriptor.units) or + ('radian' in obsDescriptor.units and + 'radian' in comparisonDescriptor.units)) and + len(obsDescriptor.lat) == len(comparisonDescriptor.lat) and + len(obsDescriptor.lon) == len(comparisonDescriptor.lon) and + numpy.all(numpy.isclose(obsDescriptor.lat, + comparisonDescriptor.lat)) and + numpy.all(numpy.isclose(obsDescriptor.lon, + comparisonDescriptor.lon))) + else: + match = False - Parameters - ---------- - ds : ``xarray.Dataset`` or ``xarray.DataArray`` object - A data set with a ``Time`` coordinate expressed as days since - 0001-01-01 or ``month`` coordinate - - calendar : ``{'gregorian', 'gregorian_noleap'}``, optional - The name of one of the calendars supported by MPAS cores, used to - determine ``month`` from ``Time`` coordinate, so must be supplied if - ``ds`` does not already have a ``month`` coordinate or data array - - maskVaries: bool, optional - If the mask (where variables in ``ds`` are ``NaN``) varies with time. - If not, the weighted average does not need make extra effort to account - for the mask. Most MPAS fields will have masks that don't vary in - time, whereas observations may sometimes be present only at some - times and not at others, requiring ``maskVaries = True``. + return match # }}} - Returns - ------- - climatology : object of same type as ``ds`` - A data set without the ``'Time'`` coordinate containing the mean - of ds over all months in monthValues, weighted by the number of days - in each month. +# }}} - Authors - ------- - Xylar Asay-Davis - Last Modified - ------------- - 04/08/2017 +class MpasClimatology(Climatology): # {{{ """ + A class for computing climatologies from an MPAS monthly mean data set. - def compute_one_month_climatology(ds): - monthValues = list(ds.month.values) - return compute_climatology(ds, monthValues, calendar, maskVaries) - - ds = add_years_months_days_in_month(ds, calendar) + Attributes + ---------- + fieldName : str + The name of the field for which the climatology is being computed, + used in constructing the names of cache files. - monthlyClimatology = \ - ds.groupby('month').apply(compute_one_month_climatology) + streamName : str + The name of the stream from which the climatology data set will + be read, used to cache times and update their bounds in the + configuration parser. - return monthlyClimatology # }}} + startYear, endYear : int + The start and end years of the climatology + climatologyDirectory : str + The directory where climatologies will be cached -def compute_climatology(ds, monthValues, calendar=None, - maskVaries=True): # {{{ - """ - Compute a monthly, seasonal or annual climatology data set from a data - set. The mean is weighted but the number of days in each month of - the data set, ignoring values masked out with NaNs. If the month - coordinate is not present, a data array ``month`` will be added based - on ``Time`` and the provided calendar. + climatologyPrefix : str + The prefix (including full path) on climatology cache files. A call + to the ``cache`` method may produce several intermediate cache files + as well as a "final" cache file with the seasonal average between + ``startYear`` and ``endYear``. - Parameters - ---------- - ds : ``xarray.Dataset`` or ``xarray.DataArray`` object - A data set with a ``Time`` coordinate expressed as days since - 0001-01-01 or ``month`` coordinate + climatologyFileName : str + The full path to final climatology cache file where ``dataSet`` + (the climatology data set before remapping) should be stored. - monthValues : int or array-like of ints - A single month or an array of months to be averaged together + remappedFileName : str + The name of the cache file where ``remappedDataSet`` (the climatology + data set after remapping) should be stored. This attribute will be + set only if a comparison grid name is supplied. - calendar : ``{'gregorian', 'gregorian_noleap'}``, optional - The name of one of the calendars supported by MPAS cores, used to - determine ``month`` from ``Time`` coordinate, so must be supplied if - ``ds`` does not already have a ``month`` coordinate or data array + sourceDescriptor : ``MpasMeshDescriptor`` object + A descriptor of the source MPAS mesh used for remapping. - maskVaries: bool, optional - If the mask (where variables in ``ds`` are ``NaN``) varies with time. - If not, the weighted average does not need make extra effort to account - for the mask. Most MPAS fields will have masks that don't vary in - time, whereas observations may sometimes be present only at some - times and not at others, requiring ``maskVaries = True``. + remapper : ``Remapper`` object + A remapper between the source and comparison grids, created + automatically if a comparison grid is supplied. + + dataSet : ``xarray.Dataset`` or ``xarray.DataArray`` object + A data set containing the climatology (before remapping, if + applicable). Available only after calling the ``compute``, + ``compute_monthly``, or ``cache`` methods. + + remappedDataSet : ``xarray.Dataset`` or ``xarray.DataArray`` object + A data set containing the remapped climatology. Loaded automatically + if a comparison grid is supplied and the cache file in + ``remappedFileName`` already exists. Otherwise, available only after + calling the ``remap_and_write`` method. + + Examples + -------- + Here is an example of how to use an ``MpasClimatology`` object within a + task. In the example, ``self._open_mpas_dataset_part`` is a function with + arguments ``inputFileNames``, ``startDate``, and ``endDate`` used to open + and return part of the data set during caching. The result of the example + code is that the xarray data set ``climatology.remappedDataSet`` is + available for computation and plotting, having either been read from a + cache file or computed for the MPAS monthly time series data set. + + >>> from ..shared.climatology import MpasClimatology + >>> restartFileName = self.runStreams.readpath('restart')[0] + >>> climatology = MpasClimatology( + task=self, + fieldName='sst', + monthNames='ANN', + streamName='timeSeriesStats', + meshFileName=restartFileName, + comparisonGrid='latlon', + mappingFileSection='climatology', + mappingFileOption='mpasMappingFile', + mappingFilePrefix='map', + method='bilinear') + >>> if climatology.remappedDataSet is None: + climatology.cache( + openDataSetFunc=self._open_mpas_dataset_part, + printProgress=True) + climatology.remap_and_write() + >>> print climatology.remappedDataSet + + Here is another example, this time without remapping: + + >>> from ..shared.climatology import MpasClimatology + >>> climatology = MpasClimatology( + task=self, + fieldName='sst', + monthNames='ANN', + streamName='timeSeriesStats') + >>> climatology.cache( + openDataSetFunc=self._open_mpas_dataset_part, + printProgress=True) + >>> print climatology.dataSet - Returns - ------- - climatology : object of same type as ``ds`` - A data set without the ``'Time'`` coordinate containing the mean - of ds over all months in monthValues, weighted by the number of days - in each month. Authors ------- Xylar Asay-Davis - - Last Modified - ------------- - 04/08/2017 """ - ds = add_years_months_days_in_month(ds, calendar) - - mask = xr.zeros_like(ds.month, bool) - - for month in monthValues: - mask = xr.ufuncs.logical_or(mask, ds.month == month) - - climatologyMonths = ds.where(mask, drop=True) - - climatology = _compute_masked_mean(climatologyMonths, maskVaries) - - return climatology # }}} - - -def cache_climatologies(ds, monthValues, config, cachePrefix, calendar, - printProgress=False): # {{{ - ''' - Cache NetCDF files for each year of an annual climatology, and then use - the cached files to compute a climatology for the full range of years. - The start and end years of the climatology are taken from ``config``, and - are updated in ``config`` if the data set ``ds`` doesn't contain this - full range. - - Note: only works with climatologies where the mask (locations of ``NaN`` - values) doesn't vary with time. - - Parameters - ---------- - ds : ``xarray.Dataset`` or ``xarray.DataArray`` object - A data set with a ``Time`` coordinate expressed as days since - 0001-01-01 + def __init__(self, task, fieldName, monthNames, streamName, + meshFileName=None, comparisonGrid=None, + mappingFileSection=None, mappingFileOption=None, + mappingFilePrefix=None, method=None): # {{{ + """ + Create a new object for computing climatologies from an MPAS monthly + mean data set. - monthValues : int or array-like of ints - A single month or an array of months to be averaged together + If a ``comparisonGrid`` is provided, the climatology object also + either loads data from a cache file containing remapped data (if + one already exists) or creates a remapper that can be used later + on (via the ``remap_and_write`` method) to compute the remapped + climatology. - config : instance of MpasAnalysisConfigParser - Contains configuration options + Parameters + ---------- + task : ``AnalysisTask`` object + An analysis task for which the climatology is needed. ``task`` is + used to get config options, define the component name, map streams, + namelists and variables, etc. - cachePrefix : str - The file prefix (including path) to which the year (or years) will be - appended as cache files are stored + fieldName : str + The name of the field for which the climatology is being computed, + used in constructing the names of cache files. - calendar : ``{'gregorian', 'gregorian_noleap'}`` - The name of one of the calendars supported by MPAS cores, used to - determine ``year`` and ``month`` from ``Time`` coordinate + monthNames : str + The months that make up the climatology, used in constructing the + names of cache files. - printProgress: bool, optional - Whether progress messages should be printed as the climatology is - computed + streamName : str + The name of the stream from which the climatology data set will + be read, used to cache times and update their bounds in the + configuration parser. - Returns - ------- - climatology : object of same type as ``ds`` - A data set without the ``'Time'`` coordinate containing the mean - of ds over all months in monthValues, weighted by the number of days - in each month. + The remaining parameters are all required if remapping is to be + performed. - Authors - ------- - Xylar Asay-Davis + meshFileName : str, optional + The name of the MPAS mesh file, used to create a descriptor of the + MPAS mesh for remapping. - Last Modified - ------------- - 04/11/2017 - ''' - startYearClimo = config.getint('climatology', 'startYear') - endYearClimo = config.getint('climatology', 'endYear') - yearsPerCacheFile = config.getint('climatology', 'yearsPerCacheFile') + comparisonGrid : {'latlon'}, optional + The name of the comparison grid to use for remapping (if any). - if printProgress: - print ' Computing and caching climatologies covering {}-year ' \ - 'spans...'.format(yearsPerCacheFile) + mappingFileSection, mappingFileOption : str, optional + Section and option in ``config`` where the name of the mapping file + may be given, or where it will be stored if a new mapping file is + created - ds = add_years_months_days_in_month(ds, calendar) + mappingFilePrefix : str, optional + A prefix to be prepended to the mapping file name - cacheInfo, cacheIndices = _setup_climatology_caching(ds, startYearClimo, - endYearClimo, - yearsPerCacheFile, - cachePrefix, - monthValues) + method : {'bilinear', 'neareststod', 'conserve'}, optional + The method of interpolation used. - ds = ds.copy() - ds.coords['cacheIndices'] = ('Time', cacheIndices) + Raises + ------ + ValueError + If comarisonGrid does not describe a known comparions grid - # compute and store each cache file with interval yearsPerCacheFile - _cache_individual_climatologies(ds, cacheInfo, printProgress, - yearsPerCacheFile, monthValues, - calendar) + Authors + ------- + Xylar Asay-Davis + """ - # compute the aggregate climatology - climatology = _cache_aggregated_climatology(startYearClimo, endYearClimo, - cachePrefix, printProgress, - monthValues, cacheInfo) + super(MpasClimatology, self).__init__(task, monthNames, comparisonGrid) - return climatology # }}} + self.fieldName = fieldName + self.streamName = streamName + climSection = 'climatology' + task.update_start_end_date(section=climSection, + streamName=streamName) + self.startYear = self.config.getint(climSection, 'startYear') + self.endYear = self.config.getint(climSection, 'endYear') -def update_start_end_year(ds, config, calendar): # {{{ - """ - Given a monthly climatology, compute a seasonal climatology weighted by - the number of days in each month (on the no-leap-year calendar). + mpasMeshName = self.config.get('input', 'mpasMeshName') - Parameters - ---------- - ds : instance of xarray.Dataset - A data set from which start and end years will be determined + self.climatologyDirectory = build_config_full_path( + self.config, 'output', 'mpasClimatologySubdirectory') - config : instance of MpasAnalysisConfigParser - Contains configuration options + make_directories(self.climatologyDirectory) - calendar : {'gregorian', 'gregorian_noleap'} - The name of one of the calendars supported by MPAS cores + self.climatologyPrefix = \ + '{}/{}_{}_{}'.format(self.climatologyDirectory, fieldName, + mpasMeshName, monthNames) - Returns - ------- - changed : bool - Whether the start and end years were changed + yearString, fileSuffix = self._get_year_string(self.startYear, + self.endYear) + self.climatologyFileName = \ + '{}_{}.nc'.format(self.climatologyPrefix, fileSuffix) - startYear, endYear : int - The start and end years of the data set - - Authors - ------- - Xylar Asay-Davis - - Last Modified - ------------- - 03/25/2017 - """ - requestedStartYear = config.getint('climatology', 'startYear') - requestedEndYear = config.getint('climatology', 'endYear') - - startYear = days_to_datetime(ds.Time.min().values, calendar=calendar).year - endYear = days_to_datetime(ds.Time.max().values, calendar=calendar).year - changed = False - if startYear != requestedStartYear or endYear != requestedEndYear: - message = "climatology start and/or end year different from " \ - "requested\n" \ - "requestd: {:04d}-{:04d}\n" \ - "actual: {:04d}-{:04d}\n".format(requestedStartYear, - requestedEndYear, - startYear, - endYear) - warnings.warn(message) - config.set('climatology', 'startYear', str(startYear)) - config.set('climatology', 'endYear', str(endYear)) - changed = True - - return changed, startYear, endYear # }}} - - -def add_years_months_days_in_month(ds, calendar=None): # {{{ - ''' - Add ``year``, ``month`` and ``daysInMonth`` as data arrays in ``ds``. - The number of days in each month of ``ds`` is computed either using the - ``startTime`` and ``endTime`` if available or assuming ``gregorian_noleap`` - calendar and ignoring leap years. ``year`` and ``month`` are computed - accounting correctly for the the calendar. - - Parameters - ---------- - ds : ``xarray.Dataset`` or ``xarray.DataArray`` object - A data set with a ``Time`` coordinate expressed as days since - 0001-01-01 + if comparisonGrid is not None: + remappedDirectory = build_config_full_path( + self.config, 'output', 'mpasRemappedClimSubdirectory') - calendar : ``{'gregorian', 'gregorian_noleap'}``, optional - The name of one of the calendars supported by MPAS cores, used to - determine ``year`` and ``month`` from ``Time`` coordinate - - Returns - ------- - ds : object of same type as ``ds`` - The data set with ``year``, ``month`` and ``daysInMonth`` data arrays - added (if not already present) + make_directories(remappedDirectory) - Authors - ------- - Xylar Asay-Davis - - Last Modified - ------------- - 04/08/2017 - ''' + self.remappedFileName = '{}/{}_{}_to_{}_{}_{}.nc'.format( + remappedDirectory, fieldName, mpasMeshName, + self.comparisonDescriptor.meshName, monthNames, fileSuffix) + + if os.path.exists(self.remappedFileName): + # no need to create the remapper, since the cached data set + # alredy exists + self.remappedDataSet = xr.open_dataset(self.remappedFileName) + else: + self.sourceDescriptor = \ + MpasMeshDescriptor(fileName=meshFileName, + meshName=mpasMeshName) + + self.create_remapper(mappingFileSection, mappingFileOption, + mappingFilePrefix, method) + + self.remappedDataSet = None + # }}} + + def cache(self, openDataSetFunc, printProgress=False): # {{{ + ''' + Cache NetCDF files for each year of an annual climatology, and then use + the cached files to compute a climatology for the full range of years. + The start and end years of the climatology are taken from ``config``, + and are updated in ``config`` if the data set ``ds`` doesn't contain + this full range. + + Note: only works with climatologies where the mask (locations of + ``NaN`` values) doesn't vary with time. + + Parameters + ---------- + openDataSetFunc : function + A function with arguments ``inputFileNames``, ``startDate``, + and ``endDate`` used to open a portion of the data set from which + the climatology is computed and cached. Typically, the function + makes a call to ``generalized_reader.open_multifile_dataset``, + perhaps performing further manipulation of the resulting data set. + + printProgress: bool, optional + Whether progress messages should be printed as the climatology is + computed + + Returns + ------- + dataSet : object of same type as ``ds`` + A data set without the ``'Time'`` coordinate containing the mean + of ds over all months in monthValues, weighted by the number of + days in each month. Also avialable through the ``dataSet`` + attribute. + + Authors + ------- + Xylar Asay-Davis + ''' + cacheInfo = self._setup_climatology_caching(printProgress) + + # compute and store each cache file with interval yearsPerCacheFile + self._cache_individual_climatologies(openDataSetFunc, cacheInfo, + printProgress) + + # compute the aggregate climatology + self.dataSet = self._cache_aggregated_climatology(cacheInfo, + printProgress) + + return self.dataSet # }}} + + def _setup_climatology_caching(self, printProgress): # {{{ + ''' + Determine which cache files already exist, which are incomplete and + which years are present in each cache file (whether existing or to be + created). + + Authors + ------- + Xylar Asay-Davis + ''' + yearsPerCacheFile = self.config.getint('climatology', + 'yearsPerCacheFile') + startDate = self.config.get('climatology', 'startDate') + endDate = self.config.get('climatology', 'endDate') + startDate = string_to_days_since_date(dateString=startDate, + calendar=self.calendar) + endDate = string_to_days_since_date(dateString=endDate, + calendar=self.calendar) + + inFileNames = self.task.get_input_file_names( + self.streamName, startAndEndDateSection='climatology') + + fullTimeCache = self.task.cache_multifile_dataset_times( + inFileNames, self.streamName, timeVariableName='Time') + + # find only those cached times between starDate and endDate + timeCache = OrderedDict() + times = [] + for fileName in fullTimeCache: + localTimes = fullTimeCache[fileName]['times'] + mask = numpy.logical_and(localTimes >= startDate, + localTimes < endDate) + if numpy.count_nonzero(mask) == 0: + continue + + times.extend(list(localTimes[mask])) + timeCache[fileName] = fullTimeCache[fileName] - if ('year' in ds.coords and 'month' in ds.coords and - 'daysInMonth' in ds.coords): - return ds + if printProgress: + print ' Computing and caching climatologies covering {}-year ' \ + 'spans...'.format(yearsPerCacheFile) + + cacheInfo = [] + + firstFile = None + lastFile = None + # figure out which files to load and which years go in each file + for firstYear in range(self.startYear, self.endYear+1, + yearsPerCacheFile): + years = numpy.arange(firstYear, firstYear+yearsPerCacheFile) + + monthsIfDone = len(self.monthValues)*len(years) + + yearString, fileSuffix = self._get_year_string(years[0], years[-1]) + outputFileName = '{}_{}.nc'.format(self.climatologyPrefix, + fileSuffix) + + done = False + if os.path.exists(outputFileName): + # already cached + dsCached = None + try: + dsCached = xr.open_dataset(outputFileName) + except IOError: + # assuming the cache file is corrupt, so deleting it. + message = 'Deleting cache file {}, which appears to have' \ + ' been corrupted.'.format(outputFileName) + warnings.warn(message) + os.remove(outputFileName) + + if ((dsCached is not None) and + (dsCached.attrs['totalMonths'] == monthsIfDone)): + # also complete, so we can move on + done = True + if dsCached is not None: + dsCached.close() + + inputFileNames = [] + for year in years: + for month in self.monthValues: + for fileName in timeCache: + timeCacheYears = timeCache[fileName]['years'] + timeCacheMonths = timeCache[fileName]['months'] + mask = numpy.logical_and(timeCacheYears == year, + timeCacheMonths == month) + if numpy.count_nonzero(mask) > 0: + inputFileNames.append(fileName) + + if len(inputFileNames) > 0: + cacheDict = {'outputFileName': outputFileName, + 'done': done, + 'years': years, + 'yearString': yearString, + 'inputFileNames': inputFileNames} + cacheInfo.append(cacheDict) + if not done: + lastFile = inputFileNames[-1] + if firstFile is None: + firstFile = inputFileNames[0] + + if printProgress and firstFile is not None: + print '\n Caching data from files:\n' \ + ' {} through\n {}'.format( + os.path.basename(firstFile), + os.path.basename(lastFile)) + + return cacheInfo # }}} + + def _cache_individual_climatologies(self, openDataSetFunc, cacheInfo, + printProgress): # {{{ + ''' + Cache individual climatologies for later aggregation. + + Authors + ------- + Xylar Asay-Davis + ''' + + startDate = self.config.get('climatology', 'startDate') + endDate = self.config.get('climatology', 'endDate') + + for info in cacheInfo: + if info['done']: + continue + outputFileName = info['outputFileName'] + yearString = info['yearString'] + years = info['years'] + inputFileNames = info['inputFileNames'] + ds = openDataSetFunc(inputFileNames, startDate, endDate) + ds = add_years_months_days_in_month(ds, self.calendar) + mask = numpy.zeros(ds.dims['Time'], bool) + for year in years: + mask = numpy.logical_or(mask, ds.year.values == year) + ds.coords['cacheMask'] = ('Time', mask) + dsYear = ds.where(ds.cacheMask, drop=True) + + if printProgress: + print ' {}'.format(yearString) + + totalDays = dsYear.daysInMonth.sum(dim='Time').values + + monthCount = dsYear.dims['Time'] + + climatology = self.compute(dsYear, maskVaries=False) + + climatology.attrs['totalDays'] = totalDays + climatology.attrs['totalMonths'] = monthCount + climatology.attrs['fingerprintClimo'] = fingerprint_generator() + + climatology.to_netcdf(outputFileName) + climatology.close() + + # }}} + + def _cache_aggregated_climatology(self, cacheInfo, printProgress): # {{{ + ''' + Cache aggregated climatology from individual climatologies. + + Authors + ------- + Xylar Asay-Davis + + ''' + yearString, fileSuffix = self._get_year_string(self.startYear, + self.endYear) + outputFileName = '{}_{}.nc'.format(self.climatologyPrefix, fileSuffix) - ds = ds.copy() + done = False + if len(cacheInfo) == 0: + climatology = None + done = True - if 'year' not in ds.coords or 'month' not in ds.coords: - if calendar is None: - raise ValueError('calendar must be provided if month and year ' - 'coordinate is not in ds.') - datetimes = days_to_datetime(ds.Time, calendar=calendar) + if os.path.exists(outputFileName): + # already cached + climatology = None + try: + climatology = xr.open_dataset(outputFileName) - if 'year' not in ds.coords: - ds.coords['year'] = ('Time', [date.year for date in datetimes]) + except IOError: + # assuming the cache file is corrupt, so deleting it. + message = 'Deleting cache file {}, which appears to have ' \ + 'been corrupted.'.format(outputFileName) + warnings.warn(message) + os.remove(outputFileName) - if 'month' not in ds.coords: - ds.coords['month'] = ('Time', [date.month for date in datetimes]) + if len(cacheInfo) == 1 and \ + outputFileName == cacheInfo[0]['outputFileName']: + # theres only one cache file and it already has the same name + # as the aggregated file so no need to aggregate + done = True - if 'daysInMonth' not in ds.coords: - if 'startTime' in ds.coords and 'endTime' in ds.coords: - ds.coords['daysInMonth'] = ds.endTime - ds.startTime + elif climatology is not None: + monthsIfDone = (self.endYear-self.startYear+1) * \ + len(self.monthValues) + if climatology.attrs['totalMonths'] == monthsIfDone: + # also complete, so we can move on + done = True + else: + climatology.close() + + if not done: + if printProgress: + print ' Computing aggregated climatology ' \ + '{}...'.format(yearString) + + first = True + for info in cacheInfo: + inputFileName = info['outputFileName'] + ds = xr.open_dataset(inputFileName) + days = ds.attrs['totalDays'] + months = ds.attrs['totalMonths'] + if first: + totalDays = days + totalMonths = months + climatology = ds * days + first = False + else: + totalDays += days + totalMonths += months + climatology = climatology + ds * days + + ds.close() + climatology = climatology / totalDays + + climatology.attrs['totalDays'] = totalDays + climatology.attrs['totalMonths'] = totalMonths + climatology.attrs['fingerprintClimo'] = fingerprint_generator() + + climatology.to_netcdf(outputFileName) + + return climatology # }}} + + def _get_year_string(self, startYear, endYear): # {{{ + if startYear == endYear: + yearString = '{:04d}'.format(startYear) + fileSuffix = 'year{}'.format(yearString) else: - if calendar == 'gregorian': - message = 'The MPAS run used the Gregorian calendar but ' \ - 'does not appear to have\n' \ - 'supplied start and end times. Climatologies ' \ - 'will be computed with\n' \ - 'month durations ignoring leap years.' - warnings.warn(message) + yearString = '{:04d}-{:04d}'.format(startYear, endYear) + fileSuffix = 'years{}'.format(yearString) - daysInMonth = numpy.array([constants.daysInMonth[month-1] for - month in ds.month.values], float) - ds.coords['daysInMonth'] = ('Time', daysInMonth) + return yearString, fileSuffix # }}} - return ds # }}} + # }}} -def remap_and_write_climatology(config, climatologyDataSet, - climatologyFileName, regriddedFileName, - remapper): # {{{ +class ObservationClimatology(Climatology): # {{{ """ - Given a field in a climatology data set, use the ``remapper`` to regrid - horizontal dimensions of all fields, write the results to an output file, - and return the regridded data set. + A class for computing climatologies from an obsevational data set. - Note that ``climatologyFileName`` and ``regriddedFileName`` will be - overwritten if they exist, so if this behavior is not desired, the calling - code should skip this call if the files exist and simply load the contents - of ``regriddedFileName``. - - Parameters + Attributes ---------- - config : instance of ``MpasAnalysisConfigParser`` - Contains configuration options - - climatologyDataSet : ``xarray.DataSet`` or ``xarray.DataArray`` object - A data set containing a climatology - fieldName : str - A field within the climatology to be remapped + The name of the field for which the climatology is being computed, + used in constructing the names of cache files. + + climatologyDirectory : str + The directory where climatologies will be cached climatologyFileName : str - The name of the output file to which the data set should be written - before regridding (if using ncremap). + The full path to final climatology cache file where ``dataSet`` + (the climatology data set before remapping) should be stored. + + remappedFileName : str + The name of the cache file where ``remappedDataSet`` (the climatology + data set after remapping) should be stored. This attribute will be + set only if a comparison grid name is supplied. - regriddedFileName : str - The name of the output file to which the regridded data set should - be written. + sourceDescriptor : ``MeshDescriptor`` object + The descriptor of the input data grid used for remapping. remapper : ``Remapper`` object - A remapper that can be used to remap files or data sets to a - comparison grid. + A remapper between the source and comparison grids, created + automatically if a comparison grid is supplied. + + dataSet : ``xarray.Dataset`` or ``xarray.DataArray`` object + A data set containing the climatology (before remapping, if + applicable). Available only after calling the ```compute`` or + ``compute_monthly`` methods. + + remappedDataSet : ``xarray.Dataset`` or ``xarray.DataArray`` object + A data set containing the remapped climatology. Loaded automatically + if a comparison grid is supplied and the cache file in + ``remappedFileName`` already exists. Otherwise, available only after + calling the ``remap_and_write`` method. + + Examples + -------- + Here is an example of how to use an ``ObservationClimatology`` object + within a task. In the example, ``obsDescriptor`` is a ``MeshDescriptor`` + object describing the input observations grid and ``dsObs`` is an xarray + data set containing the monthly mean time series of the observations. The + result of the example code is that the xarray data set + ``climatology.remappedDataSet`` is available for computation and plotting, + having either been read from a cache file or computed for the monthly + time series data set. + + >>> from ..shared.climatology import ObservationClimatology + >>> climatology = ObservationClimatology( + task=self, + fieldName='sst', + monthNames='ANN', + obsGridDescriptor=obsDescriptor, + comparisonGrid='latlon', + mappingFileSection='oceanObservations', + mappingFileOption='sstClimatologyMappingFile', + mappingFilePrefix='map_obs_sst', + method='bilinear') + >>> if climatology.remappedDataSet is None: + climatology.compute(ds=dsObs) + climatology.remap_and_write() + >>> print climatology.remappedDataSet - Returns - ------- - remappedClimatology : ``xarray.DataSet`` or ``xarray.DataArray`` object - A data set containing the remapped climatology Authors ------- Xylar Asay-Davis - - Last Modified - ------------- - 04/13/2017 """ - useNcremap = config.getboolean('climatology', 'useNcremap') - - if remapper.mappingFileName is None: - # no remapping is needed - remappedClimatology = climatologyDataSet - else: - if useNcremap: - if not os.path.exists(climatologyFileName): - climatologyDataSet.to_netcdf(climatologyFileName) - remapper.remap_file(inFileName=climatologyFileName, - outFileName=regriddedFileName, - overwrite=True) - remappedClimatology = xr.open_dataset(regriddedFileName) - else: - renormalizationThreshold = config.getfloat( - 'climatology', 'renormalizationThreshold') - - remappedClimatology = remapper.remap(climatologyDataSet, - renormalizationThreshold) - remappedClimatology.to_netcdf(regriddedFileName) - return remappedClimatology # }}} - - -def _compute_masked_mean(ds, maskVaries): # {{{ - ''' - Compute the time average of data set, masked out where the variables in ds - are NaN and, if ``maskVaries == True``, weighting by the number of days - used to compute each monthly mean time in ds. - - Authors - ------- - Xylar Asay-Davis - Last Modified - ------------- - 04/08/2017 - ''' - def ds_to_weights(ds): - # make an identical data set to ds but replacing all data arrays with - # nonnull applied to that data array - weights = ds.copy(deep=True) - if isinstance(ds, xr.core.dataarray.DataArray): - weights = ds.notnull() - elif isinstance(ds, xr.core.dataset.Dataset): - for var in ds.data_vars: - weights[var] = ds[var].notnull() - else: - raise TypeError('ds must be an instance of either xarray.Dataset ' - 'or xarray.DataArray.') + def __init__(self, task, fieldName, monthNames, + obsGridDescriptor, comparisonGrid=None, + mappingFileSection=None, mappingFileOption=None, + mappingFilePrefix=None, method=None): # {{{ + """ + Create a new object for creating climatologies from observational data + sets. - return weights + Parameters + ---------- + task : ``AnalysisTask`` object + An analysis task for which the climatology is needed. ``task`` is + used to get config options, define the component name, map streams, + namelists and variables, etc. - if maskVaries: - dsWeightedSum = (ds * ds.daysInMonth).sum(dim='Time', keep_attrs=True) + fieldName : str + The name of the field for which the climatology is being computed, + used in constructing the names of cache files. - weights = ds_to_weights(ds) + monthNames : str + The months that make up the climatology, used in constructing the + names of cache files. - weightSum = (weights * ds.daysInMonth).sum(dim='Time') + obsGridDescriptor : ``MeshDescriptor`` object + The descriptor of the input data grid - timeMean = dsWeightedSum / weightSum.where(weightSum > 0.) - else: - days = ds.daysInMonth.sum(dim='Time') + The remaining parameters are all required if remapping is to be + performed. - dsWeightedSum = (ds * ds.daysInMonth).sum(dim='Time', keep_attrs=True) + comparisonGrid : {'latlon'}, optional + The name of the comparison grid to use for remapping (if any). - timeMean = dsWeightedSum / days.where(days > 0.) + mappingFileSection, mappingFileOption : str, optional + Section and option in ``config`` where the name of the mapping file + may be given, or where it will be stored if a new mapping file is + created - return timeMean # }}} + mappingFilePrefix : str, optional + A prefix to be prepended to the mapping file name + method : {'bilinear', 'neareststod', 'conserve'}, optional + The method of interpolation used. -def _matches_comparison(obsDescriptor, comparisonDescriptor): # {{{ - ''' - Determine if the two meshes are the same + Raises + ------ + ValueError + If comarisonGrid does not describe a known comparions grid - Authors - ------- - Xylar Asay-Davis - ''' - - if isinstance(obsDescriptor, ProjectionGridDescriptor) and \ - isinstance(comparisonDescriptor, ProjectionGridDescriptor): - # pretty hard to determine if projections are the same, so we'll rely - # on the grid names - match = obsDescriptor.meshName == comparisonDescriptor.meshName and \ - len(obsDescriptor.x) == len(comparisonDescriptor.x) and \ - len(obsDescriptor.y) == len(comparisonDescriptor.y) and \ - numpy.all(numpy.isclose(obsDescriptor.x, - comparisonDescriptor.x)) and \ - numpy.all(numpy.isclose(obsDescriptor.y, - comparisonDescriptor.y)) - elif isinstance(obsDescriptor, LatLonGridDescriptor) and \ - isinstance(comparisonDescriptor, LatLonGridDescriptor): - match = ((('degree' in obsDescriptor.units and - 'degree' in comparisonDescriptor.units) or - ('radian' in obsDescriptor.units and - 'radian' in comparisonDescriptor.units)) and - len(obsDescriptor.lat) == len(comparisonDescriptor.lat) and - len(obsDescriptor.lon) == len(comparisonDescriptor.lon) and - numpy.all(numpy.isclose(obsDescriptor.lat, - comparisonDescriptor.lat)) and - numpy.all(numpy.isclose(obsDescriptor.lon, - comparisonDescriptor.lon))) - else: - match = False - - return match # }}} + Authors + ------- + Xylar Asay-Davis + """ + super(ObservationClimatology, self).__init__(task, monthNames, + comparisonGrid) -def _setup_climatology_caching(ds, startYearClimo, endYearClimo, - yearsPerCacheFile, cachePrefix, - monthValues): # {{{ - ''' - Determine which cache files already exist, which are incomplete and which - years are present in each cache file (whether existing or to be created). + self.fieldName = fieldName - Authors - ------- - Xylar Asay-Davis - - Last Modified - ------------- - 04/08/2017 - ''' - - cacheInfo = [] - - cacheIndices = -1*numpy.ones(ds.dims['Time'], int) - monthsInDs = ds.month.values - yearsInDs = ds.year.values - - # figure out which files to load and which years go in each file - for firstYear in range(startYearClimo, endYearClimo+1, yearsPerCacheFile): - years = range(firstYear, firstYear+yearsPerCacheFile) - - yearString, fileSuffix = _get_year_string(years[0], years[-1]) - outputFileClimo = '{}_{}.nc'.format(cachePrefix, fileSuffix) - - done = False - if os.path.exists(outputFileClimo): - # already cached - dsCached = None - try: - dsCached = xr.open_dataset(outputFileClimo) - except IOError: - # assuming the cache file is corrupt, so deleting it. - message = 'Deleting cache file {}, which appears to have ' \ - 'been corrupted.'.format(outputFileClimo) - warnings.warn(message) - os.remove(outputFileClimo) - - monthsIfDone = len(monthValues)*len(years) - if ((dsCached is not None) and - (dsCached.attrs['totalMonths'] == monthsIfDone)): - # also complete, so we can move on - done = True - if dsCached is not None: - dsCached.close() + obsSection = '{}Observations'.format(task.componentName) - cacheIndex = len(cacheInfo) - for year in years: - for month in monthValues: - mask = numpy.logical_and(yearsInDs == year, - monthsInDs == month) - cacheIndices[mask] = cacheIndex + climatologyDirectory = build_config_full_path( + config=self.config, section='output', + relativePathOption='climatologySubdirectory', + relativePathSection=obsSection) - if numpy.count_nonzero(cacheIndices == cacheIndex) == 0: - continue + self.sourceDescriptor = obsGridDescriptor - cacheInfo.append((outputFileClimo, done, yearString)) + obsGridName = self.sourceDescriptor.meshName + comparisonGridName = self.comparisonDescriptor.meshName - ds = ds.copy() - ds.coords['cacheIndices'] = ('Time', cacheIndices) + self.climatologyFileName = '{}/{}_{}_{}.nc'.format( + climatologyDirectory, fieldName, obsGridName, monthNames) - return cacheInfo, cacheIndices # }}} + make_directories(climatologyDirectory) + if comparisonGrid is not None: + remappedDirectory = build_config_full_path( + config=self.config, section='output', + relativePathOption='remappedClimSubdirectory', + relativePathSection=obsSection) -def _cache_individual_climatologies(ds, cacheInfo, printProgress, - yearsPerCacheFile, monthValues, - calendar): # {{{ - ''' - Cache individual climatologies for later aggregation. - - Authors - ------- - Xylar Asay-Davis - - Last Modified - ------------- - 04/19/2017 - ''' - - for cacheIndex, info in enumerate(cacheInfo): - outputFileClimo, done, yearString = info - if done: - continue - dsYear = ds.where(ds.cacheIndices == cacheIndex, drop=True) - - if printProgress: - print ' {}'.format(yearString) - - totalDays = dsYear.daysInMonth.sum(dim='Time').values - - monthCount = dsYear.dims['Time'] - - climatology = compute_climatology(dsYear, monthValues, calendar, - maskVaries=False) - - climatology.attrs['totalDays'] = totalDays - climatology.attrs['totalMonths'] = monthCount - climatology.attrs['fingerprintClimo'] = fingerprint_generator() - - climatology.to_netcdf(outputFileClimo) - climatology.close() - - # }}} - - -def _cache_aggregated_climatology(startYearClimo, endYearClimo, cachePrefix, - printProgress, monthValues, - cacheInfo): # {{{ - ''' - Cache aggregated climatology from individual climatologies. - - Authors - ------- - Xylar Asay-Davis - - Last Modified - ------------- - 04/19/2017 - ''' - - yearString, fileSuffix = _get_year_string(startYearClimo, endYearClimo) - outputFileClimo = '{}_{}.nc'.format(cachePrefix, fileSuffix) - - done = False - if len(cacheInfo) == 0: - climatology = None - done = True - - if os.path.exists(outputFileClimo): - # already cached - climatology = None - try: - climatology = xr.open_dataset(outputFileClimo) - - except IOError: - # assuming the cache file is corrupt, so deleting it. - message = 'Deleting cache file {}, which appears to have ' \ - 'been corrupted.'.format(outputFileClimo) - warnings.warn(message) - os.remove(outputFileClimo) - - if len(cacheInfo) == 1 and outputFileClimo == cacheInfo[0][0]: - # theres only one cache file and it already has the same name - # as the aggregated file so no need to aggregate - done = True - - elif climatology is not None: - monthsIfDone = (endYearClimo-startYearClimo+1)*len(monthValues) - if climatology.attrs['totalMonths'] == monthsIfDone: - # also complete, so we can move on - done = True - else: - climatology.close() - - if not done: - if printProgress: - print ' Computing aggregated climatology ' \ - '{}...'.format(yearString) - - first = True - for cacheIndex, info in enumerate(cacheInfo): - inFileClimo = info[0] - ds = xr.open_dataset(inFileClimo) - days = ds.attrs['totalDays'] - months = ds.attrs['totalMonths'] - if first: - totalDays = days - totalMonths = months - climatology = ds * days - first = False + if self._matches_comparison(self.sourceDescriptor, + self.comparisonDescriptor): + self.remappedFileName = self.climatologyFileName else: - totalDays += days - totalMonths += months - climatology = climatology + ds * days + make_directories(remappedDirectory) - ds.close() - climatology = climatology / totalDays + self.remappedFileName = '{}/{}_{}_to_{}_{}.nc'.format( + remappedDirectory, fieldName, obsGridName, + comparisonGridName, monthNames) - climatology.attrs['totalDays'] = totalDays - climatology.attrs['totalMonths'] = totalMonths - climatology.attrs['fingerprintClimo'] = fingerprint_generator() - - climatology.to_netcdf(outputFileClimo) - - return climatology # }}} + if os.path.exists(self.remappedFileName): + # no need to create the remapper, since the cached data set + # alredy exists + self.remappedDataSet = xr.open_dataset(self.remappedFileName) + else: + self.create_remapper(mappingFileSection, mappingFileOption, + mappingFilePrefix, method) -def _get_year_string(startYear, endYear): - if startYear == endYear: - yearString = '{:04d}'.format(startYear) - fileSuffix = 'year{}'.format(yearString) - else: - yearString = '{:04d}-{:04d}'.format(startYear, endYear) - fileSuffix = 'years{}'.format(yearString) + self.remappedDataSet = None + # }}} - return yearString, fileSuffix + # }}} # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/shared/generalized_reader/__init__.py b/mpas_analysis/shared/generalized_reader/__init__.py index e69de29bb..07c0a4eda 100644 --- a/mpas_analysis/shared/generalized_reader/__init__.py +++ b/mpas_analysis/shared/generalized_reader/__init__.py @@ -0,0 +1 @@ +from generalized_reader import open_multifile_dataset \ No newline at end of file diff --git a/mpas_analysis/shared/generalized_reader/generalized_reader.py b/mpas_analysis/shared/generalized_reader/generalized_reader.py index 91ed3a9d8..e2b7f15fc 100644 --- a/mpas_analysis/shared/generalized_reader/generalized_reader.py +++ b/mpas_analysis/shared/generalized_reader/generalized_reader.py @@ -39,7 +39,7 @@ def open_multifile_dataset(fileNames, calendar, config, Parameters ---------- fileNames : list of strings - A lsit of file paths to read + A list of file paths to read calendar : {'gregorian', 'gregorian_noleap'}, optional The name of one of the calendars supported by MPAS cores diff --git a/mpas_analysis/shared/grid/grid.py b/mpas_analysis/shared/grid/grid.py index 399c6c19e..2fff25ec0 100644 --- a/mpas_analysis/shared/grid/grid.py +++ b/mpas_analysis/shared/grid/grid.py @@ -15,9 +15,6 @@ ------ Xylar Asay-Davis -Last Modified -------------- -04/16/2017 ''' import netCDF4 @@ -34,10 +31,6 @@ class MeshDescriptor(object): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/13/2017 ''' def __init__(self): # {{{ @@ -50,9 +43,6 @@ def __init__(self): # {{{ ------ Xylar Asay-Davis - Last Modified - ------------- - 04/13/2017 ''' self.meshName = None # }}} @@ -70,10 +60,6 @@ def to_scrip(self, scripFileName): # {{{ Authors ------ Xylar Asay-Davis - - Last Modified - ------------- - 03/17/2017 ''' return # }}} @@ -88,10 +74,6 @@ class MpasMeshDescriptor(MeshDescriptor): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/16/2017 ''' def __init__(self, fileName, meshName=None): # {{{ @@ -112,10 +94,6 @@ def __init__(self, fileName, meshName=None): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/16/2017 ''' ds = xarray.open_dataset(fileName) @@ -153,10 +131,6 @@ def to_scrip(self, scripFileName): # {{{ Authors ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/16/2017 ''' self.scripFileName = scripFileName @@ -222,10 +196,6 @@ class LatLonGridDescriptor(MeshDescriptor): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/16/2017 ''' def __init__(self): # {{{ ''' @@ -239,22 +209,25 @@ def __init__(self): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/05/2017 ''' self.regional = False self.meshName = None # }}} - def read(self, fileName, latVarName='lat', lonVarName='lon'): # {{{ + @classmethod + def read(cls, fileName=None, ds=None, latVarName='lat', + lonVarName='lon'): # {{{ ''' Read the lat-lon grid from a file with the given lat/lon var names. Parameters ---------- - fileName : str - The path of the file containing the lat-lon grid + fileName : str, optional + The path of the file containing the lat-lon grid (if ``ds`` is not + supplied directly) + + ds : ``xarray.Dataset`` object, optional + The path of the file containing the lat-lon grid (if supplied, + ``fileName`` will be ignored) latVarName, lonVarName : str, optional The name of the latitude and longitude variables in the grid file @@ -262,38 +235,39 @@ def read(self, fileName, latVarName='lat', lonVarName='lon'): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 03/17/2017 ''' - ds = xarray.open_dataset(fileName) + if ds is None: + ds = xarray.open_dataset(fileName) - if self.meshName is None and 'meshName' in ds.attrs: - self.meshName = ds.attrs['meshName'] + descriptor = cls() + + if descriptor.meshName is None and 'meshName' in ds.attrs: + descriptor.meshName = ds.attrs['meshName'] # Get info from input file - self.lat = numpy.array(ds[latVarName].values, float) - self.lon = numpy.array(ds[lonVarName].values, float) + descriptor.lat = numpy.array(ds[latVarName].values, float) + descriptor.lon = numpy.array(ds[lonVarName].values, float) if 'degree' in ds[latVarName].units: - self.units = 'degrees' + descriptor.units = 'degrees' else: - self.units = 'radians' + descriptor.units = 'radians' - self._set_coords(latVarName, lonVarName, ds[latVarName].dims[0], - ds[lonVarName].dims[0]) + descriptor._set_coords(latVarName, lonVarName, ds[latVarName].dims[0], + ds[lonVarName].dims[0]) # interp/extrap corners - self.lonCorner = _interp_extrap_corner(self.lon) - self.latCorner = _interp_extrap_corner(self.lat) + descriptor.lonCorner = _interp_extrap_corner(descriptor.lon) + descriptor.latCorner = _interp_extrap_corner(descriptor.lat) if 'history' in ds.attrs: - self.history = '\n'.join([ds.attrs['history'], - ' '.join(sys.argv[:])]) + descriptor.history = '\n'.join([ds.attrs['history'], + ' '.join(sys.argv[:])]) else: - self.history = sys.argv[:] # }}} + descriptor.history = sys.argv[:] + return descriptor # }}} - def create(self, latCorner, lonCorner, units='degrees'): # {{{ + @classmethod + def create(cls, latCorner, lonCorner, units='degrees'): # {{{ ''' Create the lat-lon grid with the given arrays and units. @@ -309,19 +283,17 @@ def create(self, latCorner, lonCorner, units='degrees'): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 03/17/2017 ''' + descriptor = cls() - self.latCorner = latCorner - self.lonCorner = lonCorner - self.lon = 0.5*(lonCorner[0:-1] + lonCorner[1:]) - self.lat = 0.5*(latCorner[0:-1] + latCorner[1:]) - self.units = units - self.history = sys.argv[:] - self._set_coords('lat', 'lon', 'lat', 'lon') # }}} + descriptor.latCorner = latCorner + descriptor.lonCorner = lonCorner + descriptor.lon = 0.5*(lonCorner[0:-1] + lonCorner[1:]) + descriptor.lat = 0.5*(latCorner[0:-1] + latCorner[1:]) + descriptor.units = units + descriptor.history = sys.argv[:] + descriptor._set_coords('lat', 'lon', 'lat', 'lon') + return descriptor # }}} def to_scrip(self, scripFileName): # {{{ ''' @@ -335,10 +307,6 @@ def to_scrip(self, scripFileName): # {{{ Authors ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/16/2017 ''' self.scripFileName = scripFileName @@ -407,10 +375,6 @@ class ProjectionGridDescriptor(MeshDescriptor): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/16/2017 ''' def __init__(self, projection): # {{{ @@ -426,16 +390,14 @@ def __init__(self, projection): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/05/2017 ''' self.projection = projection self.latLonProjection = pyproj.Proj(proj='latlong', datum='WGS84') self.regional = True - def read(self, fileName, meshName=None, xVarName='x', yVarName='y'): # {{{ + @classmethod + def read(cls, projection, fileName, meshName=None, xVarName='x', + yVarName='y'): # {{{ ''' Given a grid file with x and y coordinates defining the axes of the logically rectangular grid, read in the x and y coordinates and @@ -443,6 +405,10 @@ def read(self, fileName, meshName=None, xVarName='x', yVarName='y'): # {{{ Parameters ---------- + projection : pyproj.Proj object + The projection used to map from grid x-y space to latitude and + longitude + fileName : str The path of the file containing the grid data @@ -457,40 +423,39 @@ def read(self, fileName, meshName=None, xVarName='x', yVarName='y'): # {{{ Authors ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/16/2017 ''' + descriptor = cls(projection) ds = xarray.open_dataset(fileName) if meshName is None: if 'meshName' not in ds.attrs: raise ValueError('No meshName provided or found in file.') - self.meshName = ds.attrs['meshName'] + descriptor.meshName = ds.attrs['meshName'] else: - self.meshName = meshName + descriptor.meshName = meshName # Get info from input file - self.x = numpy.array(ds[xVarName].values, float) - self.y = numpy.array(ds[yVarName].values, float) + descriptor.x = numpy.array(ds[xVarName].values, float) + descriptor.y = numpy.array(ds[yVarName].values, float) - self._set_coords(xVarName, yVarName, ds[xVarName].dims[0], - ds[yVarName].dims[0]) + descriptor._set_coords(xVarName, yVarName, ds[xVarName].dims[0], + ds[yVarName].dims[0]) # interp/extrap corners - self.xCorner = _interp_extrap_corner(self.x) - self.yCorner = _interp_extrap_corner(self.y) + descriptor.xCorner = _interp_extrap_corner(descriptor.x) + descriptor.yCorner = _interp_extrap_corner(descriptor.y) # Update history attribute of netCDF file if 'history' in ds.attrs: - self.history = '\n'.join([ds.attrs['history'], - ' '.join(sys.argv[:])]) + descriptor.history = '\n'.join([ds.attrs['history'], + ' '.join(sys.argv[:])]) else: - self.history = sys.argv[:] # }}} + descriptor.history = sys.argv[:] + return descriptor # }}} - def create(self, x, y, meshName): # {{{ + @classmethod + def create(cls, projection, x, y, meshName): # {{{ ''' Given x and y coordinates defining the axes of the logically rectangular grid, save the coordinates interpolate/extrapolate to @@ -508,23 +473,21 @@ def create(self, x, y, meshName): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 03/20/2017 ''' + descriptor = cls(projection) - self.meshName = meshName + descriptor.meshName = meshName - self.x = x - self.y = y + descriptor.x = x + descriptor.y = y - self._set_coords('x', 'y', 'x', 'y') + descriptor._set_coords('x', 'y', 'x', 'y') # interp/extrap corners - self.xCorner = _interp_extrap_corner(self.x) - self.yCorner = _interp_extrap_corner(self.y) - self.history = sys.argv[:] # }}} + descriptor.xCorner = _interp_extrap_corner(descriptor.x) + descriptor.yCorner = _interp_extrap_corner(descriptor.y) + descriptor.history = sys.argv[:] + return descriptor # }}} def to_scrip(self, scripFileName): # {{{ ''' @@ -538,10 +501,6 @@ def to_scrip(self, scripFileName): # {{{ Authors ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/16/2017 ''' self.scripFileName = scripFileName @@ -594,10 +553,6 @@ def project_to_lat_lon(self, X, Y): # {{{ Authors ------ Xylar Asay-Davis - - Last Modified - ------------- - 03/20/2017 ''' Lon, Lat = pyproj.transform(self.projection, self.latLonProjection, @@ -659,10 +614,6 @@ def _create_scrip(outFile, grid_size, grid_corners, grid_rank, units, Authors ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/16/2017 ''' # Write to output file # Dimensions diff --git a/mpas_analysis/shared/interpolation/remapper.py b/mpas_analysis/shared/interpolation/remapper.py index 8a4c9c979..ae30e80fa 100644 --- a/mpas_analysis/shared/interpolation/remapper.py +++ b/mpas_analysis/shared/interpolation/remapper.py @@ -11,10 +11,6 @@ Author ------ Xylar Asay-Davis - -Last Modified -------------- -04/13/2017 ''' import subprocess @@ -29,7 +25,6 @@ from ..grid import MpasMeshDescriptor, LatLonGridDescriptor, \ ProjectionGridDescriptor - class Remapper(object): ''' A class for remapping fields using a given mapping file. The weights and @@ -38,10 +33,6 @@ class Remapper(object): Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/13/2017 ''' def __init__(self, sourceDescriptor, destinationDescriptor, @@ -75,10 +66,6 @@ def __init__(self, sourceDescriptor, destinationDescriptor, Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/13/2017 ''' if not isinstance(sourceDescriptor, @@ -128,10 +115,6 @@ def build_mapping_file(self, method='bilinear', Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/13/2017 ''' if self.mappingFileName is None or \ @@ -221,10 +204,6 @@ def remap_file(self, inFileName, outFileName, Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/13/2017 ''' if self.mappingFileName is None: @@ -315,10 +294,6 @@ def remap(self, ds, renormalizationThreshold=None): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/13/2017 ''' if self.mappingFileName is None: @@ -368,10 +343,6 @@ def _load_mapping(self): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/06/2017 ''' if self.mappingLoaded: @@ -443,10 +414,6 @@ def _remap_data_array(self, dataArray, renormalizationThreshold): # {{{ Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/05/2017 ''' sourceDims = self.sourceDescriptor.dims @@ -518,10 +485,6 @@ def _remap_numpy_array(self, inField, remapAxes, Author ------ Xylar Asay-Davis - - Last Modified - ------------- - 04/05/2017 ''' # permute the dimensions of inField so the axes to remap are first, @@ -581,16 +544,6 @@ def _remap_numpy_array(self, inField, remapAxes, return outField # }}} -def _get_lock_path(fileName): # {{{ - '''Returns the name of a temporary lock file unique to a given file name''' - directory = '{}/.locks/'.format(os.path.dirname(fileName)) - try: - os.makedirs(directory) - except OSError: - pass - return '{}/{}.lock'.format(directory, os.path.basename(fileName)) # }}} - - def _get_temp_path(): # {{{ '''Returns the name of a temporary NetCDF file''' return '{}/{}.nc'.format(tempfile._get_default_tempdir(), diff --git a/mpas_analysis/shared/io/__init__.py b/mpas_analysis/shared/io/__init__.py index 0571c93b9..41d5fefdf 100644 --- a/mpas_analysis/shared/io/__init__.py +++ b/mpas_analysis/shared/io/__init__.py @@ -1,2 +1,3 @@ from .namelist_streams_interface import NameList, StreamsFile -from .utility import paths +from .utility import paths, make_directories, build_config_full_path,\ + check_path_exists diff --git a/mpas_analysis/shared/io/utility.py b/mpas_analysis/shared/io/utility.py index c36be1d90..5cb843f9d 100644 --- a/mpas_analysis/shared/io/utility.py +++ b/mpas_analysis/shared/io/utility.py @@ -12,7 +12,7 @@ import string -def paths(*args): # {{{ +def paths(*args): # {{{ """ Returns glob'd paths in list for arbitrary number of function arguments. Note, each expanded set of paths is sorted. @@ -23,21 +23,21 @@ def paths(*args): # {{{ paths = [] for aargs in args: paths += sorted(glob.glob(aargs)) - return paths # }}} + return paths # }}} def fingerprint_generator(size=12, - chars=string.ascii_uppercase + string.digits): # {{{ + chars=string.ascii_uppercase + string.digits): # {{{ """ - Returns a random string that can be used as a unique fingerprint + Returns a random string that can be used as a unique fingerprint Reference: http://stackoverflow.com/questions/2257441/random-string-generation-with-upper-case-letters-and-digits-in-python - + Phillip J. Wolfram 04/27/2017 """ - return ''.join(random.choice(chars) for _ in range(size)) # }}} + return ''.join(random.choice(chars) for _ in range(size)) # }}} def make_directories(path): # {{{ @@ -59,7 +59,7 @@ def make_directories(path): # {{{ def build_config_full_path(config, section, relativePathOption, relativePathSection=None, - defaultPath=None): # {{{ + defaultPath=None): # {{{ """ Returns a full path from a base directory and a relative path @@ -102,7 +102,7 @@ def build_config_full_path(config, section, relativePathOption, if defaultPath is not None and not os.path.exists(fullPath): fullPath = defaultPath - return fullPath # }}} + return fullPath # }}} def check_path_exists(path): # {{{ diff --git a/mpas_analysis/shared/mpas_xarray/mpas_xarray.py b/mpas_analysis/shared/mpas_xarray/mpas_xarray.py index 53f5e7150..8c9372ce0 100644 --- a/mpas_analysis/shared/mpas_xarray/mpas_xarray.py +++ b/mpas_analysis/shared/mpas_xarray/mpas_xarray.py @@ -19,7 +19,7 @@ Last modified ------------- -02/22/2017 +05/08/2017 """ @@ -119,11 +119,11 @@ def open_multifile_dataset(fileNames, calendar, def subset_variables(ds, variableList): # {{{ """ Given a data set and a list of variable names, returns a new data set that - contains only variables with those names. + contains only variables or coords with those names. Parameters ---------- - ds : xarray.DataSet object + ds : ``xarray.DataSet`` object The data set from which a subset of variables is to be extracted. variableList : string or list of strings @@ -131,9 +131,9 @@ def subset_variables(ds, variableList): # {{{ Returns ------- - ds : xarray.DataSet object - A copy of the original data set with only the variables in - variableList. + ds : ``xarray.DataSet`` object + A copy of the original data set with only the variables and/or coords + in ``variableList``. Raises ------ @@ -146,10 +146,11 @@ def subset_variables(ds, variableList): # {{{ Last modified ------------- - 02/16/2017 + 05/08/2017 """ allvars = ds.data_vars.keys() + allcoords = ds.coords.keys() # get set of variables to drop (all ds variables not in vlist) dropvars = set(allvars) - set(variableList) @@ -161,17 +162,19 @@ def subset_variables(ds, variableList): # {{{ coords = set() for avar in ds.data_vars.keys(): coords |= set(ds[avar].coords.keys()) + coords |= set(variableList) dropcoords = set(ds.coords.keys()) - coords # drop spurious coordinates ds = ds.drop(dropcoords) - if len(ds.data_vars.keys()) == 0: + if len(ds.data_vars.keys()) == 0 and len(ds.coords.keys()) == 0: raise ValueError( 'Empty dataset is returned.\n' 'Variables {}\n' 'are not found within the dataset ' - 'variables: {}.'.format(variableList, allvars)) + 'variables: {}\n' + 'or coords: {}.'.format(variableList, allvars, allcoords)) return ds # }}} diff --git a/mpas_analysis/shared/timekeeping/utility.py b/mpas_analysis/shared/timekeeping/utility.py index 2f3c4d08c..b12aa2a88 100644 --- a/mpas_analysis/shared/timekeeping/utility.py +++ b/mpas_analysis/shared/timekeeping/utility.py @@ -4,18 +4,17 @@ Author ------ Xylar Asay-Davis - -Last Modified -------------- -02/11/2017 """ import datetime import netCDF4 import numpy +import warnings from .MpasRelativeDelta import MpasRelativeDelta +from ..constants import constants + def get_simulation_start_time(streams): """ @@ -41,10 +40,6 @@ def get_simulation_start_time(streams): Author ------ Xylar Asay-Davis - - Last modified - ------------- - 02/11/2017 """ try: @@ -99,10 +94,6 @@ def string_to_datetime(dateString): # {{{ Author ------ Xylar Asay-Davis - - Last modified - ------------- - 02/04/2017 """ (year, month, day, hour, minute, second) = \ @@ -151,10 +142,6 @@ def string_to_relative_delta(dateString, calendar='gregorian'): # {{{ Author ------ Xylar Asay-Davis - - Last modified - ------------- - 02/04/2017 """ (years, months, days, hours, minutes, seconds) = \ @@ -214,10 +201,6 @@ def string_to_days_since_date(dateString, calendar='gregorian', Author ------ Xylar Asay-Davis - - Last modified - ------------- - 02/04/2017 """ isSingleString = isinstance(dateString, str) @@ -267,10 +250,6 @@ def days_to_datetime(days, calendar='gregorian', referenceDate='0001-01-01'): Author ------ Xylar Asay-Davis - - Last modified - ------------- - 02/04/2017 """ datetimes = netCDF4.num2date(days, @@ -324,10 +303,6 @@ def datetime_to_days(dates, calendar='gregorian', referenceDate='0001-01-01'): Author ------ Xylar Asay-Davis - - Last modified - ------------- - 02/11/2017 """ isSingleDate = False @@ -377,10 +352,6 @@ def date_to_days(year=1, month=1, day=1, hour=0, minute=0, second=0, Author ------ Xylar Asay-Davis - - Last modified - ------------- - 02/11/2017 """ calendar = _mpas_to_netcdf_calendar(calendar) @@ -391,6 +362,71 @@ def date_to_days(year=1, month=1, day=1, hour=0, minute=0, second=0, calendar=calendar) +def add_years_months_days_in_month(ds, calendar): # {{{ + ''' + Add ``year``, ``month`` and ``daysInMonth`` as data arrays in ``ds``. + The number of days in each month of ``ds`` is computed either using the + ``startTime`` and ``endTime`` if available or assuming + ``gregorian_noleap`` calendar and ignoring leap years. ``year`` and + ``month`` are computed accounting correctly for the the calendar. + + Parameters + ---------- + ds : ``xarray.Dataset`` or ``xarray.DataArray`` object + A data set with a ``Time`` coordinate expressed as days since + 0001-01-01 + + calendar : {'gregorian', 'gregorian_noleap'}, optional + A calendar to be used to convert days to a `datetime.datetime` object. + + Returns + ------- + ds : object of same type as ``ds`` + The data set with ``year``, ``month`` and ``daysInMonth`` data + arrays added (if not already present) + + Authors + ------- + Xylar Asay-Davis + ''' + + if ('year' in ds.coords and 'month' in ds.coords and + 'daysInMonth' in ds.coords): + return ds + + ds = ds.copy() + + if 'year' not in ds.coords or 'month' not in ds.coords: + if calendar is None: + raise ValueError('calendar must be provided if month and year ' + 'coordinate is not in ds.') + datetimes = days_to_datetime(ds.Time, calendar=calendar) + + if 'year' not in ds.coords: + ds.coords['year'] = ('Time', [date.year for date in datetimes]) + + if 'month' not in ds.coords: + ds.coords['month'] = ('Time', [date.month for date in datetimes]) + + if 'daysInMonth' not in ds.coords: + if 'startTime' in ds.coords and 'endTime' in ds.coords: + ds.coords['daysInMonth'] = ds.endTime - ds.startTime + else: + if calendar == 'gregorian': + message = 'The MPAS run used the Gregorian calendar but ' \ + 'does not appear to have\n' \ + 'supplied start and end times. Climatologies ' \ + 'will be computed with\n' \ + 'month durations ignoring leap years.' + warnings.warn(message) + + daysInMonth = numpy.array([constants.daysInMonth[month-1] for + month in ds.month.values], float) + ds.coords['daysInMonth'] = ('Time', daysInMonth) + + return ds # }}} + + def _parse_date_string(dateString, isInterval=False): # {{{ """ Given a string containing a date, returns a tuple defining a date of the @@ -432,10 +468,6 @@ def _parse_date_string(dateString, isInterval=False): # {{{ Author ------ Xylar Asay-Davis - - Last modified - ------------- - 02/04/2017 """ if isInterval: offset = 0 diff --git a/mpas_analysis/shared/variable_namelist_stream_maps/ocean_maps.py b/mpas_analysis/shared/variable_namelist_stream_maps/ocean_maps.py index 22b7738a4..4d6487412 100644 --- a/mpas_analysis/shared/variable_namelist_stream_maps/ocean_maps.py +++ b/mpas_analysis/shared/variable_namelist_stream_maps/ocean_maps.py @@ -80,12 +80,12 @@ 'time_avg_dThreshMLD_1', 'timeMonthly_avg_dThreshMLD'] -oceanVariableMap['sst'] = \ +oceanVariableMap['temperature'] = \ ['time_avg_activeTracers_temperature', 'time_avg_activeTracers_temperature_1', 'timeMonthly_avg_activeTracers_temperature'] -oceanVariableMap['sss'] = \ +oceanVariableMap['salinity'] = \ ['time_avg_activeTracers_salinity', 'time_avg_activeTracers_salinity_1', 'timeMonthly_avg_activeTracers_salinity'] diff --git a/mpas_analysis/test/test_analysis_task b/mpas_analysis/test/test_analysis_task new file mode 120000 index 000000000..9a2506922 --- /dev/null +++ b/mpas_analysis/test/test_analysis_task @@ -0,0 +1 @@ +test_climatology \ No newline at end of file diff --git a/mpas_analysis/test/test_analysis_task.py b/mpas_analysis/test/test_analysis_task.py index 2a823329a..35ba8eda6 100644 --- a/mpas_analysis/test/test_analysis_task.py +++ b/mpas_analysis/test/test_analysis_task.py @@ -5,15 +5,64 @@ """ import pytest -from mpas_analysis.test import TestCase +import shutil +import tempfile +import os + +from mpas_analysis.test import TestCase, loaddatadir from mpas_analysis.shared.analysis_task import AnalysisTask from mpas_analysis.configuration.MpasAnalysisConfigParser \ import MpasAnalysisConfigParser +@pytest.mark.usefixtures("loaddatadir") class TestAnalysisTask(TestCase): + def setUp(self): + # Create a temporary directory + self.test_dir = tempfile.mkdtemp() + + def tearDown(self): + # Remove the directory after the test + shutil.rmtree(self.test_dir) + + def setup_config(self): + config = MpasAnalysisConfigParser() + config.read('config.default') + config.set('input', 'baseDirectory', str(self.datadir)) + config.set('input', 'mpasMeshName', 'QU240') + + config.set('output', 'baseDirectory', self.test_dir) + config.set('output', 'mappingSubdirectory', '.') + + config.set('climatology', 'startYear', '2') + config.set('climatology', 'endYear', '2') + + return config - def test_checkGenerate(self): + def setup_task(self, config): + task = AnalysisTask(config=config, + taskName='genericClimatology', + componentName='ocean', + tags=['climatology', 'horizontalMap']) + task.setup_and_check() + return task + + def test_setup_and_check(self): + config = self.setup_config() + task = self.setup_task(config) + # make sure everything was set up as expected + self.assertEqual(task.calendar, 'gregorian_noleap') + self.assertEqual(task.runDirectory, '{}/.'.format(self.datadir)) + self.assertEqual(task.historyDirectory, '{}/.'.format(self.datadir)) + self.assertEqual(task.plotsDirectory, '{}/plots'.format(self.test_dir)) + assert(os.path.exists(task.plotsDirectory)) + assert(task.namelistMap is not None) + assert(task.streamMap is not None) + assert(task.variableMap is not None) + assert(config.has_option('climatology', 'startDate')) + assert(config.has_option('climatology', 'endDate')) + + def test_check_generate(self): def doTest(generate, expectedResults): config = MpasAnalysisConfigParser() @@ -145,5 +194,73 @@ def doTest(generate, expectedResults): expectedResults['timeSeriesOHC'] = False doTest("['all', 'no_timeSeriesOHC']", expectedResults) + def test_update_start_end_date(self): + config = self.setup_config() + task = self.setup_task(config) + + inputFileNames = \ + task.get_input_file_names(streamName='timeSeriesStats', + startAndEndDateSection='climatology') + + timeCache = task.cache_multifile_dataset_times( + inputFileNames, streamName='timeSeriesStats', + timeVariableName='Time') + + # make sure the times have been cached + assert(os.path.exists( + '{}/timecache/ocean_timeSeriesStats_times.pickle'.format( + self.test_dir))) + + for index, fileName in enumerate(timeCache.keys()): + assert(fileName == os.path.abspath(inputFileNames[index])) + assert(timeCache[fileName]['years'][0] == 2) + assert(timeCache[fileName]['months'][0] == index+1) + + changed = task.update_start_end_date(section='climatology', + streamName='timeSeriesStats') + + assert(not changed) + assert(config.getint('climatology', 'startYear') == 2) + assert(config.getint('climatology', 'endYear') == 2) + assert(config.get('climatology', 'startDate') == '0002-01-01_00:00:00') + assert(config.get('climatology', 'endDate') == '0002-12-31_23:59:59') + + config.set('climatology', 'endYear', '50') + + with self.assertWarns('climatology start and/or end year different ' + 'from requested'): + changed = task.update_start_end_date(section='climatology', + streamName='timeSeriesStats') + + assert(changed) + assert(config.getint('climatology', 'startYear') == 2) + assert(config.getint('climatology', 'endYear') == 2) + assert(config.get('climatology', 'startDate') == '0002-01-01_00:00:00') + assert(config.get('climatology', 'endDate') == '0002-12-31_23:59:59') + + def test_get_input_file_names(self): + config = self.setup_config() + task = self.setup_task(config) + + inputFileNames = \ + task.get_input_file_names(streamName='timeSeriesStats', + startDate='0002-01-01_00:00:00', + endDate='0002-03-01_00:00:00') + + for index, month in enumerate([1, 2]): + expectedFileName = \ + '{}/./timeSeries.0002-{:02d}-01.nc'.format(str(self.datadir), + month) + assert(inputFileNames[index] == expectedFileName) + + inputFileNames = \ + task.get_input_file_names(streamName='timeSeriesStats', + startAndEndDateSection='climatology') + + for index, month in enumerate([1, 2, 3]): + expectedFileName = \ + '{}/./timeSeries.0002-{:02d}-01.nc'.format(str(self.datadir), + month) + assert(inputFileNames[index] == expectedFileName) # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/test/test_climatology.py b/mpas_analysis/test/test_climatology.py index 55754d567..12b2a2f04 100644 --- a/mpas_analysis/test/test_climatology.py +++ b/mpas_analysis/test/test_climatology.py @@ -11,20 +11,23 @@ import os import numpy import xarray +from functools import partial from mpas_analysis.test import TestCase, loaddatadir -from mpas_analysis.shared.generalized_reader.generalized_reader \ +from mpas_analysis.shared.generalized_reader \ import open_multifile_dataset from mpas_analysis.configuration.MpasAnalysisConfigParser \ import MpasAnalysisConfigParser -from mpas_analysis.shared.climatology import \ - get_lat_lon_comparison_descriptor, get_remapper, \ - get_mpas_climatology_file_names, get_observation_climatology_file_names, \ - add_years_months_days_in_month, compute_climatology, \ - compute_monthly_climatology, update_start_end_year, cache_climatologies +from mpas_analysis.shared.climatology import Climatology, \ + MpasClimatology, ObservationClimatology from mpas_analysis.shared.grid import MpasMeshDescriptor, LatLonGridDescriptor from mpas_analysis.shared.constants import constants +from mpas_analysis.shared.timekeeping.utility import \ + add_years_months_days_in_month + +from mpas_analysis.shared.analysis_task import AnalysisTask + @pytest.mark.usefixtures("loaddatadir") class TestClimatology(TestCase): @@ -40,100 +43,105 @@ def tearDown(self): def setup_config(self, autocloseFileLimitFraction=0.5, maxChunkSize=10000): config = MpasAnalysisConfigParser() - config.add_section('input') + config.read('config.default') + config.set('input', 'baseDirectory', str(self.datadir)) config.set('input', 'autocloseFileLimitFraction', str(autocloseFileLimitFraction)) config.set('input', 'maxChunkSize', str(maxChunkSize)) config.set('input', 'mpasMeshName', 'QU240') - config.add_section('output') config.set('output', 'baseDirectory', self.test_dir) config.set('output', 'mappingSubdirectory', '.') - config.set('output', 'mpasClimatologySubdirectory', 'clim/mpas') - config.set('output', 'mpasRegriddedClimSubdirectory', - 'clim/mpas/regrid') - config.add_section('climatology') config.set('climatology', 'startYear', '2') config.set('climatology', 'endYear', '2') - config.set('climatology', 'comparisonLatResolution', '0.5') - config.set('climatology', 'comparisonLonResolution', '0.5') - - config.set('climatology', 'overwriteMapping', 'False') - config.set('climatology', 'overwriteMpasClimatology', 'False') - config.set('climatology', 'mpasInterpolationMethod', 'bilinear') - - config.add_section('oceanObservations') - config.set('oceanObservations', 'interpolationMethod', 'bilinear') - config.set('oceanObservations', 'climatologySubdirectory', 'clim/obs') - config.set('oceanObservations', 'regriddedClimSubdirectory', - 'clim/obs/regrid') return config - def setup_mpas_remapper(self, config): - mpasMeshFileName = '{}/mpasMesh.nc'.format(self.datadir) - - comparisonDescriptor = \ - get_lat_lon_comparison_descriptor(config) + def setup_task(self, config): + task = AnalysisTask(config=config, + taskName='genericClimatology', + componentName='ocean', + tags=['climatology', 'horizontalMap']) + task.setup_and_check() + return task - mpasDescriptor = MpasMeshDescriptor( - mpasMeshFileName, meshName=config.get('input', 'mpasMeshName')) + def setup_mpas_climatology(self, config, task): - remapper = get_remapper( - config=config, sourceDescriptor=mpasDescriptor, - comparisonDescriptor=comparisonDescriptor, - mappingFileSection='climatology', - mappingFileOption='mpasMappingFile', - mappingFilePrefix='map', method=config.get( - 'climatology', 'mpasInterpolationMethod')) + fieldName = 'mld' + monthNames = 'JFM' + streamName = 'timeSeriesStats' - return remapper + mpasMeshFileName = '{}/mpasMesh.nc'.format(self.datadir) - def setup_obs_remapper(self, config, fieldName): + climatology = MpasClimatology(task=task, + fieldName=fieldName, + monthNames=monthNames, + streamName=streamName, + meshFileName=mpasMeshFileName, + comparisonGrid='latlon', + mappingFileSection='climatology', + mappingFileOption='mpasMappingFile', + mappingFilePrefix='map', + method=config.get( + 'climatology', + 'mpasInterpolationMethod')) + + return climatology + + def setup_obs_climatology(self, config, task): gridFileName = '{}/obsGrid.nc'.format(self.datadir) + fieldName = 'sst' + monthNames = 'JFM' - comparisonDescriptor = \ - get_lat_lon_comparison_descriptor(config) - - obsDescriptor = LatLonGridDescriptor() - obsDescriptor.read(fileName=gridFileName, latVarName='lat', - lonVarName='lon') - - remapper = \ - get_remapper( - config=config, sourceDescriptor=obsDescriptor, - comparisonDescriptor=comparisonDescriptor, + obsDescriptor = LatLonGridDescriptor.read(fileName=gridFileName, + latVarName='lat', + lonVarName='lon') + + climatology = \ + ObservationClimatology( + task=task, + fieldName=fieldName, + monthNames=monthNames, + obsGridDescriptor=obsDescriptor, + comparisonGrid='latlon', mappingFileSection='oceanObservations', mappingFileOption='sstClimatologyMappingFile', mappingFilePrefix='map_obs_{}'.format(fieldName), method=config.get('oceanObservations', 'interpolationMethod')) - return remapper - - def open_test_ds(self, config, calendar): - fileNames = ['{}/timeSeries.0002-{:02d}-01.nc'.format(self.datadir, - month) - for month in [1, 2, 3]] + return climatology + def open_ds_part(self, task, inputFileNames, startDate, endDate): variableMap = {'mld': ['timeMonthly_avg_tThreshMLD'], 'Time': [['xtime_startMonthly', 'xtime_endMonthly']]} + variableList = ['mld'] ds = open_multifile_dataset( - fileNames=fileNames, - calendar=calendar, - config=config, + fileNames=inputFileNames, + calendar=task.calendar, + config=task.config, timeVariableName='Time', variableList=variableList, - variableMap=variableMap) + variableMap=variableMap, + startDate=startDate, + endDate=endDate) + return ds + + def open_test_ds(self, task): + fileNames = ['{}/timeSeries.0002-{:02d}-01.nc'.format(self.datadir, + month) + for month in [1, 2, 3]] + ds = self.open_ds_part(task, fileNames, None, None) assert(len(ds.Time) == 3) return ds - def test_get_mpas_remapper(self): + def test_mpas_remapping(self): config = self.setup_config() + task = self.setup_task(config) defaultMappingFileName = '{}/map_QU240_to_0.5x0.5degree_' \ 'bilinear.nc'.format(self.test_dir) @@ -145,7 +153,9 @@ def test_get_mpas_remapper(self): if setName: config.set('climatology', 'mpasMappingFile', mappingFileName) - remapper = self.setup_mpas_remapper(config) + climatology = self.setup_mpas_climatology(config, task) + + remapper = climatology.remapper assert (os.path.abspath(mappingFileName) == os.path.abspath(remapper.mappingFileName)) @@ -156,13 +166,13 @@ def test_get_mpas_remapper(self): assert isinstance(remapper.destinationDescriptor, LatLonGridDescriptor) - def test_get_observations_remapper(self): + def test_observation_remapping(self): config = self.setup_config() - fieldName = 'sst' + task = self.setup_task(config) - defaultMappingFileName = '{}/map_obs_{}_1.0x1.0degree_to_' \ + defaultMappingFileName = '{}/map_obs_sst_1.0x1.0degree_to_' \ '0.5x0.5degree_bilinear.nc'.format( - self.test_dir, fieldName) + self.test_dir) explicitMappingFileName = '{}/mapping.nc'.format(self.test_dir) @@ -173,7 +183,9 @@ def test_get_observations_remapper(self): config.set('oceanObservations', 'sstClimatologyMappingFile', mappingFileName) - remapper = self.setup_obs_remapper(config, fieldName) + climatology = self.setup_obs_climatology(config, task) + + remapper = climatology.remapper assert (os.path.abspath(mappingFileName) == os.path.abspath(remapper.mappingFileName)) @@ -184,69 +196,61 @@ def test_get_observations_remapper(self): assert isinstance(remapper.destinationDescriptor, LatLonGridDescriptor) - def test_get_mpas_climatology_file_names(self): + def test_mpas_climatology_file_names(self): config = self.setup_config() - fieldName = 'sst' - monthNames = 'JFM' + task = self.setup_task(config) - remapper = self.setup_mpas_remapper(config) + climatology = self.setup_mpas_climatology(config, task) - (climatologyFileName, climatologyPrefix, regriddedFileName) = \ - get_mpas_climatology_file_names( - config, fieldName, monthNames, - remapper.sourceDescriptor.meshName, - remapper.destinationDescriptor.meshName) - expectedClimatologyFileName = '{}/clim/mpas/sst_QU240_JFM_' \ + expectedClimatologyFileName = '{}/clim/mpas/mld_QU240_JFM_' \ 'year0002.nc'.format(self.test_dir) - self.assertEqual(climatologyFileName, expectedClimatologyFileName) + self.assertEqual(climatology.climatologyFileName, + expectedClimatologyFileName) - expectedClimatologyPrefix = '{}/clim/mpas/sst_QU240_' \ + expectedClimatologyPrefix = '{}/clim/mpas/mld_QU240_' \ 'JFM'.format(self.test_dir) - self.assertEqual(climatologyPrefix, expectedClimatologyPrefix) + self.assertEqual(climatology.climatologyPrefix, + expectedClimatologyPrefix) - expectedRegriddedFileName = '{}/clim/mpas/regrid/sst_QU240_to_' \ - '0.5x0.5degree_JFM_' \ - 'year0002.nc'.format(self.test_dir) - self.assertEqual(regriddedFileName, expectedRegriddedFileName) + expectedRemappedFileName = '{}/clim/mpas/remapped/mld_QU240_to_' \ + '0.5x0.5degree_JFM_' \ + 'year0002.nc'.format(self.test_dir) + self.assertEqual(climatology.remappedFileName, + expectedRemappedFileName) - def test_get_observation_climatology_file_names(self): + def test_observation_climatology_file_names(self): config = self.setup_config() - fieldName = 'sst' - monthNames = 'JFM' - componentName = 'ocean' + task = self.setup_task(config) - remapper = self.setup_obs_remapper(config, fieldName) + climatology = self.setup_obs_climatology(config, task) - (climatologyFileName, regriddedFileName) = \ - get_observation_climatology_file_names( - config, fieldName, monthNames, componentName, remapper) expectedClimatologyFileName = '{}/clim/obs/sst_1.0x1.0degree_' \ 'JFM.nc'.format(self.test_dir) - self.assertEqual(climatologyFileName, expectedClimatologyFileName) + self.assertEqual(climatology.climatologyFileName, + expectedClimatologyFileName) - expectedRegriddedFileName = '{}/clim/obs/regrid/sst_1.0x1.0degree_' \ - 'to_0.5x0.5degree_' \ - 'JFM.nc'.format(self.test_dir) - self.assertEqual(regriddedFileName, expectedRegriddedFileName) + expectedRemappedFileName = '{}/clim/obs/remapped/sst_1.0x1.0degree_' \ + 'to_0.5x0.5degree_' \ + 'JFM.nc'.format(self.test_dir) + self.assertEqual(climatology.remappedFileName, + expectedRemappedFileName) - def test_compute_climatology(self): + def test_climatology_compute(self): config = self.setup_config() - calendar = 'gregorian_noleap' - ds = self.open_test_ds(config, calendar) + task = self.setup_task(config) + ds = self.open_test_ds(task) assert('month' not in ds.coords.keys()) assert('daysInMonth' not in ds.coords.keys()) # test add_months_and_days_in_month - ds = add_years_months_days_in_month(ds, calendar) + ds = add_years_months_days_in_month(ds, task.calendar) self.assertArrayEqual(ds.month.values, [1, 2, 3]) self.assertArrayEqual(numpy.round(ds.daysInMonth.values), [31, 28, 31]) - # test compute_climatology on a data set - monthNames = 'JFM' - monthValues = constants.monthDictionary[monthNames] - dsClimatology = compute_climatology(ds, monthValues, calendar) + climatology = self.setup_mpas_climatology(config, task) + dsClimatology = climatology.compute(ds, maskVaries=False) assert('Time' not in dsClimatology.dims.keys()) @@ -258,7 +262,7 @@ def test_compute_climatology(self): refClimatology.mld.values) # test compute_climatology on a data array - mldClimatology = compute_climatology(ds.mld, monthValues, calendar) + mldClimatology = climatology.compute(ds.mld, maskVaries=False) assert('Time' not in mldClimatology.dims) @@ -271,64 +275,27 @@ def test_compute_climatology(self): def test_compute_monthly_climatology(self): config = self.setup_config() - calendar = 'gregorian_noleap' - ds = self.open_test_ds(config, calendar) + task = self.setup_task(config) + ds = self.open_test_ds(task) - monthlyClimatology = compute_monthly_climatology(ds, calendar) + climatology = Climatology(task) + dsMonthly = climatology.compute_monthly(ds, maskVaries=False) - assert(len(monthlyClimatology.month) == 3) + assert(len(dsMonthly.month) == 3) - self.assertEqual(monthlyClimatology.data_vars.keys(), ['mld']) + self.assertEqual(dsMonthly.data_vars.keys(), ['mld']) climFileName = '{}/refMonthlyClim.nc'.format(self.datadir) refClimatology = xarray.open_dataset(climFileName) - self.assertArrayApproxEqual(monthlyClimatology.mld.values, + self.assertArrayApproxEqual(dsMonthly.mld.values, refClimatology.mld.values) - self.assertArrayApproxEqual(monthlyClimatology.month.values, + self.assertArrayApproxEqual(dsMonthly.month.values, refClimatology.month.values) - def test_update_start_end_year(self): - config = self.setup_config() - calendar = 'gregorian_noleap' - ds = self.open_test_ds(config, calendar) - - changed, startYear, endYear = \ - update_start_end_year(ds, config, calendar) - - assert(not changed) - assert(startYear == 2) - assert(endYear == 2) - - config.set('climatology', 'endYear', '50') - ds = self.open_test_ds(config, calendar) - - with self.assertWarns('climatology start and/or end year different ' - 'from requested'): - changed, startYear, endYear = \ - update_start_end_year(ds, config, calendar) - - assert(changed) - assert(startYear == 2) - assert(endYear == 2) - - def cache_climatologies_setup(self): - config = self.setup_config() - calendar = 'gregorian_noleap' - ds = self.open_test_ds(config, calendar) - fieldName = 'mld' - climFileName = '{}/refSeasonalClim.nc'.format(self.datadir) - refClimatology = xarray.open_dataset(climFileName) - - remapper = self.setup_mpas_remapper(config) - - return {'config': config, 'calendar': calendar, 'ds': ds, - 'fieldName': fieldName, 'climFileName': climFileName, - 'refClimatology': refClimatology, 'remapper': remapper} - def test_jan_1yr_climo_test1(self): - setup = self.cache_climatologies_setup() + task, refClimatology = self.cache_climatologies_setup() # test1: Just January, 1-year climatologies are cached; only one file # is produced with suffix year0002; a second run of # cache_climatologies doesn't modify any files @@ -341,10 +308,10 @@ def test_jan_1yr_climo_test1(self): 'expectedDays': 30.958333, 'expectedMonths': 1, 'refClimatology': None} - self.cache_climatologies_driver(test1, **setup) + self.cache_climatologies_driver(test1, task) def test_jfm_1yr_climo_test2(self): - setup = self.cache_climatologies_setup() + task, refClimatology = self.cache_climatologies_setup() # same as test1 but with JFM test2 = {'monthNames': 'JFM', 'monthValues': constants.monthDictionary['JFM'], @@ -354,11 +321,11 @@ def test_jfm_1yr_climo_test2(self): # weird value because first time step of Jan. missing in ds 'expectedDays': 89.958333, 'expectedMonths': 3, - 'refClimatology': setup['refClimatology']} - self.cache_climatologies_driver(test2, **setup) + 'refClimatology': refClimatology} + self.cache_climatologies_driver(test2, task) def test_jan_2yr_climo_test3(self): - setup = self.cache_climatologies_setup() + task, refClimatology = self.cache_climatologies_setup() # test3: 2-year climatologies are cached; 2 files are produced # with suffix years0002-0003 (the "individual" climatology # file) and year0002 (the "aggregated" climatology file); @@ -376,10 +343,10 @@ def test_jan_2yr_climo_test3(self): 'expectedDays': 30.958333, 'expectedMonths': 1, 'refClimatology': None} - self.cache_climatologies_driver(test3, **setup) + self.cache_climatologies_driver(test3, task) def test_jfm_2yr_climo_test4(self): - setup = self.cache_climatologies_setup() + task, refClimatology = self.cache_climatologies_setup() # test4: same as test3 but with JFM test4 = {'monthNames': 'JFM', 'monthValues': constants.monthDictionary['JFM'], @@ -389,13 +356,18 @@ def test_jfm_2yr_climo_test4(self): # weird value because first time step of Jan. missing in ds 'expectedDays': 89.958333, 'expectedMonths': 3, - 'refClimatology': setup['refClimatology']} - self.cache_climatologies_driver(test4, **setup) + 'refClimatology': refClimatology} + self.cache_climatologies_driver(test4, task) + + def cache_climatologies_setup(self): + config = self.setup_config() + task = self.setup_task(config) + climFileName = '{}/refSeasonalClim.nc'.format(self.datadir) + refClimatology = xarray.open_dataset(climFileName) + return task, refClimatology - def cache_climatologies_driver(self, test, config, fieldName, - ds, remapper, calendar, **kwargs): + def cache_climatologies_driver(self, test, task): monthNames = test['monthNames'] - monthValues = test['monthValues'] yearsPerCacheFile = test['yearsPerCacheFile'] expectedSuffixes = test['expectedSuffixes'] expectedModified = test['expectedModified'] @@ -403,17 +375,31 @@ def cache_climatologies_driver(self, test, config, fieldName, expectedMonths = test['expectedMonths'] refClimatology = test['refClimatology'] - (climatologyFileName, climatologyPrefix) = \ - get_mpas_climatology_file_names( - config, fieldName, monthNames, - remapper.sourceDescriptor.meshName) + task.config.set('climatology', 'yearsPerCacheFile', + str(yearsPerCacheFile)) + + fieldName = 'mld' + streamName = 'timeSeriesStats' + + mpasMeshFileName = '{}/mpasMesh.nc'.format(self.datadir) + + climatology = MpasClimatology(task=task, + fieldName=fieldName, + monthNames=monthNames, + streamName=streamName, + meshFileName=mpasMeshFileName, + comparisonGrid='latlon', + mappingFileSection='climatology', + mappingFileOption='mpasMappingFile', + mappingFilePrefix='map', + method=task.config.get( + 'climatology', + 'mpasInterpolationMethod')) - config.set('climatology', 'yearsPerCacheFile', - str(yearsPerCacheFile)) # once without cache files - dsClimatology = cache_climatologies( - ds, monthValues, config, climatologyPrefix, calendar, - printProgress=True) + openDataSetFunc = partial(self.open_ds_part, task) + dsClimatology = climatology.cache(openDataSetFunc=openDataSetFunc, + printProgress=True) if refClimatology is not None: self.assertArrayApproxEqual(dsClimatology.mld.values, @@ -437,9 +423,8 @@ def cache_climatologies_driver(self, test, config, fieldName, fingerprints.append(dsClimatology.fingerprintClimo) # try it again with cache files saved - dsClimatology = cache_climatologies( - ds, monthValues, config, climatologyPrefix, calendar, - printProgress=True) + dsClimatology = climatology.cache(openDataSetFunc=openDataSetFunc, + printProgress=True) if refClimatology is not None: self.assertArrayApproxEqual(dsClimatology.mld.values, @@ -468,5 +453,4 @@ def cache_climatologies_driver(self, test, config, fieldName, # remove the cache file for the next try os.remove(expectedClimatologyFileName) - # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/test/test_climatology/mpas-o_in b/mpas_analysis/test/test_climatology/mpas-o_in new file mode 100644 index 000000000..8a779b161 --- /dev/null +++ b/mpas_analysis/test/test_climatology/mpas-o_in @@ -0,0 +1,1083 @@ +&run_modes + config_ocean_run_mode = 'forward' +/ +&time_management + config_calendar_type = 'gregorian_noleap' + config_do_restart = .true. + config_restart_timestamp_name = 'rpointer.ocn' + config_start_time = 'file' +/ +&io + config_pio_num_iotasks = 0 + config_pio_stride = 1 + config_write_output_on_startup = .true. +/ +&decomposition + config_block_decomp_file_prefix = '/lustre/scratch3/turquoise/jonbob/ACME/input_data/ocn/mpas-o/oQU240/mpas-o.graph.info.151209.part.' + config_explicit_proc_decomp = .false. + config_num_halos = 3 + config_number_of_blocks = 0 + config_proc_decomp_file_prefix = 'graph.info.part.' +/ +&init_setup + config_expand_sphere = .false. + config_init_configuration = 'none' + config_realistic_coriolis_parameter = .false. + config_vert_levels = -1 + config_vertical_grid = 'uniform' + config_write_cull_cell_mask = .true. +/ +&cvtgenerator + config_1dcvtgenerator_dzseed = 1.2 + config_1dcvtgenerator_stretch1 = 1.0770 + config_1dcvtgenerator_stretch2 = 1.0275 +/ +&init_ssh_and_landicepressure + config_iterative_init_variable = 'landIcePressure' +/ +&time_integration + config_dt = '01:00:00' + config_time_integrator = 'split_explicit' +/ +&ale_vertical_grid + config_dzdk_positive = .false. + config_max_thickness_factor = 6.0 + config_min_thickness = 1.0 + config_use_min_max_thickness = .false. + config_vert_coord_movement = 'uniform_stretching' +/ +&ale_frequency_filtered_thickness + config_highfreqthick_del2 = 100.0 + config_highfreqthick_restore_time = 30.0 + config_thickness_filter_timescale = 5.0 + config_use_freq_filtered_thickness = .false. + config_use_highfreqthick_del2 = .false. + config_use_highfreqthick_restore = .false. +/ +&partial_bottom_cells + config_alter_ics_for_pbcs = .false. + config_min_pbc_fraction = 0.10 + config_pbc_alteration_type = 'full_cell' +/ +&hmix + config_apvm_scale_factor = 0.0 + config_hmix_scalewithmesh = .false. + config_maxmeshdensity = -1.0 +/ +&hmix_del2 + config_mom_del2 = 10.0 + config_tracer_del2 = 10.0 + config_use_mom_del2 = .false. + config_use_tracer_del2 = .false. +/ +&hmix_del4 + config_mom_del4 = 2.0e14 + config_mom_del4_div_factor = 1.0 + config_tracer_del4 = 0.0 + config_use_mom_del4 = .true. + config_use_tracer_del4 = .false. +/ +&hmix_leith + config_leith_dx = 15000.0 + config_leith_parameter = 1.0 + config_leith_visc2_max = 2.5e3 + config_use_leith_del2 = .false. +/ +&mesoscale_eddy_parameterization + config_gravwavespeed_trunc = 0.3 + config_max_relative_slope = 0.01 + config_redi_bottom_layer_tapering_depth = 0.0 + config_redi_kappa = 0.0 + config_redi_surface_layer_tapering_extent = 0.0 + config_standardgm_tracer_kappa = 600.0 + config_use_redi_bottom_layer_tapering = .false. + config_use_redi_surface_layer_tapering = .false. + config_use_standardgm = .true. +/ +&hmix_del2_tensor + config_mom_del2_tensor = 10.0 + config_use_mom_del2_tensor = .false. +/ +&hmix_del4_tensor + config_mom_del4_tensor = 5.0e13 + config_use_mom_del4_tensor = .false. +/ +&rayleigh_damping + config_rayleigh_damping_coeff = 0.0 + config_rayleigh_friction = .false. +/ +&vmix + config_convective_diff = 1.0 + config_convective_visc = 1.0 +/ +&vmix_const + config_use_const_diff = .false. + config_use_const_visc = .false. + config_vert_diff = 1.0e-5 + config_vert_visc = 1.0e-4 +/ +&vmix_rich + config_bkrd_vert_diff = 1.0e-5 + config_bkrd_vert_visc = 1.0e-4 + config_rich_mix = 0.005 + config_use_rich_diff = .false. + config_use_rich_visc = .false. +/ +&vmix_tanh + config_max_diff_tanh = 2.5e-2 + config_max_visc_tanh = 2.5e-1 + config_min_diff_tanh = 1.0e-5 + config_min_visc_tanh = 1.0e-4 + config_use_tanh_diff = .false. + config_use_tanh_visc = .false. + config_zmid_tanh = -100 + config_zwidth_tanh = 100 +/ +&cvmix + config_cvmix_background_diffusion = 1.0e-5 + config_cvmix_background_viscosity = 1.0e-4 + config_cvmix_convective_basedonbvf = .true. + config_cvmix_convective_diffusion = 1.0 + config_cvmix_convective_triggerbvf = 0.0 + config_cvmix_convective_viscosity = 1.0 + config_cvmix_kpp_boundary_layer_depth = 30.0 + config_cvmix_kpp_criticalbulkrichardsonnumber = 0.25 + config_cvmix_kpp_ekmanobl = .false. + config_cvmix_kpp_interpolationomltype = 'quadratic' + config_cvmix_kpp_matching = 'SimpleShapes' + config_cvmix_kpp_monobobl = .false. + config_cvmix_kpp_stop_obl_search = 100.0 + config_cvmix_kpp_surface_layer_averaging = 5.0 + config_cvmix_kpp_surface_layer_extent = 0.1 + config_cvmix_kpp_use_enhanced_diff = .true. + config_cvmix_prandtl_number = 1.0 + config_cvmix_shear_kpp_exp = 3 + config_cvmix_shear_kpp_nu_zero = 0.005 + config_cvmix_shear_kpp_ri_zero = 0.7 + config_cvmix_shear_mixing_scheme = 'KPP' + config_cvmix_shear_pp_alpha = 5.0 + config_cvmix_shear_pp_exp = 2.0 + config_cvmix_shear_pp_nu_zero = 0.005 + config_use_cvmix = .true. + config_use_cvmix_background = .true. + config_use_cvmix_convection = .true. + config_use_cvmix_double_diffusion = .false. + config_use_cvmix_fixed_boundary_layer = .false. + config_use_cvmix_kpp = .true. + config_use_cvmix_shear = .true. + config_use_cvmix_tidal_mixing = .false. + configure_cvmix_kpp_minimum_obl_under_sea_ice = 10.0 +/ +&forcing + config_flux_attenuation_coefficient = 0.001 + config_flux_attenuation_coefficient_runoff = 10.0 + config_use_bulk_thickness_flux = .true. + config_use_bulk_wind_stress = .true. +/ +&shortwaveradiation + config_forcing_restart_file = 'Restart_forcing_time_stamp' + config_jerlov_water_type = 3 + config_surface_buoyancy_depth = 1 + config_sw_absorption_type = 'jerlov' +/ +&frazil_ice + config_frazil_fractional_thickness_limit = 0.1 + config_frazil_heat_of_fusion = 3.337e5 + config_frazil_ice_density = 1000.0 + config_frazil_in_open_ocean = .true. + config_frazil_land_ice_reference_salinity = 0.0 + config_frazil_maximum_depth = 100.0 + config_frazil_maximum_freezing_temperature = 0.0 + config_frazil_sea_ice_reference_salinity = 4.0 + config_frazil_under_land_ice = .true. + config_frazil_use_surface_pressure = .false. + config_specific_heat_sea_water = 3.996e3 + config_use_frazil_ice_formation = .true. +/ +&land_ice_fluxes + config_land_ice_flux_attenuation_coefficient = 10.0 + config_land_ice_flux_boundarylayerneighborweight = 0.0 + config_land_ice_flux_boundarylayerthickness = 10.0 + config_land_ice_flux_cp_ice = 2.009e3 + config_land_ice_flux_formulation = 'Jenkins' + config_land_ice_flux_isomip_gammat = 1e-4 + config_land_ice_flux_jenkins_heat_transfer_coefficient = 0.011 + config_land_ice_flux_jenkins_salt_transfer_coefficient = 3.1e-4 + config_land_ice_flux_mode = 'off' + config_land_ice_flux_rho_ice = 918 + config_land_ice_flux_rms_tidal_velocity = 5e-2 + config_land_ice_flux_topdragcoeff = 2.5e-3 + config_land_ice_flux_usehollandjenkinsadvdiff = .false. +/ +&advection + config_coef_3rd_order = 0.25 + config_horiz_tracer_adv_order = 3 + config_monotonic = .true. + config_vert_tracer_adv = 'stencil' + config_vert_tracer_adv_order = 3 +/ +&bottom_drag + config_bottom_drag_coeff = 1.0e-3 +/ +&ocean_constants + config_density0 = 1026.0 +/ +&pressure_gradient + config_common_level_weight = 0.5 + config_pressure_gradient_type = 'Jacobian_from_TS' +/ +&eos + config_eos_type = 'jm' + config_land_ice_cavity_freezing_temperature_coeff_0 = 6.22e-2 + config_land_ice_cavity_freezing_temperature_coeff_p = -7.43e-8 + config_land_ice_cavity_freezing_temperature_coeff_ps = -1.74e-10 + config_land_ice_cavity_freezing_temperature_coeff_s = -5.63e-2 + config_land_ice_cavity_freezing_temperature_reference_pressure = 0.0 + config_open_ocean_freezing_temperature_coeff_0 = -1.8 + config_open_ocean_freezing_temperature_coeff_p = 0.0 + config_open_ocean_freezing_temperature_coeff_ps = 0.0 + config_open_ocean_freezing_temperature_coeff_s = 0.0 + config_open_ocean_freezing_temperature_reference_pressure = 0.0 +/ +&eos_linear + config_eos_linear_alpha = 0.2 + config_eos_linear_beta = 0.8 + config_eos_linear_densityref = 1000.0 + config_eos_linear_sref = 35.0 + config_eos_linear_tref = 5.0 +/ +&split_explicit_ts + config_btr_dt = '0000_00:03:00' + config_btr_gam1_velwt1 = 0.5 + config_btr_gam2_sshwt1 = 1.0 + config_btr_gam3_velwt2 = 1.0 + config_btr_solve_ssh2 = .false. + config_btr_subcycle_loop_factor = 2 + config_n_bcl_iter_beg = 1 + config_n_bcl_iter_end = 2 + config_n_bcl_iter_mid = 2 + config_n_btr_cor_iter = 2 + config_n_ts_iter = 2 + config_vel_correction = .true. +/ +&testing + config_conduct_tests = .false. + config_tensor_test_function = 'sph_uCosCos' + config_test_tensors = .false. +/ +&debug + config_check_ssh_consistency = .true. + config_check_tracer_monotonicity = .false. + config_check_zlevel_consistency = .false. + config_disable_redi_horizontal_term1 = .false. + config_disable_redi_horizontal_term2 = .false. + config_disable_redi_horizontal_term3 = .false. + config_disable_redi_k33 = .false. + config_disable_thick_all_tend = .false. + config_disable_thick_hadv = .false. + config_disable_thick_sflux = .false. + config_disable_thick_vadv = .false. + config_disable_tr_adv = .false. + config_disable_tr_all_tend = .false. + config_disable_tr_hmix = .false. + config_disable_tr_nonlocalflux = .false. + config_disable_tr_sflux = .false. + config_disable_tr_vmix = .false. + config_disable_vel_all_tend = .false. + config_disable_vel_coriolis = .false. + config_disable_vel_hmix = .false. + config_disable_vel_pgrad = .false. + config_disable_vel_surface_stress = .false. + config_disable_vel_vadv = .false. + config_disable_vel_vmix = .false. + config_filter_btr_mode = .false. + config_include_ke_vertex = .false. + config_prescribe_thickness = .false. + config_prescribe_velocity = .false. + config_read_nearest_restart = .false. +/ +&constrain_haney_number + config_rx1_horiz_smooth_open_ocean_cells = 20 + config_rx1_horiz_smooth_weight = 1.0 + config_rx1_init_inner_weight = 0.1 + config_rx1_inner_iter_count = 10 + config_rx1_max = 5.0 + config_rx1_min_layer_thickness = 1.0 + config_rx1_min_levels = 3 + config_rx1_outer_iter_count = 20 + config_rx1_slope_weight = 1e-1 + config_rx1_vert_smooth_weight = 1.0 + config_rx1_zstar_weight = 1.0 + config_use_rx1_constraint = .false. +/ +&baroclinic_channel + config_baroclinic_channel_bottom_depth = 1000.0 + config_baroclinic_channel_bottom_temperature = 10.1 + config_baroclinic_channel_coriolis_parameter = -1.2e-4 + config_baroclinic_channel_gradient_width_dist = 40e3 + config_baroclinic_channel_gradient_width_frac = 0.08 + config_baroclinic_channel_salinity = 35.0 + config_baroclinic_channel_surface_temperature = 13.1 + config_baroclinic_channel_temperature_difference = 1.2 + config_baroclinic_channel_use_distances = .false. + config_baroclinic_channel_vert_levels = 20 +/ +&lock_exchange + config_lock_exchange_bottom_depth = 20.0 + config_lock_exchange_cold_temperature = 5.0 + config_lock_exchange_direction = 'y' + config_lock_exchange_isopycnal_min_thickness = 0.01 + config_lock_exchange_layer_type = 'z-level' + config_lock_exchange_salinity = 35.0 + config_lock_exchange_vert_levels = 20 + config_lock_exchange_warm_temperature = 30.0 +/ +&internal_waves + config_internal_waves_amplitude_width_dist = 50e3 + config_internal_waves_amplitude_width_frac = 0.33 + config_internal_waves_bottom_depth = 500.0 + config_internal_waves_bottom_temperature = 10.1 + config_internal_waves_isopycnal_displacement = 125.0 + config_internal_waves_layer_type = 'z-level' + config_internal_waves_salinity = 35.0 + config_internal_waves_surface_temperature = 20.1 + config_internal_waves_temperature_difference = 2.0 + config_internal_waves_use_distances = false + config_internal_waves_vert_levels = 20 +/ +&overflow + config_overflow_bottom_depth = 2000.0 + config_overflow_domain_temperature = 20.0 + config_overflow_isopycnal_min_thickness = 0.01 + config_overflow_layer_type = 'z-level' + config_overflow_plug_temperature = 10.0 + config_overflow_plug_width_dist = 20e3 + config_overflow_plug_width_frac = 0.10 + config_overflow_ridge_depth = 500.0 + config_overflow_salinity = 35.0 + config_overflow_slope_center_dist = 40e3 + config_overflow_slope_center_frac = 0.20 + config_overflow_slope_width_dist = 7e3 + config_overflow_slope_width_frac = 0.05 + config_overflow_use_distances = false + config_overflow_vert_levels = 100 +/ +&global_ocean + config_global_ocean_chlorophyll_varname = 'none' + config_global_ocean_clearsky_varname = 'none' + config_global_ocean_cull_inland_seas = .true. + config_global_ocean_deepen_critical_passages = .true. + config_global_ocean_depress_by_land_ice = .false. + config_global_ocean_depth_conversion_factor = 1.0 + config_global_ocean_depth_dimname = 'none' + config_global_ocean_depth_file = 'none' + config_global_ocean_depth_varname = 'none' + config_global_ocean_ecosys_depth_conversion_factor = 1.0 + config_global_ocean_ecosys_depth_varname = 'none' + config_global_ocean_ecosys_file = 'unknown' + config_global_ocean_ecosys_forcing_file = 'unknown' + config_global_ocean_ecosys_forcing_time_dimname = 'none' + config_global_ocean_ecosys_lat_varname = 'none' + config_global_ocean_ecosys_latlon_degrees = .true. + config_global_ocean_ecosys_lon_varname = 'none' + config_global_ocean_ecosys_method = 'bilinear_interpolation' + config_global_ocean_ecosys_ndepth_dimname = 'none' + config_global_ocean_ecosys_nlat_dimname = 'none' + config_global_ocean_ecosys_nlon_dimname = 'none' + config_global_ocean_ecosys_vert_levels = -1 + config_global_ocean_interior_restore_rate = 1.0e-7 + config_global_ocean_land_ice_topo_draft_varname = 'none' + config_global_ocean_land_ice_topo_file = 'none' + config_global_ocean_land_ice_topo_grounded_frac_varname = 'none' + config_global_ocean_land_ice_topo_ice_frac_varname = 'none' + config_global_ocean_land_ice_topo_lat_varname = 'none' + config_global_ocean_land_ice_topo_latlon_degrees = .true. + config_global_ocean_land_ice_topo_lon_varname = 'none' + config_global_ocean_land_ice_topo_nlat_dimname = 'none' + config_global_ocean_land_ice_topo_nlon_dimname = 'none' + config_global_ocean_land_ice_topo_thickness_varname = 'none' + config_global_ocean_minimum_depth = 15 + config_global_ocean_piston_velocity = 5.0e-5 + config_global_ocean_salinity_file = 'none' + config_global_ocean_salinity_varname = 'none' + config_global_ocean_smooth_ecosys_iterations = 0 + config_global_ocean_smooth_topography = .true. + config_global_ocean_smooth_ts_iterations = 0 + config_global_ocean_swdata_file = 'none' + config_global_ocean_swdata_lat_varname = 'none' + config_global_ocean_swdata_latlon_degrees = .true. + config_global_ocean_swdata_lon_varname = 'none' + config_global_ocean_swdata_method = 'bilinear_interpolation' + config_global_ocean_swdata_nlat_dimname = 'none' + config_global_ocean_swdata_nlon_dimname = 'none' + config_global_ocean_temperature_file = 'none' + config_global_ocean_temperature_varname = 'none' + config_global_ocean_topography_file = 'none' + config_global_ocean_topography_has_ocean_frac = .false. + config_global_ocean_topography_lat_varname = 'none' + config_global_ocean_topography_latlon_degrees = .true. + config_global_ocean_topography_lon_varname = 'none' + config_global_ocean_topography_method = 'bilinear_interpolation' + config_global_ocean_topography_nlat_dimname = 'none' + config_global_ocean_topography_nlon_dimname = 'none' + config_global_ocean_topography_ocean_frac_varname = 'none' + config_global_ocean_topography_varname = 'none' + config_global_ocean_tracer_depth_conversion_factor = 1.0 + config_global_ocean_tracer_depth_varname = 'none' + config_global_ocean_tracer_lat_varname = 'none' + config_global_ocean_tracer_latlon_degrees = .true. + config_global_ocean_tracer_lon_varname = 'none' + config_global_ocean_tracer_method = 'bilinear_interpolation' + config_global_ocean_tracer_ndepth_dimname = 'none' + config_global_ocean_tracer_nlat_dimname = 'none' + config_global_ocean_tracer_nlon_dimname = 'none' + config_global_ocean_tracer_vert_levels = -1 + config_global_ocean_windstress_conversion_factor = 1 + config_global_ocean_windstress_file = 'none' + config_global_ocean_windstress_lat_varname = 'none' + config_global_ocean_windstress_latlon_degrees = .true. + config_global_ocean_windstress_lon_varname = 'none' + config_global_ocean_windstress_meridional_varname = 'none' + config_global_ocean_windstress_method = 'bilinear_interpolation' + config_global_ocean_windstress_nlat_dimname = 'none' + config_global_ocean_windstress_nlon_dimname = 'none' + config_global_ocean_windstress_zonal_varname = 'none' + config_global_ocean_zenithangle_varname = 'none' +/ +&cvmix_wswsbf + config_cvmix_wswsbf_bottom_depth = 400.0 + config_cvmix_wswsbf_coriolis_parameter = 1.0e-4 + config_cvmix_wswsbf_evaporation_flux = 0.0 + config_cvmix_wswsbf_interior_salinity_restoring_rate = 1.0e-6 + config_cvmix_wswsbf_interior_temperature_restoring_rate = 1.0e-6 + config_cvmix_wswsbf_latent_heat_flux = 0.0 + config_cvmix_wswsbf_max_windstress = 0.10 + config_cvmix_wswsbf_mixed_layer_depth_salinity = 0.0 + config_cvmix_wswsbf_mixed_layer_depth_temperature = 0.0 + config_cvmix_wswsbf_mixed_layer_salinity_change = 0.0 + config_cvmix_wswsbf_mixed_layer_temperature_change = 0.0 + config_cvmix_wswsbf_rain_flux = 0.0 + config_cvmix_wswsbf_salinity_gradient = 0.0 + config_cvmix_wswsbf_salinity_gradient_mixed_layer = 0.0 + config_cvmix_wswsbf_salinity_piston_velocity = 4.0e-6 + config_cvmix_wswsbf_sensible_heat_flux = 0.0 + config_cvmix_wswsbf_shortwave_heat_flux = 0.0 + config_cvmix_wswsbf_surface_restoring_salinity = 35.0 + config_cvmix_wswsbf_surface_restoring_temperature = 15.0 + config_cvmix_wswsbf_surface_salinity = 35.0 + config_cvmix_wswsbf_surface_temperature = 15.0 + config_cvmix_wswsbf_temperature_gradient = 0.01 + config_cvmix_wswsbf_temperature_gradient_mixed_layer = 0.0 + config_cvmix_wswsbf_temperature_piston_velocity = 4.0e-6 + config_cvmix_wswsbf_vert_levels = 100 + config_cvmix_wswsbf_vertical_grid = 'uniform' +/ +&iso + config_iso_acc_wind = 0.2 + config_iso_asf_wind = -0.05 + config_iso_cont_slope_flag = .true. + config_iso_depression_center_lon = 60 + config_iso_depression_depth = 800 + config_iso_depression_flag = .true. + config_iso_depression_north_lat = -65 + config_iso_depression_south_lat = -72 + config_iso_depression_width = 480000 + config_iso_embayment_center_lat = -71 + config_iso_embayment_center_lon = 60 + config_iso_embayment_depth = 2000 + config_iso_embayment_flag = .true. + config_iso_embayment_radius = 500000 + config_iso_heat_flux_lat_mn = -53 + config_iso_heat_flux_lat_sm = -65 + config_iso_heat_flux_lat_ss = -70 + config_iso_heat_flux_middle = 10 + config_iso_heat_flux_north = -5 + config_iso_heat_flux_region1 = -5 + config_iso_heat_flux_region1_flag = false + config_iso_heat_flux_region1_radius = 300000 + config_iso_heat_flux_region2 = -5 + config_iso_heat_flux_region2_flag = false + config_iso_heat_flux_region2_radius = 240000 + config_iso_heat_flux_south = -5 + config_iso_initial_temp_h0 = 1200 + config_iso_initial_temp_h1 = 500 + config_iso_initial_temp_latn = -50 + config_iso_initial_temp_lats = -75 + config_iso_initial_temp_mt = 0.000075 + config_iso_initial_temp_t1 = 3.5 + config_iso_initial_temp_t2 = 4.0 + config_iso_main_channel_depth = 4000.0 + config_iso_max_cont_slope = 0.01 + config_iso_north_wall_lat = -50 + config_iso_plateau_center_lat = -58 + config_iso_plateau_center_lon = 300 + config_iso_plateau_flag = .true. + config_iso_plateau_height = 2000 + config_iso_plateau_radius = 200000 + config_iso_plateau_slope_width = 1000000 + config_iso_region1_center_lat = -75 + config_iso_region1_center_lon = 60 + config_iso_region2_center_lat = -71 + config_iso_region2_center_lon = 150 + config_iso_region3_center_lat = -71 + config_iso_region3_center_lon = 240 + config_iso_region4_center_lat = -71 + config_iso_region4_center_lon = 330 + config_iso_ridge_center_lon = 180 + config_iso_ridge_flag = .true. + config_iso_ridge_height = 2000.0 + config_iso_ridge_width = 2000000 + config_iso_salinity = 35.0 + config_iso_shelf_depth = 500 + config_iso_shelf_flag = .true. + config_iso_shelf_width = 120000 + config_iso_south_wall_lat = -70 + config_iso_surface_temperature_piston_velocity = 5.787e-5 + config_iso_temperature_restore_lcx1 = 600000 + config_iso_temperature_restore_lcx2 = 600000 + config_iso_temperature_restore_lcx3 = 600000 + config_iso_temperature_restore_lcx4 = 600000 + config_iso_temperature_restore_lcy1 = 600000 + config_iso_temperature_restore_lcy2 = 250000 + config_iso_temperature_restore_lcy3 = 250000 + config_iso_temperature_restore_lcy4 = 250000 + config_iso_temperature_restore_region1_flag = .true. + config_iso_temperature_restore_region2_flag = .true. + config_iso_temperature_restore_region3_flag = .true. + config_iso_temperature_restore_region4_flag = .true. + config_iso_temperature_restore_t1 = -1 + config_iso_temperature_restore_t2 = -1 + config_iso_temperature_restore_t3 = -1 + config_iso_temperature_restore_t4 = -1 + config_iso_temperature_sponge_h1 = 1000 + config_iso_temperature_sponge_l1 = 120000 + config_iso_temperature_sponge_t1 = 10 + config_iso_temperature_sponge_tau1 = 10.0 + config_iso_vert_levels = 100 + config_iso_wind_stress_max = 0.01 + config_iso_wind_trans = -65 +/ +&soma + config_soma_bottom_depth = 2500.0 + config_soma_center_latitude = 35.0 + config_soma_center_longitude = 0.0 + config_soma_density_difference = 4.0 + config_soma_density_difference_linear = 0.05 + config_soma_domain_width = 1.25e6 + config_soma_phi = 0.1 + config_soma_ref_density = 1000.0 + config_soma_shelf_depth = 100.0 + config_soma_shelf_width = -0.4 + config_soma_surface_salinity = 33.0 + config_soma_surface_temperature = 20.0 + config_soma_thermocline_depth = 300.0 + config_soma_vert_levels = 100 +/ +&ziso + config_ziso_add_easterly_wind_stress_asf = false + config_ziso_antarctic_shelf_front_width = 600000 + config_ziso_bottom_depth = 2500.0 + config_ziso_coriolis_gradient = 1e-11 + config_ziso_frazil_enable = false + config_ziso_frazil_temperature_anomaly = -3.0 + config_ziso_initial_temp_h1 = 300.0 + config_ziso_initial_temp_mt = 7.5e-5 + config_ziso_initial_temp_t1 = 6.0 + config_ziso_initial_temp_t2 = 3.6 + config_ziso_mean_restoring_temp = 3.0 + config_ziso_meridional_extent = 2.0e6 + config_ziso_reference_coriolis = -1e-4 + config_ziso_restoring_sponge_l = 8.0e4 + config_ziso_restoring_temp_dev_ta = 2.0 + config_ziso_restoring_temp_dev_tb = 2.0 + config_ziso_restoring_temp_piston_vel = 1.93e-5 + config_ziso_restoring_temp_tau = 30.0 + config_ziso_restoring_temp_ze = 1250.0 + config_ziso_shelf_depth = 500.0 + config_ziso_slope_center_position = 5.0e5 + config_ziso_slope_half_width = 1.0e5 + config_ziso_use_slopping_bathymetry = false + config_ziso_vert_levels = 100 + config_ziso_wind_stress_max = 0.2 + config_ziso_wind_stress_shelf_front_max = -0.05 + config_ziso_wind_transition_position = 800000.0 + config_ziso_zonal_extent = 1.0e6 +/ +&sub_ice_shelf_2d + config_sub_ice_shelf_2d_bottom_depth = 2000.0 + config_sub_ice_shelf_2d_bottom_salinity = 34.7 + config_sub_ice_shelf_2d_cavity_thickness = 25.0 + config_sub_ice_shelf_2d_edge_width = 15.0e3 + config_sub_ice_shelf_2d_slope_height = 500.0 + config_sub_ice_shelf_2d_surface_salinity = 34.5 + config_sub_ice_shelf_2d_temperature = 1.0 + config_sub_ice_shelf_2d_vert_levels = 20 + config_sub_ice_shelf_2d_y1 = 30.0e3 + config_sub_ice_shelf_2d_y2 = 60.0e3 +/ +&periodic_planar + config_periodic_planar_bottom_depth = 2500.0 + config_periodic_planar_velocity_strength = 1.0 + config_periodic_planar_vert_levels = 100 +/ +&ecosys_column + config_ecosys_column_bottom_depth = 6000.0 + config_ecosys_column_ecosys_filename = 'unknown' + config_ecosys_column_ts_filename = 'unknown' + config_ecosys_column_vert_levels = 100 + config_ecosys_column_vertical_grid = '100layerACMEv1' +/ +&sea_mount + config_sea_mount_bottom_depth = 5000.0 + config_sea_mount_coriolis_parameter = -1.0e-4 + config_sea_mount_density_alpha = 0.2 + config_sea_mount_density_coef_exp = 1028 + config_sea_mount_density_coef_linear = 1024 + config_sea_mount_density_depth_exp = 500 + config_sea_mount_density_depth_linear = 4500 + config_sea_mount_density_gradient_exp = 3.0 + config_sea_mount_density_gradient_linear = 0.1 + config_sea_mount_density_ref = 1028 + config_sea_mount_density_tref = 5.0 + config_sea_mount_height = 4500.0 + config_sea_mount_layer_type = 'sigma' + config_sea_mount_radius = 10.0e3 + config_sea_mount_salinity = 35.0 + config_sea_mount_stratification_type = 'exponential' + config_sea_mount_vert_levels = 10 + config_sea_mount_width = 40.0e3 +/ +&isomip + config_isomip_bottom_depth = -900.0 + config_isomip_coriolis_parameter = -1.4e-4 + config_isomip_eastern_boundary = 500e3 + config_isomip_ice_fraction1 = 1.0 + config_isomip_ice_fraction2 = 1.0 + config_isomip_ice_fraction3 = 1.0 + config_isomip_northern_boundary = 1000e3 + config_isomip_restoring_salinity = 34.4 + config_isomip_restoring_temperature = -1.9 + config_isomip_salinity = 34.4 + config_isomip_salinity_piston_velocity = 1.157e-5 + config_isomip_southern_boundary = 0.0 + config_isomip_temperature = -1.9 + config_isomip_temperature_piston_velocity = 1.157e-5 + config_isomip_vert_levels = 30 + config_isomip_vertical_level_distribution = 'constant' + config_isomip_western_boundary = 0.0 + config_isomip_y1 = 0.0 + config_isomip_y2 = 400e3 + config_isomip_y3 = 1000e3 + config_isomip_z1 = -700.0 + config_isomip_z2 = -200.0 + config_isomip_z3 = -200.0 +/ +&isomip_plus + config_isomip_plus_coriolis_parameter = -1.409e-4 + config_isomip_plus_effective_density = 1026. + config_isomip_plus_init_bot_sal = 34.5 + config_isomip_plus_init_bot_temp = -1.9 + config_isomip_plus_init_top_sal = 33.8 + config_isomip_plus_init_top_temp = -1.9 + config_isomip_plus_max_bottom_depth = -720.0 + config_isomip_plus_min_column_thickness = 10.0 + config_isomip_plus_min_ocean_fraction = 0.5 + config_isomip_plus_minimum_levels = 3 + config_isomip_plus_restore_bot_sal = 34.7 + config_isomip_plus_restore_bot_temp = 1.0 + config_isomip_plus_restore_evap_rate = 200 + config_isomip_plus_restore_rate = 10.0 + config_isomip_plus_restore_top_sal = 33.8 + config_isomip_plus_restore_top_temp = -1.9 + config_isomip_plus_restore_xmax = 800.0e3 + config_isomip_plus_restore_xmin = 790.0e3 + config_isomip_plus_topography_file = 'input_geometry_processed.nc' + config_isomip_plus_vert_levels = 36 + config_isomip_plus_vertical_level_distribution = 'constant' +/ +&tracer_forcing_activetracers + config_salinity_restoring_constant_piston_velocity = 0.0 + config_salinity_restoring_max_difference = 0.5 + config_use_activetracers = .true. + config_use_activetracers_exponential_decay = .false. + config_use_activetracers_idealage_forcing = .false. + config_use_activetracers_interior_restoring = .false. + config_use_activetracers_surface_bulk_forcing = .true. + config_use_activetracers_surface_restoring = .false. + config_use_activetracers_ttd_forcing = .false. + config_use_surface_salinity_monthly_restoring = .false. +/ +&tracer_forcing_debugtracers + config_use_debugtracers = .false. + config_use_debugtracers_exponential_decay = .false. + config_use_debugtracers_idealage_forcing = .false. + config_use_debugtracers_interior_restoring = .false. + config_use_debugtracers_surface_bulk_forcing = .false. + config_use_debugtracers_surface_restoring = .false. + config_use_debugtracers_ttd_forcing = .false. +/ +&tracer_forcing_ecosystracers + config_use_ecosystracers = .false. + config_use_ecosystracers_exponential_decay = .false. + config_use_ecosystracers_idealage_forcing = .false. + config_use_ecosystracers_interior_restoring = .false. + config_use_ecosystracers_sea_ice_coupling = .false. + config_use_ecosystracers_surface_bulk_forcing = .false. + config_use_ecosystracers_surface_restoring = .false. + config_use_ecosystracers_surface_value = .false. + config_use_ecosystracers_ttd_forcing = .false. +/ +&tracer_forcing_dmstracers + config_use_dmstracers = .false. + config_use_dmstracers_exponential_decay = .false. + config_use_dmstracers_idealage_forcing = .false. + config_use_dmstracers_interior_restoring = .false. + config_use_dmstracers_sea_ice_coupling = .false. + config_use_dmstracers_surface_bulk_forcing = .false. + config_use_dmstracers_surface_restoring = .false. + config_use_dmstracers_surface_value = .false. + config_use_dmstracers_ttd_forcing = .false. +/ +&tracer_forcing_macromoleculestracers + config_use_macromoleculestracers = .false. + config_use_macromoleculestracers_exponential_decay = .false. + config_use_macromoleculestracers_idealage_forcing = .false. + config_use_macromoleculestracers_interior_restoring = .false. + config_use_macromoleculestracers_sea_ice_coupling = .false. + config_use_macromoleculestracers_surface_bulk_forcing = .false. + config_use_macromoleculestracers_surface_restoring = .false. + config_use_macromoleculestracers_surface_value = .false. + config_use_macromoleculestracers_ttd_forcing = .false. +/ +&am_globalstats + config_am_globalstats_compute_interval = 'output_interval' + config_am_globalstats_compute_on_startup = .true. + config_am_globalstats_directory = 'analysis_members' + config_am_globalstats_enable = .true. + config_am_globalstats_output_stream = 'globalStatsOutput' + config_am_globalstats_text_file = .false. + config_am_globalstats_write_on_startup = .true. +/ +&am_surfaceareaweightedaverages + config_am_surfaceareaweightedaverages_compute_interval = '0000-00-00_01:00:00' + config_am_surfaceareaweightedaverages_compute_on_startup = .true. + config_am_surfaceareaweightedaverages_enable = .true. + config_am_surfaceareaweightedaverages_output_stream = 'surfaceAreaWeightedAveragesOutput' + config_am_surfaceareaweightedaverages_write_on_startup = .true. +/ +&am_watermasscensus + config_am_watermasscensus_compute_interval = '0000-00-00_01:00:00' + config_am_watermasscensus_compute_on_startup = .true. + config_am_watermasscensus_enable = .false. + config_am_watermasscensus_maxsalinity = 37.0 + config_am_watermasscensus_maxtemperature = 30.0 + config_am_watermasscensus_minsalinity = 32.0 + config_am_watermasscensus_mintemperature = -2.0 + config_am_watermasscensus_output_stream = 'waterMassCensusOutput' + config_am_watermasscensus_write_on_startup = .true. +/ +&am_layervolumeweightedaverage + config_am_layervolumeweightedaverage_compute_interval = '0000-00-00_01:00:00' + config_am_layervolumeweightedaverage_compute_on_startup = .true. + config_am_layervolumeweightedaverage_enable = .true. + config_am_layervolumeweightedaverage_output_stream = 'layerVolumeWeightedAverageOutput' + config_am_layervolumeweightedaverage_write_on_startup = .true. +/ +&am_zonalmean + config_am_zonalmean_compute_interval = '0000-00-00_01:00:00' + config_am_zonalmean_compute_on_startup = .true. + config_am_zonalmean_enable = .false. + config_am_zonalmean_max_bin = -1.0e34 + config_am_zonalmean_min_bin = -1.0e34 + config_am_zonalmean_num_bins = 180 + config_am_zonalmean_output_stream = 'zonalMeanOutput' + config_am_zonalmean_write_on_startup = .true. +/ +&am_okuboweiss + config_am_okuboweiss_compute_eddy_census = .true. + config_am_okuboweiss_compute_interval = '0000-00-00_01:00:00' + config_am_okuboweiss_compute_on_startup = .true. + config_am_okuboweiss_directory = 'analysis_members' + config_am_okuboweiss_eddy_min_cells = 20 + config_am_okuboweiss_enable = .false. + config_am_okuboweiss_lambda2_normalization = 1e-10 + config_am_okuboweiss_normalization = 1e-10 + config_am_okuboweiss_output_stream = 'okuboWeissOutput' + config_am_okuboweiss_threshold_value = -0.2 + config_am_okuboweiss_use_lat_lon_coords = .true. + config_am_okuboweiss_write_on_startup = .true. +/ +&am_meridionalheattransport + config_am_meridionalheattransport_compute_interval = '0000-00-00_01:00:00' + config_am_meridionalheattransport_compute_on_startup = .true. + config_am_meridionalheattransport_enable = .true. + config_am_meridionalheattransport_max_bin = -1.0e34 + config_am_meridionalheattransport_min_bin = -1.0e34 + config_am_meridionalheattransport_num_bins = 180 + config_am_meridionalheattransport_output_stream = 'meridionalHeatTransportOutput' + config_am_meridionalheattransport_write_on_startup = .true. +/ +&am_testcomputeinterval + config_am_testcomputeinterval_compute_interval = '00-00-01_00:00:00' + config_am_testcomputeinterval_compute_on_startup = .true. + config_am_testcomputeinterval_enable = .false. + config_am_testcomputeinterval_output_stream = 'testComputeIntervalOutput' + config_am_testcomputeinterval_write_on_startup = .true. +/ +&am_highfrequencyoutput + config_am_highfrequencyoutput_compute_interval = 'output_interval' + config_am_highfrequencyoutput_compute_on_startup = .true. + config_am_highfrequencyoutput_enable = .false. + config_am_highfrequencyoutput_output_stream = 'highFrequencyOutput' + config_am_highfrequencyoutput_write_on_startup = .true. +/ +&am_timefilters + config_am_timefilters_compute_cell_centered_values = .true. + config_am_timefilters_compute_interval = 'dt' + config_am_timefilters_compute_on_startup = .true. + config_am_timefilters_enable = .false. + config_am_timefilters_initialize_filters = .true. + config_am_timefilters_output_stream = 'timeFiltersOutput' + config_am_timefilters_restart_stream = 'timeFiltersRestart' + config_am_timefilters_tau = '90_00:00:00' + config_am_timefilters_write_on_startup = .true. +/ +&am_lagrparttrack + config_am_lagrparttrack_compute_interval = 'dt' + config_am_lagrparttrack_compute_on_startup = .false. + config_am_lagrparttrack_enable = .false. + config_am_lagrparttrack_filter_number = 0 + config_am_lagrparttrack_input_stream = 'lagrPartTrackInput' + config_am_lagrparttrack_output_stream = 'lagrPartTrackOutput' + config_am_lagrparttrack_region_stream = 'lagrPartTrackRegions' + config_am_lagrparttrack_reset_criteria = 'none' + config_am_lagrparttrack_reset_global_timestamp = '0000_00:00:00' + config_am_lagrparttrack_reset_if_inside_region = .false. + config_am_lagrparttrack_reset_if_outside_region = .false. + config_am_lagrparttrack_restart_stream = 'lagrPartTrackRestart' + config_am_lagrparttrack_write_on_startup = .true. +/ +&am_eliassenpalm + config_am_eliassenpalm_compute_interval = 'output_interval' + config_am_eliassenpalm_compute_on_startup = .true. + config_am_eliassenpalm_debug = .false. + config_am_eliassenpalm_enable = .false. + config_am_eliassenpalm_nbuoyancylayers = 45 + config_am_eliassenpalm_output_stream = 'eliassenPalmOutput' + config_am_eliassenpalm_restart_stream = 'eliassenPalmRestart' + config_am_eliassenpalm_rhomax_buoycoor = 1080 + config_am_eliassenpalm_rhomin_buoycoor = 900 + config_am_eliassenpalm_write_on_startup = .true. +/ +&am_mixedlayerdepths + config_am_mixedlayerdepths_compute_interval = '0000-00-00_01:00:00' + config_am_mixedlayerdepths_compute_on_startup = .true. + config_am_mixedlayerdepths_crit_dens_threshold = 0.03 + config_am_mixedlayerdepths_crit_temp_threshold = 0.2 + config_am_mixedlayerdepths_den_gradient_threshold = 5E-8 + config_am_mixedlayerdepths_dgradient = .true. + config_am_mixedlayerdepths_dthreshold = .true. + config_am_mixedlayerdepths_enable = .true. + config_am_mixedlayerdepths_interp_method = 1 + config_am_mixedlayerdepths_output_stream = 'mixedLayerDepthsOutput' + config_am_mixedlayerdepths_reference_pressure = 1.0E5 + config_am_mixedlayerdepths_temp_gradient_threshold = 5E-7 + config_am_mixedlayerdepths_tgradient = .true. + config_am_mixedlayerdepths_tthreshold = .true. + config_am_mixedlayerdepths_write_on_startup = .true. +/ +&am_regionalstatsdaily + config_am_regionalstatsdaily_1d_weighting_field = 'areaCell' + config_am_regionalstatsdaily_1d_weighting_function = 'mul' + config_am_regionalstatsdaily_2d_weighting_field = 'volumeCell' + config_am_regionalstatsdaily_2d_weighting_function = 'mul' + config_am_regionalstatsdaily_compute_interval = 'output_interval' + config_am_regionalstatsdaily_compute_on_startup = .false. + config_am_regionalstatsdaily_enable = .false. + config_am_regionalstatsdaily_input_stream = 'regionalMasksInput' + config_am_regionalstatsdaily_operation = 'avg' + config_am_regionalstatsdaily_output_stream = 'regionalStatsDailyOutput' + config_am_regionalstatsdaily_region_group = 'all' + config_am_regionalstatsdaily_region_type = 'cell' + config_am_regionalstatsdaily_restart_stream = 'regionalMasksInput' + config_am_regionalstatsdaily_vertical_dimension = 'nVertLevels' + config_am_regionalstatsdaily_vertical_mask = 'cellMask' + config_am_regionalstatsdaily_write_on_startup = .false. +/ +&am_regionalstatsweekly + config_am_regionalstatsweekly_1d_weighting_field = 'areaCell' + config_am_regionalstatsweekly_1d_weighting_function = 'mul' + config_am_regionalstatsweekly_2d_weighting_field = 'volumeCell' + config_am_regionalstatsweekly_2d_weighting_function = 'mul' + config_am_regionalstatsweekly_compute_interval = 'output_interval' + config_am_regionalstatsweekly_compute_on_startup = .false. + config_am_regionalstatsweekly_enable = .false. + config_am_regionalstatsweekly_input_stream = 'regionalMasksInput' + config_am_regionalstatsweekly_operation = 'avg' + config_am_regionalstatsweekly_output_stream = 'regionalStatsWeeklyOutput' + config_am_regionalstatsweekly_region_group = 'all' + config_am_regionalstatsweekly_region_type = 'cell' + config_am_regionalstatsweekly_restart_stream = 'regionalMasksInput' + config_am_regionalstatsweekly_vertical_dimension = 'nVertLevels' + config_am_regionalstatsweekly_vertical_mask = 'cellMask' + config_am_regionalstatsweekly_write_on_startup = .false. +/ +&am_regionalstatsmonthly + config_am_regionalstatsmonthly_1d_weighting_field = 'areaCell' + config_am_regionalstatsmonthly_1d_weighting_function = 'mul' + config_am_regionalstatsmonthly_2d_weighting_field = 'volumeCell' + config_am_regionalstatsmonthly_2d_weighting_function = 'mul' + config_am_regionalstatsmonthly_compute_interval = 'output_interval' + config_am_regionalstatsmonthly_compute_on_startup = .false. + config_am_regionalstatsmonthly_enable = .false. + config_am_regionalstatsmonthly_input_stream = 'regionalMasksInput' + config_am_regionalstatsmonthly_operation = 'avg' + config_am_regionalstatsmonthly_output_stream = 'regionalStatsMonthlyOutput' + config_am_regionalstatsmonthly_region_group = 'all' + config_am_regionalstatsmonthly_region_type = 'cell' + config_am_regionalstatsmonthly_restart_stream = 'regionalMasksInput' + config_am_regionalstatsmonthly_vertical_dimension = 'nVertLevels' + config_am_regionalstatsmonthly_vertical_mask = 'cellMask' + config_am_regionalstatsmonthly_write_on_startup = .false. +/ +&am_regionalstatscustom + config_am_regionalstatscustom_1d_weighting_field = 'areaCell' + config_am_regionalstatscustom_1d_weighting_function = 'mul' + config_am_regionalstatscustom_2d_weighting_field = 'volumeCell' + config_am_regionalstatscustom_2d_weighting_function = 'mul' + config_am_regionalstatscustom_compute_interval = 'output_interval' + config_am_regionalstatscustom_compute_on_startup = .false. + config_am_regionalstatscustom_enable = .false. + config_am_regionalstatscustom_input_stream = 'regionalMasksInput' + config_am_regionalstatscustom_operation = 'avg' + config_am_regionalstatscustom_output_stream = 'regionalStatsCustomOutput' + config_am_regionalstatscustom_region_group = 'all' + config_am_regionalstatscustom_region_type = 'cell' + config_am_regionalstatscustom_restart_stream = 'regionalMasksInput' + config_am_regionalstatscustom_vertical_dimension = 'nVertLevels' + config_am_regionalstatscustom_vertical_mask = 'cellMask' + config_am_regionalstatscustom_write_on_startup = .false. +/ +&am_timeseriesstatsdaily + config_am_timeseriesstatsdaily_backward_output_offset = '00-00-01_00:00:00' + config_am_timeseriesstatsdaily_compute_interval = '00-00-00_01:00:00' + config_am_timeseriesstatsdaily_compute_on_startup = .false. + config_am_timeseriesstatsdaily_duration_intervals = 'repeat_interval' + config_am_timeseriesstatsdaily_enable = .false. + config_am_timeseriesstatsdaily_operation = 'avg' + config_am_timeseriesstatsdaily_output_stream = 'timeSeriesStatsDailyOutput' + config_am_timeseriesstatsdaily_reference_times = 'initial_time' + config_am_timeseriesstatsdaily_repeat_intervals = 'reset_interval' + config_am_timeseriesstatsdaily_reset_intervals = '00-00-01_00:00:00' + config_am_timeseriesstatsdaily_restart_stream = 'timeSeriesStatsDailyRestart' + config_am_timeseriesstatsdaily_write_on_startup = .false. +/ +&am_timeseriesstatsmonthly + config_am_timeseriesstatsmonthly_backward_output_offset = '00-01-00_00:00:00' + config_am_timeseriesstatsmonthly_compute_interval = '00-00-00_01:00:00' + config_am_timeseriesstatsmonthly_compute_on_startup = .false. + config_am_timeseriesstatsmonthly_duration_intervals = 'repeat_interval' + config_am_timeseriesstatsmonthly_enable = .true. + config_am_timeseriesstatsmonthly_operation = 'avg' + config_am_timeseriesstatsmonthly_output_stream = 'timeSeriesStatsMonthlyOutput' + config_am_timeseriesstatsmonthly_reference_times = 'initial_time' + config_am_timeseriesstatsmonthly_repeat_intervals = 'reset_interval' + config_am_timeseriesstatsmonthly_reset_intervals = '00-01-00_00:00:00' + config_am_timeseriesstatsmonthly_restart_stream = 'timeSeriesStatsMonthlyRestart' + config_am_timeseriesstatsmonthly_write_on_startup = .true. +/ +&am_timeseriesstatsclimatology + config_am_timeseriesstatsclimatology_backward_output_offset = '00-03-00_00:00:00' + config_am_timeseriesstatsclimatology_compute_interval = '00-00-00_01:00:00' + config_am_timeseriesstatsclimatology_compute_on_startup = .false. + config_am_timeseriesstatsclimatology_duration_intervals = '00-03-00_00:00:00;00-03-00_00:00:00;00-03-00_00:00:00;00-03-00_00:00:00' + config_am_timeseriesstatsclimatology_enable = .false. + config_am_timeseriesstatsclimatology_operation = 'avg' + config_am_timeseriesstatsclimatology_output_stream = 'timeSeriesStatsClimatologyOutput' + config_am_timeseriesstatsclimatology_reference_times = '00-03-01_00:00:00;00-06-01_00:00:00;00-09-01_00:00:00;00-12-01_00:00:00' + config_am_timeseriesstatsclimatology_repeat_intervals = '01-00-00_00:00:00;01-00-00_00:00:00;01-00-00_00:00:00;01-00-00_00:00:00' + config_am_timeseriesstatsclimatology_reset_intervals = '1000-00-00_00:00:00;1000-00-00_00:00:00;1000-00-00_00:00:00;1000-00-00_00:00:00' + config_am_timeseriesstatsclimatology_restart_stream = 'timeSeriesStatsClimatologyRestart' + config_am_timeseriesstatsclimatology_write_on_startup = .false. +/ +&am_timeseriesstatscustom + config_am_timeseriesstatscustom_backward_output_offset = '00-00-01_00:00:00' + config_am_timeseriesstatscustom_compute_interval = '00-00-00_01:00:00' + config_am_timeseriesstatscustom_compute_on_startup = .false. + config_am_timeseriesstatscustom_duration_intervals = 'repeat_interval' + config_am_timeseriesstatscustom_enable = .false. + config_am_timeseriesstatscustom_operation = 'avg' + config_am_timeseriesstatscustom_output_stream = 'timeSeriesStatsCustomOutput' + config_am_timeseriesstatscustom_reference_times = 'initial_time' + config_am_timeseriesstatscustom_repeat_intervals = 'reset_interval' + config_am_timeseriesstatscustom_reset_intervals = '00-00-07_00:00:00' + config_am_timeseriesstatscustom_restart_stream = 'timeSeriesStatsCustomRestart' + config_am_timeseriesstatscustom_write_on_startup = .false. +/ +&am_pointwisestats + config_am_pointwisestats_compute_interval = 'output_interval' + config_am_pointwisestats_compute_on_startup = .true. + config_am_pointwisestats_enable = .false. + config_am_pointwisestats_output_stream = 'pointwiseStatsOutput' + config_am_pointwisestats_write_on_startup = .true. +/ +&am_debugdiagnostics + config_am_debugdiagnostics_check_state = .true. + config_am_debugdiagnostics_compute_interval = 'dt' + config_am_debugdiagnostics_compute_on_startup = .true. + config_am_debugdiagnostics_enable = .false. + config_am_debugdiagnostics_output_stream = 'debugDiagnosticsOutput' + config_am_debugdiagnostics_write_on_startup = .false. +/ +&am_rpncalculator + config_am_rpncalculator_compute_interval = '0010-00-00_00:00:00' + config_am_rpncalculator_compute_on_startup = .true. + config_am_rpncalculator_enable = .false. + config_am_rpncalculator_expression_1 = 'a b *' + config_am_rpncalculator_expression_2 = 'none' + config_am_rpncalculator_expression_3 = 'none' + config_am_rpncalculator_expression_4 = 'none' + config_am_rpncalculator_output_name_1 = 'volumeCell' + config_am_rpncalculator_output_name_2 = 'none' + config_am_rpncalculator_output_name_3 = 'none' + config_am_rpncalculator_output_name_4 = 'none' + config_am_rpncalculator_output_stream = 'none' + config_am_rpncalculator_variable_a = 'layerThickness' + config_am_rpncalculator_variable_b = 'areaCell' + config_am_rpncalculator_variable_c = 'none' + config_am_rpncalculator_variable_d = 'none' + config_am_rpncalculator_variable_e = 'none' + config_am_rpncalculator_variable_f = 'none' + config_am_rpncalculator_variable_g = 'none' + config_am_rpncalculator_variable_h = 'none' + config_am_rpncalculator_write_on_startup = .false. +/ +&am_transecttransport + config_am_transecttransport_compute_interval = 'output_interval' + config_am_transecttransport_compute_on_startup = .true. + config_am_transecttransport_enable = .false. + config_am_transecttransport_output_stream = 'transectTransportOutput' + config_am_transecttransport_transect_group = 'all' + config_am_transecttransport_write_on_startup = .true. +/ +&am_eddyproductvariables + config_am_eddyproductvariables_compute_interval = 'dt' + config_am_eddyproductvariables_compute_on_startup = .true. + config_am_eddyproductvariables_enable = .false. + config_am_eddyproductvariables_output_stream = 'eddyProductVariablesOutput' + config_am_eddyproductvariables_write_on_startup = .false. +/ +&am_mocstreamfunction + config_am_mocstreamfunction_compute_interval = 'output_interval' + config_am_mocstreamfunction_compute_on_startup = .true. + config_am_mocstreamfunction_enable = .false. + config_am_mocstreamfunction_max_bin = -1.0e34 + config_am_mocstreamfunction_min_bin = -1.0e34 + config_am_mocstreamfunction_normal_velocity_value = 'normalVelocity' + config_am_mocstreamfunction_num_bins = 180 + config_am_mocstreamfunction_output_stream = 'mocStreamfunctionOutput' + config_am_mocstreamfunction_region_group = 'all' + config_am_mocstreamfunction_transect_group = 'all' + config_am_mocstreamfunction_vertical_velocity_value = 'vertVelocityTop' + config_am_mocstreamfunction_write_on_startup = .true. +/ diff --git a/mpas_analysis/test/test_climatology/streams.ocean b/mpas_analysis/test/test_climatology/streams.ocean new file mode 100644 index 000000000..599918fe1 --- /dev/null +++ b/mpas_analysis/test/test_climatology/streams.ocean @@ -0,0 +1,16 @@ + + + + + + + + + diff --git a/mpas_analysis/test/test_interpolate.py b/mpas_analysis/test/test_interpolate.py index 0008459c4..b1fb23976 100644 --- a/mpas_analysis/test/test_interpolate.py +++ b/mpas_analysis/test/test_interpolate.py @@ -42,8 +42,9 @@ def get_mpas_descriptor(self): def get_latlon_file_descriptor(self): latLonGridFileName = str(self.datadir.join('SST_annual_1870-1900.nc')) - descriptor = LatLonGridDescriptor() - descriptor.read(latLonGridFileName, latVarName='lat', lonVarName='lon') + descriptor = LatLonGridDescriptor.read(latLonGridFileName, + latVarName='lat', + lonVarName='lon') return (descriptor, latLonGridFileName) @@ -57,8 +58,7 @@ def get_latlon_array_descriptor(self): lon = numpy.array(config.getExpression('interpolate', 'lon', usenumpyfunc=True)) - descriptor = LatLonGridDescriptor() - descriptor.create(lat, lon, units='degrees') + descriptor = LatLonGridDescriptor.create(lat, lon, units='degrees') return descriptor def get_stereographic_array_descriptor(self): @@ -73,9 +73,9 @@ def get_stereographic_array_descriptor(self): res = 100e3 nx = 2*int(xMax/res)+1 x = numpy.linspace(-xMax, xMax, nx) - descriptor = ProjectionGridDescriptor(projection) meshName = '{}km_Antarctic_stereo'.format(int(res*1e-3)) - descriptor.create(x, x, meshName) + descriptor = ProjectionGridDescriptor.create(projection, x, x, + meshName) return descriptor def get_file_names(self, suffix): diff --git a/run_analysis.py b/run_analysis.py index 2905c5b61..ed38f4d3b 100755 --- a/run_analysis.py +++ b/run_analysis.py @@ -17,6 +17,7 @@ import warnings import subprocess import time +from collections import OrderedDict from mpas_analysis.configuration.MpasAnalysisConfigParser \ import MpasAnalysisConfigParser @@ -25,12 +26,139 @@ make_directories +def build_analysis_list(config, isSubtask): # {{{ + """ + Build a list of analysis tasks based on the 'generate' config option. + + Authors + ------- + Xylar Asay-Davis + """ + + # choose the right rendering backend, depending on whether we're displaying + # to the screen + if not config.getboolean('plot', 'displayToScreen'): + mpl.use('Agg') + + # analysis can only be imported after the right MPL renderer is selected + from mpas_analysis import ocean + from mpas_analysis import sea_ice + from mpas_analysis.shared.cache_dataset_times_task \ + import CacheDatasetTimesTask + + # analyses will be a list of analysis classes + analyses = [] + + # add cacheOceanTimeSeriesStatsTimes task for caching the times in + # MPAS-Ocean timeSeriesStatsMonthly output files. This is a prerequisite + # for all ocean analysis. + analyses.append(CacheDatasetTimesTask( + config=config, + componentName='ocean', + streamName='timeSeriesStats', + startAndEndDateSections=['climatology', 'timeSeries', 'index'], + namelistOption='config_am_timeseriesstatsmonthly_enable')) + + # add cacheSeaIceTimeSeriesStatsTimes task for caching the times in + # MPAS-Ocean timeSeriesStatsMonthly output files. This is a prerequisite + # for all sea ice analysis. + analyses.append(CacheDatasetTimesTask( + config=config, + componentName='seaIce', + streamName='timeSeriesStats', + startAndEndDateSections=['climatology', 'timeSeries'], + namelistOption='config_am_timeseriesstatsmonthly_enable')) + + # Ocean Analyses + analyses.append(ocean.TimeSeriesOHC(config)) + analyses.append(ocean.TimeSeriesSST(config)) + analyses.append(ocean.IndexNino34(config)) + analyses.append(ocean.MeridionalHeatTransport(config)) + analyses.append(ocean.StreamfunctionMOC(config)) + + analyses.append(ocean.ClimatologyMapSST(config)) + analyses.append(ocean.ClimatologyMapMLD(config)) + analyses.append(ocean.ClimatologyMapSSS(config)) + + # Sea Ice Analyses + analyses.append(sea_ice.TimeSeriesSeaIce(config)) + analyses.append(sea_ice.ClimatologyMapSeaIce(config)) + + possibleAnalyses = OrderedDict() + for analysisTask in analyses: + possibleAnalyses[analysisTask.taskName] = analysisTask + + # check which analysis we actually want to generate and only keep those + analysesToGenerate = OrderedDict() + for analysisTask in possibleAnalyses.itervalues(): + # update the dictionary with this task and perhaps its prerequisites + analysesToAdd = add_task_and_prereqisites(analysisTask, + possibleAnalyses, + analysesToGenerate, + isPrerequisite=False, + isSubtask=isSubtask) + analysesToGenerate.update(analysesToAdd) + + return analysesToGenerate # }}} + + +def add_task_and_prereqisites(analysisTask, possibleAnalyses, + analysesToGenerate, isPrerequisite, + isSubtask): # {{{ + """ + If a task has been requested through the generate config option or + if it is a prerequisite of a requested task, add it to the dictionary of + tasks to generate. + + Authors + ------- + Xylar Asay-Davis + """ + + analysesToAdd = OrderedDict() + # for each anlaysis task, check if we want to generate this task + # and if the analysis task has a valid configuration + if isPrerequisite or analysisTask.check_generate(): + add = False + try: + analysisTask.setup_and_check() + add = True + except: + traceback.print_exc(file=sys.stdout) + print "ERROR: analysis task {} failed during check and " \ + "will not be run".format(analysisTask.taskName) + if add and not isSubtask: + # first, we should try to add the prerequisites + prereqs = analysisTask.prerequisiteTasks + if prereqs is not None: + for prereq in prereqs: + if prereq not in analysesToGenerate.keys(): + prereqToAdd = add_task_and_prereqisites( + possibleAnalyses[prereq], possibleAnalyses, + analysesToGenerate, isPrerequisite=True) + if len(prereqToAdd.keys()) == 0: + # a prerequisite failed setup_and_check + print "ERROR: a prerequisite of analysis task {}" \ + " failed during check and will not be" \ + " run".format(analysisTask.taskName) + add = False + break + # the prerequisite (and its prerequisites) should be + # added + analysesToAdd.update(prereqToAdd) + if add: + analysesToAdd[analysisTask.taskName] = analysisTask + + return analysesToAdd # }}} + + def update_generate(config, generate): # {{{ """ Update the 'generate' config option using a string from the command line. - Author: Xylar Asay-Davis - Last Modified: 03/07/2017 + Authors + ------- + Xylar Asay-Davis """ # overwrite the 'generate' in config with a string that parses to @@ -42,8 +170,7 @@ def update_generate(config, generate): # {{{ config.set('output', 'generate', generateString) # }}} -def run_parallel_tasks(config, analyses, configFiles, taskCount): - # {{{ +def run_analysis(config, analyses, configFiles, isSubtask): # {{{ """ Run this script once each for several parallel tasks. @@ -52,40 +179,99 @@ def run_parallel_tasks(config, analyses, configFiles, taskCount): Xylar Asay-Davis """ - taskNames = [analysisTask.taskName for analysisTask in analyses] + taskCount = config.getWithDefault('execute', 'parallelTaskCount', + default=1) - taskCount = min(taskCount, len(taskNames)) + isParallel = not isSubtask and taskCount > 1 and len(analyses) > 1 - (processes, logs) = launch_tasks(taskNames[0:taskCount], config, - configFiles) - remainingTasks = taskNames[taskCount:] - while len(processes) > 0: - (taskName, process) = wait_for_task(processes) - if process.returncode == 0: - print "Task {} has finished successfully.".format(taskName) + for analysisTask in analyses.itervalues(): + if analysisTask.prerequisiteTasks is None or isSubtask: + analysisTask.status = 'ready' else: - print "ERROR in task {}. See log file {} for details".format( - taskName, logs[taskName].name) - logs[taskName].close() - # remove the process from the process dictionary (no need to bother) - processes.pop(taskName) - - if len(remainingTasks) > 0: - (process, log) = launch_tasks(remainingTasks[0:1], config, - configFiles) - # merge the new process and log into these dictionaries - processes.update(process) - logs.update(log) - remainingTasks = remainingTasks[1:] + analysisTask.status = 'blocked' + + processes = {} + logs = {} + + # run each analysis task + lastException = None + + runningCount = 0 + while True: + # we still have tasks to run + for analysisTask in analyses.itervalues(): + if analysisTask.status == 'blocked': + prereqStatus = [analyses[prereq].status for prereq in + analysisTask.prerequisiteTasks] + if any([status == 'fail' for status in prereqStatus]): + # a prerequisite failed so this task cannot succeed + analysisTask.status = 'fail' + if all([status == 'success' for status in prereqStatus]): + # no unfinished prerequisites so we can run this task + analysisTask.status = 'ready' + + unfinishedCount = 0 + for analysisTask in analyses.itervalues(): + if analysisTask.status not in ['success', 'fail']: + unfinishedCount += 1 + + if unfinishedCount <= 0: + # we're done + break + + # launch new tasks + for taskName, analysisTask in analyses.items(): + if analysisTask.status == 'ready': + if isParallel: + process, logFile = launch_task(taskName, config, + configFiles) + processes[taskName] = process + logs[taskName] = logFile + analysisTask.status = 'running' + runningCount += 1 + if runningCount >= taskCount: + break + else: + exception = run_task(config, analysisTask) + if exception is None: + analysisTask.status = 'success' + else: + lastException = exception + analysisTask.status = 'fail' + + if isParallel: + # wait for a task to finish + (taskName, process) = wait_for_task(processes) + analysisTask = analyses[taskName] + runningCount -= 1 + processes.pop(taskName) + if process.returncode == 0: + print "Task {} has finished successfully.".format(taskName) + analysisTask.status = 'success' + else: + print "ERROR in task {}. See log file {} for details".format( + taskName, logs[taskName].name) + analysisTask.status = 'fail' + logs[taskName].close() + + if not isParallel and config.getboolean('plot', 'displayToScreen'): + import matplotlib.pyplot as plt + plt.show() + + # raise the last exception so the process exits with an error + if lastException is not None: + raise lastException + # }}} -def launch_tasks(taskNames, config, configFiles): # {{{ +def launch_task(taskName, config, configFiles): # {{{ """ - Launch one or more tasks + Launch a parallel tasks - Author: Xylar Asay-Davis - Last Modified: 03/08/2017 + Authors + ------- + Xylar Asay-Davis """ thisFile = os.path.realpath(__file__) @@ -96,25 +282,21 @@ def launch_tasks(taskNames, config, configFiles): # {{{ else: commandPrefix = commandPrefix.split(' ') - processes = {} - logs = {} - for taskName in taskNames: - args = commandPrefix + [thisFile, '--generate', taskName] + configFiles + args = commandPrefix + [thisFile, '--subtask', '--generate', taskName] \ + + configFiles - logFileName = '{}/{}.log'.format(logsDirectory, taskName) + logFileName = '{}/{}.log'.format(logsDirectory, taskName) - # write the command to the log file - logFile = open(logFileName, 'w') - logFile.write('Command: {}\n'.format(' '.join(args))) - # make sure the command gets written before the rest of the log - logFile.flush() - print 'Running {}'.format(taskName) - process = subprocess.Popen(args, stdout=logFile, - stderr=subprocess.STDOUT) - processes[taskName] = process - logs[taskName] = logFile + # write the command to the log file + logFile = open(logFileName, 'w') + logFile.write('Command: {}\n'.format(' '.join(args))) + # make sure the command gets written before the rest of the log + logFile.flush() + print 'Running {}'.format(taskName) + process = subprocess.Popen(args, stdout=logFile, + stderr=subprocess.STDOUT) - return (processes, logs) # }}} + return (process, logFile) # }}} def wait_for_task(processes): # {{{ @@ -122,8 +304,9 @@ def wait_for_task(processes): # {{{ Wait for the next process to finish and check its status. Returns both the task name and the process that finished. - Author: Xylar Asay-Davis - Last Modified: 03/08/2017 + Authors + ------- + Xylar Asay-Davis """ # first, check if any process has already finished @@ -144,8 +327,9 @@ def is_running(process): # {{{ """ Returns whether a given process is currently running - Author: Xylar Asay-Davis - Last Modified: 03/08/2017 + Authors + ------- + Xylar Asay-Davis """ try: @@ -156,108 +340,56 @@ def is_running(process): # {{{ return True # }}} -def build_analysis_list(config): # {{{ +def run_task(config, analysisTask): # {{{ """ - Build a list of analysis modules based on the 'generate' config option. + Run a single analysis task, time the task, write out the config file + (including any modifications specific to the task) and return the exception + raised (if any) - Author: Xylar Asay-Davis - Last Modified: 03/07/2017 + Authors + ------- + Xylar Asay-Davis """ - # choose the right rendering backend, depending on whether we're displaying - # to the screen - if not config.getboolean('plot', 'displayToScreen'): - mpl.use('Agg') - - # analysis can only be imported after the right MPL renderer is selected - from mpas_analysis import ocean - from mpas_analysis import sea_ice - - # analyses will be a list of analysis classes - analyses = [] - - # Ocean Analyses - analyses.append(ocean.TimeSeriesOHC(config)) - analyses.append(ocean.TimeSeriesSST(config)) - analyses.append(ocean.IndexNino34(config)) - analyses.append(ocean.MeridionalHeatTransport(config)) - analyses.append(ocean.StreamfunctionMOC(config)) - - analyses.append(ocean.ClimatologyMapSST(config)) - analyses.append(ocean.ClimatologyMapMLD(config)) - analyses.append(ocean.ClimatologyMapSSS(config)) - - # Sea Ice Analyses - analyses.append(sea_ice.TimeSeriesSeaIce(config)) - analyses.append(sea_ice.ClimatologyMapSeaIce(config)) - - # check which analysis we actually want to generate and only keep those - analysesToGenerate = [] - for analysisTask in analyses: - # for each anlaysis module, check if we want to generate this task - # and if the analysis task has a valid configuration - if analysisTask.check_generate(): - add = False - try: - analysisTask.setup_and_check() - add = True - except: - traceback.print_exc(file=sys.stdout) - print "ERROR: analysis module {} failed during check and " \ - "will not be run".format(analysisTask.taskName) - if add: - analysesToGenerate.append(analysisTask) - - return analysesToGenerate # }}} - - -def run_analysis(config, analyses): # {{{ - - # run each analysis task - lastException = None - for analysisTask in analyses: - # write out a copy of the configuration to document the run - logsDirectory = build_config_full_path(config, 'output', - 'logsSubdirectory') - try: - startTime = time.clock() - analysisTask.run() - runDuration = time.clock() - startTime - m, s = divmod(runDuration, 60) - h, m = divmod(int(m), 60) - print 'Execution time: {}:{:02d}:{:05.2f}'.format(h, m, s) - except (Exception, BaseException) as e: - if isinstance(e, KeyboardInterrupt): - raise e - traceback.print_exc(file=sys.stdout) - print "ERROR: analysis module {} failed during run".format( - analysisTask.taskName) - lastException = e - - configFileName = '{}/configs/config.{}'.format(logsDirectory, - analysisTask.taskName) - configFile = open(configFileName, 'w') - config.write(configFile) - configFile.close() - - if config.getboolean('plot', 'displayToScreen'): - import matplotlib.pyplot as plt - plt.show() - - # raise the last exception so the process exits with an error - if lastException is not None: - raise lastException - - return # }}} + # write out a copy of the configuration to document the run + logsDirectory = build_config_full_path(config, 'output', + 'logsSubdirectory') + exception = None + try: + startTime = time.clock() + analysisTask.run() + runDuration = time.clock() - startTime + m, s = divmod(runDuration, 60) + h, m = divmod(int(m), 60) + print 'Execution time: {}:{:02d}:{:05.2f}'.format(h, m, s) + except (Exception, BaseException) as e: + if isinstance(e, KeyboardInterrupt): + raise e + traceback.print_exc(file=sys.stdout) + print "ERROR: analysis task {} failed during run".format( + analysisTask.taskName) + exception = e + + configFileName = '{}/configs/config.{}'.format(logsDirectory, + analysisTask.taskName) + configFile = open(configFileName, 'w') + config.write(configFile) + configFile.close() + + return exception # }}} if __name__ == "__main__": parser = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.RawTextHelpFormatter) + parser.add_argument("--subtask", dest="subtask", action='store_true', + help="If this is a subtask when running parallel " + "tasks") parser.add_argument("-g", "--generate", dest="generate", help="A list of analysis modules to generate " - "(nearly identical generate option in config file).", + "(nearly identical generate option in config " + "file).", metavar="ANALYSIS1[,ANALYSIS2,ANALYSIS3,...]") parser.add_argument('configFiles', metavar='CONFIG', type=str, nargs='+', help='config file') @@ -286,14 +418,8 @@ def run_analysis(config, analyses): # {{{ make_directories(logsDirectory) make_directories('{}/configs/'.format(logsDirectory)) - analyses = build_analysis_list(config) - - parallelTaskCount = config.getWithDefault('execute', 'parallelTaskCount', - default=1) + analyses = build_analysis_list(config, args.subtask) - if parallelTaskCount <= 1 or len(analyses) == 1: - run_analysis(config, analyses) - else: - run_parallel_tasks(config, analyses, configFiles, parallelTaskCount) + run_analysis(config, analyses, configFiles, args.subtask) # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python