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