diff --git a/configs/anvil/config.20170926.FCT2.A_WCYCL1850S.ne30_oECv3.anvil b/configs/anvil/config.20170926.FCT2.A_WCYCL1850S.ne30_oECv3.anvil index 7892f47c2..2768885a2 100644 --- a/configs/anvil/config.20170926.FCT2.A_WCYCL1850S.ne30_oECv3.anvil +++ b/configs/anvil/config.20170926.FCT2.A_WCYCL1850S.ne30_oECv3.anvil @@ -11,6 +11,13 @@ mainRunName = 20170926.FCT2.A_WCYCL1850S.ne30_oECv3.anvil # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [input] ## options related to reading in the results to be analyzed diff --git a/configs/cori/config.20171201.default.GMPAS-IAF.T62_oRRS30to10v3wLI.cori-knl b/configs/cori/config.20171201.default.GMPAS-IAF.T62_oRRS30to10v3wLI.cori-knl index d21913142..d7e31c660 100644 --- a/configs/cori/config.20171201.default.GMPAS-IAF.T62_oRRS30to10v3wLI.cori-knl +++ b/configs/cori/config.20171201.default.GMPAS-IAF.T62_oRRS30to10v3wLI.cori-knl @@ -11,6 +11,13 @@ mainRunName = 20171201.default.GMPAS-IAF.T62_oRRS30to10v3wLI.cori-knl # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [execute] ## options related to executing parallel tasks diff --git a/configs/cori/job_script.cori-haswell.bash b/configs/cori/job_script.cori-haswell.bash index ec4d6b79c..aefa4fb46 100644 --- a/configs/cori/job_script.cori-haswell.bash +++ b/configs/cori/job_script.cori-haswell.bash @@ -32,10 +32,10 @@ if [ ! -f $run_config_file ]; then echo "File $run_config_file not found!" exit 1 fi -if [ ! -f ./run_analysis.py ]; then - echo "run_analysis.py not found in current directory!" +if [ ! -f ./run_mpas_analysis ]; then + echo "run_mpas_analysis not found in current directory!" exit 1 fi -srun -N 1 -n 1 ./run_analysis.py $run_config_file +srun -N 1 -n 1 ./run_mpas_analysis $run_config_file diff --git a/configs/cori/job_script.cori-knl.bash b/configs/cori/job_script.cori-knl.bash index 212af40cc..b3936933b 100644 --- a/configs/cori/job_script.cori-knl.bash +++ b/configs/cori/job_script.cori-knl.bash @@ -32,10 +32,10 @@ if [ ! -f $run_config_file ]; then echo "File $run_config_file not found!" exit 1 fi -if [ ! -f ./run_analysis.py ]; then - echo "run_analysis.py not found in current directory!" +if [ ! -f ./run_mpas_analysis ]; then + echo "run_mpas_analysis not found in current directory!" exit 1 fi -srun -N 1 -n 1 ./run_analysis.py $run_config_file +srun -N 1 -n 1 ./run_mpas_analysis $run_config_file diff --git a/configs/edison/config.20170807.beta1.G_oQU240.edison b/configs/edison/config.20170807.beta1.G_oQU240.edison index 161c18819..15176e615 100644 --- a/configs/edison/config.20170807.beta1.G_oQU240.edison +++ b/configs/edison/config.20170807.beta1.G_oQU240.edison @@ -11,6 +11,13 @@ mainRunName = 20170807.beta1.G_oQU240.edison # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [input] ## options related to reading in the results to be analyzed diff --git a/configs/edison/config.20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison b/configs/edison/config.20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison index fb89b8531..47faf9573 100644 --- a/configs/edison/config.20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison +++ b/configs/edison/config.20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison @@ -12,6 +12,13 @@ mainRunName = 20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [input] ## options related to reading in the results to be analyzed diff --git a/configs/edison/config.20171102.beta3rc02_1850.ne30_oECv3_ICG.edison b/configs/edison/config.20171102.beta3rc02_1850.ne30_oECv3_ICG.edison index 18fbffb70..c98f76515 100644 --- a/configs/edison/config.20171102.beta3rc02_1850.ne30_oECv3_ICG.edison +++ b/configs/edison/config.20171102.beta3rc02_1850.ne30_oECv3_ICG.edison @@ -12,6 +12,13 @@ mainRunName = 20171102.beta3rc02_1850.ne30_oECv3_ICG.edison # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [input] ## options related to reading in the results to be analyzed diff --git a/configs/edison/config.B_low_res_ice_shelves_1696_JWolfe_layout_Edison b/configs/edison/config.B_low_res_ice_shelves_1696_JWolfe_layout_Edison index 7d4bbeaed..8a76e16bd 100644 --- a/configs/edison/config.B_low_res_ice_shelves_1696_JWolfe_layout_Edison +++ b/configs/edison/config.B_low_res_ice_shelves_1696_JWolfe_layout_Edison @@ -11,6 +11,13 @@ mainRunName = B_low_res_ice_shelves_1696_JWolfe_layout_Edison # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [execute] ## options related to executing parallel tasks diff --git a/configs/lanl/config.20170207.MPAS-SeaIce.QU60km_polar.wolf b/configs/lanl/config.20170207.MPAS-SeaIce.QU60km_polar.wolf index 3e5d25354..51cf6628f 100644 --- a/configs/lanl/config.20170207.MPAS-SeaIce.QU60km_polar.wolf +++ b/configs/lanl/config.20170207.MPAS-SeaIce.QU60km_polar.wolf @@ -11,6 +11,13 @@ mainRunName = MPAS-SeaIce.QU60km_polar # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [input] ## options related to reading in the results to be analyzed diff --git a/configs/lanl/config.MatchBoth_orig b/configs/lanl/config.MatchBoth_orig index 7ec34e297..2c04fc127 100644 --- a/configs/lanl/config.MatchBoth_orig +++ b/configs/lanl/config.MatchBoth_orig @@ -11,6 +11,13 @@ mainRunName = MatchBoth_orig # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [input] ## options related to reading in the results to be analyzed diff --git a/configs/olcf/config.20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison b/configs/olcf/config.20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison index b57a1eef4..c4b567970 100644 --- a/configs/olcf/config.20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison +++ b/configs/olcf/config.20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison @@ -11,6 +11,13 @@ mainRunName = 20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [input] ## options related to reading in the results to be analyzed diff --git a/configs/olcf/config.20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison b/configs/olcf/config.20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison index 49c4b749c..ceb34c0e2 100644 --- a/configs/olcf/config.20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison +++ b/configs/olcf/config.20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison @@ -11,6 +11,13 @@ mainRunName = 20170915.beta2.A_WCYCL1850S.ne30_oECv3_ICG.edison # MPAS-Analysis) preprocessedReferenceRunName = B1850C5_ne30_v0.4 +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [input] ## options related to reading in the results to be analyzed diff --git a/configs/olcf/config.GMPAS-IAF_oRRS18to6v3.titan b/configs/olcf/config.GMPAS-IAF_oRRS18to6v3.titan index 582c94e86..b38ae1983 100644 --- a/configs/olcf/config.GMPAS-IAF_oRRS18to6v3.titan +++ b/configs/olcf/config.GMPAS-IAF_oRRS18to6v3.titan @@ -4,10 +4,7 @@ # mainRunName is a name that identifies the simulation being analyzed. mainRunName = GMPAS-IAF_oRRS18to6v3 -# referenceRunName is the name of a reference run to compare against (or None -# to turn off comparison with a reference, e.g. if no reference case is -# available) -referenceRunName = None + # preprocessedReferenceRunName is the name of a reference run that has been # preprocessed to compare against (or None to turn off comparison). Reference # runs of this type would have preprocessed results because they were not @@ -15,6 +12,13 @@ referenceRunName = None # MPAS-Analysis) preprocessedReferenceRunName = None +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [input] ## options related to reading in the results to be analyzed diff --git a/configs/theta/config.20171031.tenYearTest.GMPAS-IAF.T62_oEC60to30v3wLI.60layer.theta b/configs/theta/config.20171031.tenYearTest.GMPAS-IAF.T62_oEC60to30v3wLI.60layer.theta index 328e09307..2d87a8b84 100644 --- a/configs/theta/config.20171031.tenYearTest.GMPAS-IAF.T62_oEC60to30v3wLI.60layer.theta +++ b/configs/theta/config.20171031.tenYearTest.GMPAS-IAF.T62_oEC60to30v3wLI.60layer.theta @@ -11,6 +11,13 @@ mainRunName = 20171031.tenYearTest.GMPAS-IAF.T62_oEC60to30v3wLI.60layer.theta # MPAS-Analysis) preprocessedReferenceRunName = None +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [execute] ## options related to executing parallel tasks diff --git a/docs/api.rst b/docs/api.rst index f5e9bbe01..5a240eff0 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -59,6 +59,7 @@ Ocean tasks TimeSeriesTemperatureAnomaly TimeSeriesSalinityAnomaly TimeSeriesSST + TimeSeriesAntarcticMelt .. currentmodule:: mpas_analysis.ocean.compute_anomaly_subtask @@ -166,13 +167,18 @@ Climatology compute_monthly_climatology compute_climatology add_years_months_days_in_month + get_unmasked_mpas_climatology_directory + get_unmasked_mpas_climatology_file_name + get_masked_mpas_climatology_file_name + get_remapped_mpas_climatology_file_name MpasClimatologyTask MpasClimatologyTask.add_variables MpasClimatologyTask.get_file_name RemapMpasClimatologySubtask - RemapMpasClimatologySubtask.get_file_name + RemapMpasClimatologySubtask.get_masked_file_name + RemapMpasClimatologySubtask.get_remapped_file_name RemapObservedClimatologySubtask RemapObservedClimatologySubtask.get_observation_descriptor @@ -187,7 +193,10 @@ Time Series :toctree: generated/ cache_time_series + compute_moving_avg_anomaly_from_start + compute_moving_avg + MpasTimeSeriesTask Interpolation ------------- diff --git a/mpas_analysis/config.default b/mpas_analysis/config.default index e751ebc1c..a96053572 100644 --- a/mpas_analysis/config.default +++ b/mpas_analysis/config.default @@ -22,10 +22,7 @@ # mainRunName is a name that identifies the simulation being analyzed. mainRunName = runName -# referenceRunName is the name of a reference run to compare against (or None -# to turn off comparison with a reference, e.g. if no reference case is -# available) -referenceRunName = None + # preprocessedReferenceRunName is the name of a reference run that has been # preprocessed to compare against (or None to turn off comparison). Reference # runs of this type would have preprocessed results because they were not @@ -33,6 +30,13 @@ referenceRunName = None # MPAS-Analysis) preprocessedReferenceRunName = None +# config file for a reference run to which this run will be compared. The +# analysis should have already been run to completion once with this config +# file, so that the relevant MPAS climatologies already exist and have been +# remapped to the comparison grid. Leave this option commented out if no +# reference run is desired. +# referenceRunConfigFile = /path/to/config/file + [execute] ## options related to executing parallel tasks @@ -73,7 +77,7 @@ oceanStreamsFileName = streams.ocean seaIceNamelistFileName = mpas-cice_in seaIceStreamsFileName = streams.cice -# names of ocean and sea ice meshes (e.g. EC60to30, QU240, RRS30to10, etc.) +# names of ocean and sea ice meshes (e.g. oEC60to30, oQU240, oRRS30to10, etc.) mpasMeshName = mesh # The system has a limit to how many files can be open at one time. By @@ -96,6 +100,7 @@ maxChunkSize = 10000 # placed in the output mappingSubdirectory # mappingDirectory = /dir/for/mapping/files + [output] ## options related to writing out plots, intermediate cached data sets, logs, ## etc. @@ -109,7 +114,6 @@ scratchSubdirectory = scratch plotsSubdirectory = plots logsSubdirectory = logs mpasClimatologySubdirectory = clim/mpas -mpasRemappedClimSubdirectory = clim/mpas/remapped mappingSubdirectory = mapping timeSeriesSubdirectory = timeseries timeCacheSubdirectory = timecache diff --git a/mpas_analysis/ocean/climatology_map_antarctic_melt.py b/mpas_analysis/ocean/climatology_map_antarctic_melt.py index 37a0b2f8b..47e4a82af 100644 --- a/mpas_analysis/ocean/climatology_map_antarctic_melt.py +++ b/mpas_analysis/ocean/climatology_map_antarctic_melt.py @@ -28,18 +28,22 @@ class ClimatologyMapAntarcticMelt(AnalysisTask): # {{{ ------- Xylar Asay-Davis """ - def __init__(self, config, mpasClimatologyTask): # {{{ + def __init__(self, config, mpasClimatologyTask, + refConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -57,17 +61,6 @@ def __init__(self, config, mpasClimatologyTask): # {{{ mpasFieldName = 'timeMonthly_avg_landIceFreshwaterFlux' iselValues = None - observationTitleLabel = \ - 'Observations (Rignot et al, 2013)' - - observationsDirectory = build_config_full_path( - config, 'oceanObservations', 'meltSubdirectory') - - obsFileName = \ - '{}/Rignot_2013_melt_rates_6000.0x6000.0km_10.0km_' \ - 'Antarctic_stereo.nc'.format(observationsDirectory) - obsFieldName = 'meltRate' - # read in what seasons we want to plot seasons = config.getExpression(sectionName, 'seasons') @@ -93,31 +86,60 @@ def __init__(self, config, mpasClimatologyTask): # {{{ seasons=seasons, iselValues=iselValues) - remapObservationsSubtask = RemapObservedAntarcticMeltClimatology( - parentTask=self, seasons=seasons, fileName=obsFileName, - outFilePrefix=obsFieldName, - comparisonGridNames=comparisonGridNames) - self.add_subtask(remapObservationsSubtask) + if refConfig is None: + + refTitleLabel = \ + 'Observations (Rignot et al, 2013)' + + observationsDirectory = build_config_full_path( + config, 'oceanObservations', 'meltSubdirectory') + + obsFileName = \ + '{}/Rignot_2013_melt_rates_6000.0x6000.0km_10.0km_' \ + 'Antarctic_stereo.nc'.format(observationsDirectory) + refFieldName = 'meltRate' + outFileLabel = 'meltRignot' + galleryName = 'Observations: Rignot et al. (2013)' + + remapObservationsSubtask = RemapObservedAntarcticMeltClimatology( + parentTask=self, seasons=seasons, fileName=obsFileName, + outFilePrefix=refFieldName, + comparisonGridNames=comparisonGridNames) + self.add_subtask(remapObservationsSubtask) + diffTitleLabel = 'Model - Observations' + + else: + remapObservationsSubtask = None + refRunName = refConfig.get('runs', 'mainRunName') + galleryName = None + refTitleLabel = 'Ref: {}'.format(refRunName) + + refFieldName = mpasFieldName + outFileLabel = 'melt' + diffTitleLabel = 'Main - Reference' + for comparisonGridName in comparisonGridNames: for season in seasons: # make a new subtask for this season and comparison grid subtask = PlotClimatologyMapSubtask(self, season, comparisonGridName, remapClimatologySubtask, - remapObservationsSubtask) + remapObservationsSubtask, + refConfig) subtask.set_plot_info( - outFileLabel='meltRignot', + outFileLabel=outFileLabel, fieldNameInTitle='Melt Rate', mpasFieldName=mpasFieldName, - obsFieldName=obsFieldName, - observationTitleLabel=observationTitleLabel, + refFieldName=refFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel=diffTitleLabel, unitsLabel=r'm a$^{-1}$', imageCaption='Antarctic Melt Rate', galleryGroup='Melt Rate', groupSubtitle=None, groupLink='antarctic_melt', - galleryName='Observations: Rignot et al. (2013)') + galleryName=galleryName) self.add_subtask(subtask) # }}} diff --git a/mpas_analysis/ocean/climatology_map_mld.py b/mpas_analysis/ocean/climatology_map_mld.py index 762a665e4..174f065a6 100644 --- a/mpas_analysis/ocean/climatology_map_mld.py +++ b/mpas_analysis/ocean/climatology_map_mld.py @@ -27,18 +27,22 @@ class ClimatologyMapMLD(AnalysisTask): # {{{ ------- Luke Van Roekel, Xylar Asay-Davis, Milena Veneziani """ - def __init__(self, config, mpasClimatologyTask): # {{{ + def __init__(self, config, mpasClimatologyTask, + refConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -55,17 +59,6 @@ def __init__(self, config, mpasClimatologyTask): # {{{ mpasFieldName = 'timeMonthly_avg_dThreshMLD' iselValues = None - observationTitleLabel = \ - 'Observations (HolteTalley density threshold MLD)' - - observationsDirectory = build_config_full_path( - config, 'oceanObservations', - '{}Subdirectory'.format(fieldName)) - - obsFileName = \ - "{}/holtetalley_mld_climatology.nc".format(observationsDirectory) - obsFieldName = 'mld' - # read in what seasons we want to plot seasons = config.getExpression(sectionName, 'seasons') @@ -91,31 +84,60 @@ def __init__(self, config, mpasClimatologyTask): # {{{ seasons=seasons, iselValues=iselValues) - remapObservationsSubtask = RemapObservedMLDClimatology( - parentTask=self, seasons=seasons, fileName=obsFileName, - outFilePrefix=obsFieldName, - comparisonGridNames=comparisonGridNames) - self.add_subtask(remapObservationsSubtask) + if refConfig is None: + + observationsDirectory = build_config_full_path( + config, 'oceanObservations', + '{}Subdirectory'.format(fieldName)) + + obsFileName = "{}/holtetalley_mld_climatology.nc".format( + observationsDirectory) + + refFieldName = 'mld' + outFileLabel = 'mldHolteTalleyARGO' + + remapObservationsSubtask = RemapObservedMLDClimatology( + parentTask=self, seasons=seasons, fileName=obsFileName, + outFilePrefix=refFieldName, + comparisonGridNames=comparisonGridNames) + self.add_subtask(remapObservationsSubtask) + galleryName = 'Observations: Holte-Talley ARGO' + refTitleLabel = \ + 'Observations (HolteTalley density threshold MLD)' + diffTitleLabel = 'Model - Observations' + + else: + remapObservationsSubtask = None + refRunName = refConfig.get('runs', 'mainRunName') + galleryName = None + refTitleLabel = 'Ref: {}'.format(refRunName) + + refFieldName = mpasFieldName + outFileLabel = 'mld' + diffTitleLabel = 'Main - Reference' + for comparisonGridName in comparisonGridNames: for season in seasons: # make a new subtask for this season and comparison grid subtask = PlotClimatologyMapSubtask(self, season, comparisonGridName, remapClimatologySubtask, - remapObservationsSubtask) + remapObservationsSubtask, + refConfig) subtask.set_plot_info( - outFileLabel='mldHolteTalleyARGO', + outFileLabel=outFileLabel, fieldNameInTitle='MLD', mpasFieldName=mpasFieldName, - obsFieldName=obsFieldName, - observationTitleLabel=observationTitleLabel, + refFieldName=refFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel=diffTitleLabel, unitsLabel=r'm', imageCaption='Mean Mixed-Layer Depth', galleryGroup='Mixed-Layer Depth', groupSubtitle=None, groupLink='mld', - galleryName='Observations: Holte-Talley ARGO') + galleryName=galleryName) self.add_subtask(subtask) # }}} diff --git a/mpas_analysis/ocean/climatology_map_sose.py b/mpas_analysis/ocean/climatology_map_sose.py index 5fa894443..894b8ebd6 100644 --- a/mpas_analysis/ocean/climatology_map_sose.py +++ b/mpas_analysis/ocean/climatology_map_sose.py @@ -35,18 +35,22 @@ class ClimatologyMapSoseTemperature(AnalysisTask): # {{{ Xylar Asay-Davis """ - def __init__(self, config, mpasClimatologyTask): # {{{ + def __init__(self, config, mpasClimatologyTask, + refConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -63,16 +67,6 @@ def __init__(self, config, mpasClimatologyTask): # {{{ mpasFieldName = 'timeMonthly_avg_activeTracers_temperature' iselValues = None - observationTitleLabel = 'State Estimate (SOSE)' - - observationsDirectory = build_config_full_path( - config, 'oceanObservations', 'soseSubdirectory') - - obsFileName = \ - '{}/SOSE_2005-2010_monthly_pot_temp_6000.0x' \ - '6000.0km_10.0km_Antarctic_stereo.nc'.format(observationsDirectory) - obsFieldName = 'theta' - # read in what seasons we want to plot seasons = config.getExpression(sectionName, 'seasons') @@ -105,14 +99,42 @@ def __init__(self, config, mpasClimatologyTask): # {{{ comparisonGridNames=comparisonGridNames, iselValues=iselValues) - remapObservationsSubtask = RemapSoseClimatology( - parentTask=self, seasons=seasons, fileName=obsFileName, - outFilePrefix=obsFieldName, - fieldName=obsFieldName, - botFieldName='botTheta', - depths=depths, - comparisonGridNames=comparisonGridNames) - self.add_subtask(remapObservationsSubtask) + if refConfig is None: + + refTitleLabel = 'State Estimate (SOSE)' + + observationsDirectory = build_config_full_path( + config, 'oceanObservations', 'soseSubdirectory') + + obsFileName = \ + '{}/SOSE_2005-2010_monthly_pot_temp_6000.0x' \ + '6000.0km_10.0km_Antarctic_stereo.nc'.format( + observationsDirectory) + refFieldName = 'theta' + outFileLabel = 'tempSOSE' + galleryName = 'State Estimate: SOSE' + diffTitleLabel = 'Model - State Estimate' + + remapObservationsSubtask = RemapSoseClimatology( + parentTask=self, seasons=seasons, fileName=obsFileName, + outFilePrefix=refFieldName, + fieldName=refFieldName, + botFieldName='botTheta', + depths=depths, + comparisonGridNames=comparisonGridNames) + + self.add_subtask(remapObservationsSubtask) + + else: + remapObservationsSubtask = None + refRunName = refConfig.get('runs', 'mainRunName') + galleryName = 'Ref: {}'.format(refRunName) + refTitleLabel = galleryName + + refFieldName = mpasFieldName + outFileLabel = 'temp' + diffTitleLabel = 'Main - Reference' + for comparisonGridName in comparisonGridNames: for season in seasons: for depth in depths: @@ -122,21 +144,22 @@ def __init__(self, config, mpasClimatologyTask): # {{{ comparisonGridName=comparisonGridName, remapMpasClimatologySubtask=remapClimatologySubtask, remapObsClimatologySubtask=remapObservationsSubtask, + refConfig=refConfig, depth=depth) subtask.set_plot_info( - outFileLabel='tempSOSE', + outFileLabel=outFileLabel, fieldNameInTitle='Temperature', mpasFieldName=mpasFieldName, - obsFieldName=obsFieldName, - observationTitleLabel=observationTitleLabel, - diffTitleLabel='Model - State Estimate', + refFieldName=refFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel=diffTitleLabel, unitsLabel=r'$^\circ$C', imageCaption='Temperature', galleryGroup='Temperature', groupSubtitle=None, groupLink='temp', - galleryName='State Estimate: SOSE') + galleryName=galleryName) self.add_subtask(subtask) # }}} @@ -154,18 +177,22 @@ class ClimatologyMapSoseSalinity(AnalysisTask): # {{{ Xylar Asay-Davis """ - def __init__(self, config, mpasClimatologyTask): # {{{ + def __init__(self, config, mpasClimatologyTask, + refConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -182,16 +209,6 @@ def __init__(self, config, mpasClimatologyTask): # {{{ mpasFieldName = 'timeMonthly_avg_activeTracers_salinity' iselValues = None - observationTitleLabel = 'State Estimate (SOSE)' - - observationsDirectory = build_config_full_path( - config, 'oceanObservations', 'soseSubdirectory') - - obsFileName = \ - '{}/SOSE_2005-2010_monthly_salinity_6000.0x' \ - '6000.0km_10.0km_Antarctic_stereo.nc'.format(observationsDirectory) - obsFieldName = 'salinity' - # read in what seasons we want to plot seasons = config.getExpression(sectionName, 'seasons') @@ -224,14 +241,42 @@ def __init__(self, config, mpasClimatologyTask): # {{{ comparisonGridNames=comparisonGridNames, iselValues=iselValues) - remapObservationsSubtask = RemapSoseClimatology( - parentTask=self, seasons=seasons, fileName=obsFileName, - outFilePrefix=obsFieldName, - fieldName=obsFieldName, - botFieldName='botSalinity', - depths=depths, - comparisonGridNames=comparisonGridNames) - self.add_subtask(remapObservationsSubtask) + if refConfig is None: + + refTitleLabel = 'State Estimate (SOSE)' + + observationsDirectory = build_config_full_path( + config, 'oceanObservations', 'soseSubdirectory') + + obsFileName = \ + '{}/SOSE_2005-2010_monthly_salinity_6000.0x' \ + '6000.0km_10.0km_Antarctic_stereo.nc'.format( + observationsDirectory) + refFieldName = 'salinity' + outFileLabel = 'salinSOSE' + galleryName = 'State Estimate: SOSE' + diffTitleLabel = 'Model - State Estimate' + + remapObservationsSubtask = RemapSoseClimatology( + parentTask=self, seasons=seasons, fileName=obsFileName, + outFilePrefix=refFieldName, + fieldName=refFieldName, + botFieldName='botSalinity', + depths=depths, + comparisonGridNames=comparisonGridNames) + + self.add_subtask(remapObservationsSubtask) + + else: + remapObservationsSubtask = None + refRunName = refConfig.get('runs', 'mainRunName') + galleryName = None + refTitleLabel = 'Ref: {}'.format(refRunName) + + refFieldName = mpasFieldName + outFileLabel = 'salin' + diffTitleLabel = 'Main - Reference' + for comparisonGridName in comparisonGridNames: for season in seasons: for depth in depths: @@ -241,21 +286,22 @@ def __init__(self, config, mpasClimatologyTask): # {{{ comparisonGridName=comparisonGridName, remapMpasClimatologySubtask=remapClimatologySubtask, remapObsClimatologySubtask=remapObservationsSubtask, + refConfig=refConfig, depth=depth) subtask.set_plot_info( - outFileLabel='salinSOSE', + outFileLabel=outFileLabel, fieldNameInTitle='Salinity', mpasFieldName=mpasFieldName, - obsFieldName=obsFieldName, - observationTitleLabel=observationTitleLabel, - diffTitleLabel='Model - State Estimate', + refFieldName=refFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel=diffTitleLabel, unitsLabel=r'PSU', imageCaption='Salinity', galleryGroup='Salinity', groupSubtitle=None, groupLink='salin', - galleryName='State Estimate: SOSE') + galleryName=galleryName) self.add_subtask(subtask) # }}} @@ -325,7 +371,8 @@ def __init__(self, parentTask, seasons, fileName, outFilePrefix, self.botFieldName = botFieldName self.depths = depths - # call the constructor from the base class (AnalysisTask) + # call the constructor from the base class + # (RemapObservedClimatologySubtask) super(RemapSoseClimatology, self).__init__( parentTask, seasons, fileName, outFilePrefix, comparisonGridNames, subtaskName) diff --git a/mpas_analysis/ocean/climatology_map_sss.py b/mpas_analysis/ocean/climatology_map_sss.py index 168a3d59c..0fdc0877d 100644 --- a/mpas_analysis/ocean/climatology_map_sss.py +++ b/mpas_analysis/ocean/climatology_map_sss.py @@ -25,46 +25,38 @@ class ClimatologyMapSSS(AnalysisTask): # {{{ ------- Luke Van Roekel, Xylar Asay-Davis, Milena Veneziani """ - def __init__(self, config, mpasClimatologyTask): # {{{ + def __init__(self, config, mpasClimatologyTask, + refConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis """ - self.fieldName = 'sss' + fieldName = 'sss' # call the constructor from the base class (AnalysisTask) super(ClimatologyMapSSS, self).__init__( config=config, taskName='climatologyMapSSS', componentName='ocean', - tags=['climatology', 'horizontalMap', self.fieldName]) + tags=['climatology', 'horizontalMap', fieldName]) mpasFieldName = 'timeMonthly_avg_activeTracers_salinity' iselValues = {'nVertLevels': 0} sectionName = self.taskName - observationTitleLabel = \ - 'Observations (Aquarius, 2011-2014)' - - observationsDirectory = build_config_full_path( - config, 'oceanObservations', - '{}Subdirectory'.format(self.fieldName)) - - obsFileName = \ - "{}/Aquarius_V3_SSS_Monthly.nc".format( - observationsDirectory) - obsFieldName = 'sss' - # read in what seasons we want to plot seasons = config.getExpression(sectionName, 'seasons') @@ -84,37 +76,67 @@ def __init__(self, config, mpasClimatologyTask): # {{{ remapClimatologySubtask = RemapMpasClimatologySubtask( mpasClimatologyTask=mpasClimatologyTask, parentTask=self, - climatologyName=self.fieldName, + climatologyName=fieldName, variableList=[mpasFieldName], comparisonGridNames=comparisonGridNames, seasons=seasons, iselValues=iselValues) - remapObservationsSubtask = RemapObservedSSSClimatology( - parentTask=self, seasons=seasons, fileName=obsFileName, - outFilePrefix=obsFieldName, - comparisonGridNames=comparisonGridNames) - self.add_subtask(remapObservationsSubtask) + if refConfig is None: + + refTitleLabel = \ + 'Observations (Aquarius, 2011-2014)' + + observationsDirectory = build_config_full_path( + config, 'oceanObservations', + '{}Subdirectory'.format(fieldName)) + + obsFileName = \ + "{}/Aquarius_V3_SSS_Monthly.nc".format( + observationsDirectory) + refFieldName = 'sss' + outFileLabel = 'sssAquarius' + galleryName = 'Observations: Aquarius' + + remapObservationsSubtask = RemapObservedSSSClimatology( + parentTask=self, seasons=seasons, fileName=obsFileName, + outFilePrefix=refFieldName, + comparisonGridNames=comparisonGridNames) + self.add_subtask(remapObservationsSubtask) + diffTitleLabel = 'Model - Observations' + + else: + remapObservationsSubtask = None + refRunName = refConfig.get('runs', 'mainRunName') + galleryName = None + refTitleLabel = 'Ref: {}'.format(refRunName) + + refFieldName = mpasFieldName + outFileLabel = 'sss' + diffTitleLabel = 'Main - Reference' + for comparisonGridName in comparisonGridNames: for season in seasons: # make a new subtask for this season and comparison grid subtask = PlotClimatologyMapSubtask(self, season, comparisonGridName, remapClimatologySubtask, - remapObservationsSubtask) + remapObservationsSubtask, + refConfig) subtask.set_plot_info( - outFileLabel='sssAquarius', + outFileLabel=outFileLabel, fieldNameInTitle='SSS', mpasFieldName=mpasFieldName, - obsFieldName=obsFieldName, - observationTitleLabel=observationTitleLabel, + refFieldName=refFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel=diffTitleLabel, unitsLabel=r'PSU', imageCaption='Mean Sea Surface Salinity', galleryGroup='Sea Surface Salinity', groupSubtitle=None, groupLink='sss', - galleryName='Observations: Aquarius') + galleryName=galleryName) self.add_subtask(subtask) # }}} diff --git a/mpas_analysis/ocean/climatology_map_sst.py b/mpas_analysis/ocean/climatology_map_sst.py index db4bcca5a..0a5b1d699 100644 --- a/mpas_analysis/ocean/climatology_map_sst.py +++ b/mpas_analysis/ocean/climatology_map_sst.py @@ -25,28 +25,32 @@ class ClimatologyMapSST(AnalysisTask): # {{{ ------- Luke Van Roekel, Xylar Asay-Davis, Milena Veneziani """ - def __init__(self, config, mpasClimatologyTask): # {{{ + def __init__(self, config, mpasClimatologyTask, + refConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis """ - self.fieldName = 'sst' + fieldName = 'sst' # call the constructor from the base class (AnalysisTask) super(ClimatologyMapSST, self).__init__( config=config, taskName='climatologyMapSST', componentName='ocean', - tags=['climatology', 'horizontalMap', self.fieldName]) + tags=['climatology', 'horizontalMap', fieldName]) mpasFieldName = 'timeMonthly_avg_activeTracers_temperature' iselValues = {'nVertLevels': 0} @@ -58,25 +62,6 @@ def __init__(self, config, mpasClimatologyTask): # {{{ climEndYear = config.getint('oceanObservations', 'sstClimatologyEndYear') - if climStartYear < 1925: - period = 'pre-industrial' - else: - period = 'present-day' - - observationTitleLabel = \ - 'Observations (Hadley/OI, {} {:04d}-{:04d})'.format(period, - climStartYear, - climEndYear) - - observationsDirectory = build_config_full_path( - config, 'oceanObservations', - '{}Subdirectory'.format(self.fieldName)) - - obsFileName = \ - "{}/MODEL.SST.HAD187001-198110.OI198111-201203.nc".format( - observationsDirectory) - obsFieldName = 'sst' - # read in what seasons we want to plot seasons = config.getExpression(sectionName, 'seasons') @@ -91,42 +76,77 @@ def __init__(self, config, mpasClimatologyTask): # {{{ raise ValueError('config section {} does not contain valid list ' 'of comparison grids'.format(sectionName)) - # the variable self.mpasFieldName will be added to mpasClimatologyTask + # the variable mpasFieldName will be added to mpasClimatologyTask # along with the seasons. remapClimatologySubtask = RemapMpasClimatologySubtask( mpasClimatologyTask=mpasClimatologyTask, parentTask=self, - climatologyName=self.fieldName, + climatologyName=fieldName, variableList=[mpasFieldName], comparisonGridNames=comparisonGridNames, seasons=seasons, iselValues=iselValues) - remapObservationsSubtask = RemapObservedSSTClimatology( - parentTask=self, seasons=seasons, fileName=obsFileName, - outFilePrefix=obsFieldName, - comparisonGridNames=comparisonGridNames) - self.add_subtask(remapObservationsSubtask) + if refConfig is None: + if climStartYear < 1925: + period = 'pre-industrial' + else: + period = 'present-day' + + refTitleLabel = \ + 'Observations (Hadley/OI, {} {:04d}-{:04d})'.format( + period, climStartYear, climEndYear) + + observationsDirectory = build_config_full_path( + config, 'oceanObservations', + '{}Subdirectory'.format(fieldName)) + + obsFileName = \ + "{}/MODEL.SST.HAD187001-198110.OI198111-201203.nc".format( + observationsDirectory) + refFieldName = 'sst' + outFileLabel = 'sstHADOI' + galleryName = 'Observations: Hadley-NOAA-OI' + + remapObservationsSubtask = RemapObservedSSTClimatology( + parentTask=self, seasons=seasons, fileName=obsFileName, + outFilePrefix=refFieldName, + comparisonGridNames=comparisonGridNames) + self.add_subtask(remapObservationsSubtask) + diffTitleLabel = 'Model - Observations' + + else: + remapObservationsSubtask = None + refRunName = refConfig.get('runs', 'mainRunName') + galleryName = None + refTitleLabel = 'Ref: {}'.format(refRunName) + + refFieldName = mpasFieldName + outFileLabel = 'sst' + diffTitleLabel = 'Main - Reference' + for comparisonGridName in comparisonGridNames: for season in seasons: # make a new subtask for this season and comparison grid subtask = PlotClimatologyMapSubtask(self, season, comparisonGridName, remapClimatologySubtask, - remapObservationsSubtask) + remapObservationsSubtask, + refConfig) subtask.set_plot_info( - outFileLabel='sstHADOI', + outFileLabel=outFileLabel, fieldNameInTitle='SST', mpasFieldName=mpasFieldName, - obsFieldName=obsFieldName, - observationTitleLabel=observationTitleLabel, + refFieldName=refFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel=diffTitleLabel, unitsLabel=r'$^o$C', imageCaption='Mean Sea Surface Temperature', galleryGroup='Sea Surface Temperature', groupSubtitle=None, groupLink='sst', - galleryName='Observations: Hadley-NOAA-OI') + galleryName=galleryName) self.add_subtask(subtask) # }}} diff --git a/mpas_analysis/ocean/index_nino34.py b/mpas_analysis/ocean/index_nino34.py index 5b8fec209..608621328 100644 --- a/mpas_analysis/ocean/index_nino34.py +++ b/mpas_analysis/ocean/index_nino34.py @@ -36,23 +36,30 @@ class IndexNino34(AnalysisTask): # {{{ mpasTimeSeriesTask : ``MpasTimeSeriesTask`` The task that extracts the time series from MPAS monthly output + refConfig : ``MpasAnalysisConfigParser`` + Configuration options for a reference run (if any) + Authors ------- Luke Van Roekel, Xylar Asay-Davis ''' - def __init__(self, config, mpasTimeSeriesTask): # {{{ + def __init__(self, config, mpasTimeSeriesTask, refConfig=None): + # {{{ ''' Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasTimeSeriesTask : ``MpasTimeSeriesTask`` The task that extracts the time series from MPAS monthly output + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -66,6 +73,7 @@ def __init__(self, config, mpasTimeSeriesTask): # {{{ tags=['timeSeries', 'index', 'nino']) self.mpasTimeSeriesTask = mpasTimeSeriesTask + self.refConfig = refConfig self.run_after(mpasTimeSeriesTask) @@ -87,9 +95,6 @@ def setup_and_check(self): # {{{ # self.calendar super(IndexNino34, self).setup_and_check() - self.startDate = self.config.get('index', 'startDate') - self.endDate = self.config.get('index', 'endDate') - self.variableList = \ ['timeMonthly_avg_avgValueWithinOceanRegion_avgSurfaceTemperature'] self.mpasTimeSeriesTask.add_variables(variableList=self.variableList) @@ -125,6 +130,12 @@ def run_task(self): # {{{ config = self.config calendar = self.calendar + startDate = self.config.get('index', 'startDate') + endDate = self.config.get('index', 'endDate') + + startYear = self.config.getint('index', 'startYear') + endYear = self.config.getint('index', 'endYear') + dataSource = config.get('indexNino34', 'observationData') observationsDirectory = build_config_full_path( @@ -153,8 +164,8 @@ def run_task(self): # {{{ ds = open_mpas_dataset(fileName=self.inputFile, calendar=calendar, variableList=self.variableList, - startDate=self.startDate, - endDate=self.endDate) + startDate=startDate, + endDate=endDate) # Observations have been processed to the nino34Index prior to reading dsObs = xr.open_dataset(dataPath, decode_cf=False, decode_times=False) @@ -180,8 +191,12 @@ def run_task(self): # {{{ # Compute the observational spectra over the last 30 years for # comparison. Only saving the spectra - subsetStartYear = 1976 subsetEndYear = 2016 + if self.refConfig is None: + subsetStartYear = 1976 + else: + # make the subset the same length as the input data set + subsetStartYear = subsetEndYear - (endYear - startYear) time_start = datetime_to_days(datetime.datetime(subsetStartYear, 1, 1), calendar=calendar) time_end = datetime_to_days(datetime.datetime(subsetEndYear, 12, 31), @@ -189,21 +204,50 @@ def run_task(self): # {{{ nino34Subset = nino34Obs.sel(Time=slice(time_start, time_end)) spectraSubset = self._compute_nino34_spectra(nino34Subset) + if self.refConfig is None: + nino34s = [nino34Obs[2:-3], nino34Subset, nino34Main[2:-3]] + titles = ['{} (Full Record)'.format(obsTitle), + '{} ({} - {})'.format(obsTitle, subsetStartYear, + subsetEndYear), + mainRunName] + spectra = [spectraObs, spectraSubset, spectraMain] + else: + baseDirectory = build_config_full_path( + self.refConfig, 'output', 'timeSeriesSubdirectory') + + refFileName = '{}/{}.nc'.format( + baseDirectory, self.mpasTimeSeriesTask.fullTaskName) + + dsRef = open_mpas_dataset( + fileName=refFileName, + calendar=calendar, + variableList=self.variableList) + + regionSSTRef = dsRef[varName].isel(nOceanRegions=regionIndex) + nino34Ref = self._compute_nino34_index(regionSSTRef, calendar) + + nino34s = [nino34Subset, nino34Main[2:-3], nino34Ref[2:-3]] + refRunName = self.refConfig.get('runs', 'mainRunName') + + spectraRef = self._compute_nino34_spectra(nino34Ref) + + titles = ['{} ({} - {})'.format(obsTitle, subsetStartYear, + subsetEndYear), + mainRunName, + 'Ref: {}'.format(refRunName)] + spectra = [spectraSubset, spectraMain, spectraRef] + # Convert frequencies to period in years - for spectra in [spectraMain, spectraObs, spectraSubset]: - spectra['period'] = \ - 1.0 / (constants.eps + spectra['f']*constants.sec_per_year) + for s in spectra: + s['period'] = \ + 1.0 / (constants.eps + s['f']*constants.sec_per_year) self.logger.info(' Plot NINO3.4 index and spectra...') outFileName = '{}/NINO34_{}.png'.format(self.plotsDirectory, mainRunName) - titles = ['{} (Full Record)'.format(obsTitle), - '{} ({} - {})'.format(obsTitle, subsetStartYear, - subsetEndYear), - mainRunName] self._nino34_timeseries_plot( - nino34s=[nino34Obs[2:-3], nino34Subset, nino34Main[2:-3]], + nino34s=nino34s, title='NINO 3.4 Index', panelTitles=titles, outFileName=outFileName) @@ -214,7 +258,7 @@ def run_task(self): # {{{ outFileName = '{}/NINO34_spectra_{}.png'.format(self.plotsDirectory, mainRunName) self._nino34_spectra_plot( - spectra=[spectraObs, spectraSubset, spectraMain], + spectra=spectra, title='NINO3.4 power spectrum', panelTitles=titles, outFileName=outFileName) diff --git a/mpas_analysis/ocean/meridional_heat_transport.py b/mpas_analysis/ocean/meridional_heat_transport.py index ef6d871da..e6b139e76 100644 --- a/mpas_analysis/ocean/meridional_heat_transport.py +++ b/mpas_analysis/ocean/meridional_heat_transport.py @@ -9,11 +9,14 @@ from ..shared.plot.plotting import plot_vertical_section,\ setup_colormap, plot_1D -from ..shared.io.utility import build_config_full_path +from ..shared.io.utility import build_config_full_path, make_directories +from ..shared.io import write_netcdf, subset_variables from ..shared import AnalysisTask from ..shared.html import write_image_xml +from ..shared.climatology import get_unmasked_mpas_climatology_file_name + class MeridionalHeatTransport(AnalysisTask): # {{{ ''' @@ -25,23 +28,29 @@ class MeridionalHeatTransport(AnalysisTask): # {{{ mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted + refConfig : ``MpasAnalysisConfigParser`` + Configuration options for a reference run (if any) + Authors ------- Mark Petersen, Milena Veneziani, Xylar Asay-Davis ''' - def __init__(self, config, mpasClimatologyTask): # {{{ + def __init__(self, config, mpasClimatologyTask, refConfig=None): # {{{ ''' Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -57,16 +66,14 @@ def __init__(self, config, mpasClimatologyTask): # {{{ self.mpasClimatologyTask = mpasClimatologyTask self.run_after(mpasClimatologyTask) + self.refConfig = refConfig + # }}} def setup_and_check(self): # {{{ ''' Perform steps to set up the analysis and check for errors in the setup. - Raises - ------ - ValueError: if myArg has an invalid value - Authors ------- Mark Petersen, Milena Veneziani, Xylar Asay-Davis @@ -148,77 +155,106 @@ def run_task(self): # {{{ config = self.config + mainRunName = config.get('runs', 'mainRunName') + + depthLimGlobal = config.getExpression(self.sectionName, + 'depthLimGlobal') + xLimGlobal = config.getExpression(self.sectionName, 'xLimGlobal') movingAveragePoints = config.getint('meridionalHeatTransport', 'movingAveragePoints') - # Read in depth and MHT latitude points - # Latitude is from binBoundaryMerHeatTrans - try: - restartFileName = self.runStreams.readpath('restart')[0] - except ValueError: - raise IOError('No MPAS-O restart file found: need at least one ' - 'for MHT calcuation') - - with xr.open_dataset(restartFileName) as dsRestart: - refBottomDepth = dsRestart.refBottomDepth.values - - nVertLevels = len(refBottomDepth) - refLayerThickness = np.zeros(nVertLevels) - refLayerThickness[0] = refBottomDepth[0] - refLayerThickness[1:nVertLevels] = (refBottomDepth[1:nVertLevels] - - refBottomDepth[0:nVertLevels-1]) - - refZMid = -refBottomDepth + 0.5*refLayerThickness - - binBoundaryMerHeatTrans = None - # first try timeSeriesStatsMonthly for bin boundaries, then try - # meridionalHeatTranspor steram as a backup option - for streamName in ['timeSeriesStatsMonthlyOutput', - 'meridionalHeatTransportOutput']: - try: - inputFile = self.historyStreams.readpath(streamName)[0] - except ValueError: - raise IOError('At least one file from stream {} is needed to ' - 'compute MHT'.format(streamName)) - - with xr.open_dataset(inputFile) as ds: - if 'binBoundaryMerHeatTrans' in ds.data_vars: - binBoundaryMerHeatTrans = ds.binBoundaryMerHeatTrans.values - break + outputDirectory = build_config_full_path(config, 'output', + 'mpasClimatologySubdirectory') - if binBoundaryMerHeatTrans is None: - raise ValueError('Could not find binBoundaryMerHeatTrans in either' - 'timeSeriesStatsMonthlyOutput or ' - 'meridionalHeatTransportOutput streams') + make_directories(outputDirectory) - binBoundaryMerHeatTrans = np.rad2deg(binBoundaryMerHeatTrans) + outFileName = \ + '{}/meridionalHeatTransport_years{:04d}-{:04d}.nc'.format( + outputDirectory, self.startYear, self.endYear) - ###################################################################### - # Mark P Note: Currently only supports global MHT. - # Need to add variables merHeatTransLatRegion and - # merHeatTransLatZRegion - # These are not computed by default in ACME right now. - # Then we will need to add another section for regions with a loop - # over number of regions. - ###################################################################### - self.logger.info('\n Plotting global meridional heat transport') - - self.logger.info(' Load data...') - - climatologyFileName = self.mpasClimatologyTask.get_file_name( - season='ANN') + if os.path.exists(outFileName): + self.logger.info(' Reading results from previous analysis run...') + annualClimatology = xr.open_dataset(outFileName) + refZMid = annualClimatology.refZMid.values + binBoundaryMerHeatTrans = \ + annualClimatology.binBoundaryMerHeatTrans.values + else: - annualClimatology = xr.open_dataset(climatologyFileName) - annualClimatology = annualClimatology.isel(Time=0) + # Read in depth and MHT latitude points + # Latitude is from binBoundaryMerHeatTrans + try: + restartFileName = self.runStreams.readpath('restart')[0] + except ValueError: + raise IOError('No MPAS-O restart file found: need at least ' + 'one for MHT calcuation') + + with xr.open_dataset(restartFileName) as dsRestart: + refBottomDepth = dsRestart.refBottomDepth.values + + nVertLevels = len(refBottomDepth) + refLayerThickness = np.zeros(nVertLevels) + refLayerThickness[0] = refBottomDepth[0] + refLayerThickness[1:nVertLevels] = \ + refBottomDepth[1:nVertLevels] - refBottomDepth[0:nVertLevels-1] + + refZMid = -refBottomDepth + 0.5*refLayerThickness + + binBoundaryMerHeatTrans = None + # first try timeSeriesStatsMonthly for bin boundaries, then try + # meridionalHeatTranspor steram as a backup option + for streamName in ['timeSeriesStatsMonthlyOutput', + 'meridionalHeatTransportOutput']: + try: + inputFile = self.historyStreams.readpath(streamName)[0] + except ValueError: + raise IOError('At least one file from stream {} is needed ' + 'to compute MHT'.format(streamName)) + + with xr.open_dataset(inputFile) as ds: + if 'binBoundaryMerHeatTrans' in ds.data_vars: + binBoundaryMerHeatTrans = \ + ds.binBoundaryMerHeatTrans.values + break + + if binBoundaryMerHeatTrans is None: + raise ValueError('Could not find binBoundaryMerHeatTrans in ' + 'either timeSeriesStatsMonthlyOutput or ' + 'meridionalHeatTransportOutput streams') + + binBoundaryMerHeatTrans = np.rad2deg(binBoundaryMerHeatTrans) + + ################################################################### + # Mark P Note: Currently only supports global MHT. + # Need to add variables merHeatTransLatRegion and + # merHeatTransLatZRegion + # These are not computed by default in ACME right now. + # Then we will need to add another section for regions with a loop + # over number of regions. + ################################################################### + + self.logger.info('\n Plotting global meridional heat transport') + + self.logger.info(' Load data...') + + climatologyFileName = self.mpasClimatologyTask.get_file_name( + season='ANN') + + variableList = ['timeMonthly_avg_meridionalHeatTransportLat', + 'timeMonthly_avg_meridionalHeatTransportLatZ'] + + annualClimatology = xr.open_dataset(climatologyFileName) + annualClimatology = subset_variables(annualClimatology, + variableList) + annualClimatology = annualClimatology.isel(Time=0) + + annualClimatology.coords['refZMid'] = (('nVertLevels',), refZMid) + annualClimatology.coords['binBoundaryMerHeatTrans'] = \ + (('nMerHeatTransBinsP1',), binBoundaryMerHeatTrans) + + write_netcdf(annualClimatology, outFileName) # **** Plot MHT **** - # Define plotting variables - mainRunName = config.get('runs', 'mainRunName') - xLimGlobal = config.getExpression(self.sectionName, 'xLimGlobal') - depthLimGlobal = config.getExpression(self.sectionName, - 'depthLimGlobal') - self.logger.info(' Plot global MHT...') # Plot 1D MHT (zonally averaged, depth integrated) x = binBoundaryMerHeatTrans @@ -229,6 +265,12 @@ def run_task(self): # {{{ self.startYear, self.endYear, mainRunName) filePrefix = self.filePrefixes['mht'] figureName = '{}/{}.png'.format(self.plotsDirectory, filePrefix) + lineColors = ['r'] + lineWidths = [1.6] + legendText = [mainRunName] + xArrays = [x] + fieldArrays = [y] + errArrays = [None] if self.observationsFile is not None: # Load in observations dsObs = xr.open_dataset(self.observationsFile) @@ -238,23 +280,43 @@ def run_task(self): # {{{ ecmwfGlobal = dsObs.GLOBALECMWF_ADJUSTED ecmwfErrGlobal = dsObs.GLOBALECMWF_ERR - lineColors = ['r', 'b', 'g'] - lineWidths = [1.6, 1.2, 1.2] - legendText = ['model', 'NCEP', 'ECMWF'] - plot_1D(config, [x, xObs, xObs], - [y, ncepGlobal, ecmwfGlobal], - [None, ncepErrGlobal, ecmwfErrGlobal], - lineColors, lineWidths, legendText, - title, xLabel, yLabel, figureName, - xLim=xLimGlobal) - else: - lineColors = ['r'] - lineWidths = [1.6] + lineColors.extend(['b', 'g']) + lineWidths.extend([1.2, 1.2]) + legendText.extend(['NCEP', 'ECMWF']) + xArrays.extend([xObs, xObs]) + fieldArrays.extend([ncepGlobal, ecmwfGlobal]) + errArrays.extend([ncepErrGlobal, ecmwfErrGlobal]) + + if self.refConfig is not None: + + refStartYear = self.refConfig.getint('climatology', 'startYear') + refEndYear = self.refConfig.getint('climatology', 'endYear') + refDirectory = build_config_full_path( + self.refConfig, 'output', 'mpasClimatologySubdirectory') + + refFileName = \ + '{}/meridionalHeatTransport_years{:04d}-{:04d}.nc'.format( + refDirectory, refStartYear, refEndYear) + + dsRef = xr.open_dataset(refFileName) + refRunName = self.refConfig.get('runs', 'mainRunName') + + lineColors.append('k') + lineWidths.append(1.2) + legendText.append(refRunName) + xArrays.append(dsRef.binBoundaryMerHeatTrans) + fieldArrays.append( + dsRef.timeMonthly_avg_meridionalHeatTransportLat) + errArrays.append(None) + + if len(legendText) == 1: + # no need for a legend legendText = [None] - plot_1D(config, [x], [y], [None], - lineColors, lineWidths, legendText, - title, xLabel, yLabel, figureName, - xLim=xLimGlobal) + + plot_1D(config, xArrays, fieldArrays, errArrays, + lineColors, lineWidths, legendText, + title, xLabel, yLabel, figureName, + xLim=xLimGlobal) self._write_xml(filePrefix) diff --git a/mpas_analysis/ocean/plot_climatology_map_subtask.py b/mpas_analysis/ocean/plot_climatology_map_subtask.py index dd717f522..826046a0f 100644 --- a/mpas_analysis/ocean/plot_climatology_map_subtask.py +++ b/mpas_analysis/ocean/plot_climatology_map_subtask.py @@ -22,6 +22,8 @@ from ..shared.grid import interp_extrap_corner +from ..shared.climatology import get_remapped_mpas_climatology_file_name + def nans_to_numpy_mask(field): field = np.ma.masked_array( @@ -50,6 +52,10 @@ class PlotClimatologyMapSubtask(AnalysisTask): # {{{ The subtask for remapping the observational climatology that this subtask will plot + remapMpasRefClimatologySubtask : ``RemapMpasReferenceClimatologySubtask`` + The subtask for remapping the MPAS climatology for the reference + run that this subtask will plot + outFileLabel : str The prefix on each plot and associated XML file @@ -98,8 +104,8 @@ class PlotClimatologyMapSubtask(AnalysisTask): # {{{ """ def __init__(self, parentTask, season, comparisonGridName, - remapMpasClimatologySubtask, remapObsClimatologySubtask, - depth=None): + remapMpasClimatologySubtask, remapObsClimatologySubtask=None, + refConfig=None, depth=None): # {{{ ''' Construct one analysis subtask for each plot (i.e. each season and @@ -121,10 +127,13 @@ def __init__(self, parentTask, season, comparisonGridName, The subtask for remapping the MPAS climatology that this subtask will plot - remapObsClimatologySubtask : ``RemapObservedClimatologySubtask`` + remapObsClimatologySubtask : ``RemapObservedClimatologySubtask``, optional The subtask for remapping the observational climatology that this subtask will plot + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + depth : {float, 'top', 'bot'}, optional Depth the data is being plotted, 'top' for the sea surface 'bot' for the sea floor @@ -140,6 +149,7 @@ def __init__(self, parentTask, season, comparisonGridName, self.comparisonGridName = comparisonGridName self.remapMpasClimatologySubtask = remapMpasClimatologySubtask self.remapObsClimatologySubtask = remapObsClimatologySubtask + self.refConfig = refConfig subtaskName = 'plot{}_{}'.format(season, comparisonGridName) if depth is None: @@ -160,11 +170,12 @@ def __init__(self, parentTask, season, comparisonGridName, # this task should not run until the remapping subtasks are done, since # it relies on data from those subtasks self.run_after(remapMpasClimatologySubtask) - self.run_after(remapObsClimatologySubtask) + if remapObsClimatologySubtask is not None: + self.run_after(remapObsClimatologySubtask) # }}} def set_plot_info(self, outFileLabel, fieldNameInTitle, mpasFieldName, - obsFieldName, observationTitleLabel, unitsLabel, + refFieldName, refTitleLabel, unitsLabel, imageCaption, galleryGroup, groupSubtitle, groupLink, galleryName, diffTitleLabel='Model - Observations'): # {{{ @@ -182,11 +193,12 @@ def set_plot_info(self, outFileLabel, fieldNameInTitle, mpasFieldName, mpasFieldName : str The name of the variable in the MPAS timeSeriesStatsMonthly output - obsFieldName : str - The name of the variable to use from the observations file + refFieldName : str + The name of the variable to use from the observations or reference + file - observationTitleLabel : str - the title of the observations subplot + refTitleLabel : str + the title of the observations or reference subplot unitsLabel : str the units of the plotted field, to be displayed on color bars @@ -218,8 +230,8 @@ def set_plot_info(self, outFileLabel, fieldNameInTitle, mpasFieldName, self.outFileLabel = outFileLabel self.fieldNameInTitle = fieldNameInTitle self.mpasFieldName = mpasFieldName - self.obsFieldName = obsFieldName - self.observationTitleLabel = observationTitleLabel + self.refFieldName = refFieldName + self.refTitleLabel = refTitleLabel self.diffTitleLabel = diffTitleLabel self.unitsLabel = unitsLabel @@ -290,6 +302,7 @@ def setup_and_check(self): # {{{ def run_task(self): # {{{ """ Plots a comparison of ACME/MPAS output to SST, MLD or SSS observations + or a reference run Authors ------- @@ -304,9 +317,9 @@ def run_task(self): # {{{ season, comparisonGridName)) # first read the model climatology - remappedFileName = self.remapMpasClimatologySubtask.get_file_name( - season=season, stage='remapped', - comparisonGridName=comparisonGridName) + remappedFileName = \ + self.remapMpasClimatologySubtask.get_remapped_file_name( + season=season, comparisonGridName=comparisonGridName) remappedModelClimatology = xr.open_dataset(remappedFileName) @@ -320,31 +333,49 @@ def run_task(self): # {{{ remappedModelClimatology = remappedModelClimatology.sel( depthSlice=str(depth), drop=True) - # now the observations - remappedFileName = self.remapObsClimatologySubtask.get_file_name( - stage='remapped', season=season, - comparisonGridName=comparisonGridName) + # now the observations or reference run + if self.remapObsClimatologySubtask is not None: + remappedFileName = self.remapObsClimatologySubtask.get_file_name( + stage='remapped', season=season, + comparisonGridName=comparisonGridName) + + remappedRefClimatology = xr.open_dataset(remappedFileName) + elif self.refConfig is not None: + climatologyName = self.remapMpasClimatologySubtask.climatologyName + remappedFileName = \ + get_remapped_mpas_climatology_file_name( + self.refConfig, season=season, + componentName=self.componentName, + climatologyName=climatologyName, + comparisonGridName=comparisonGridName) + remappedRefClimatology = xr.open_dataset(remappedFileName) + refStartYear = self.refConfig.getint('climatology', 'startYear') + refEndYear = self.refConfig.getint('climatology', 'endYear') + if refStartYear != self.startYear or refEndYear != self.endYear: + self.refTitleLabel = '{}\n(years {:04d}-{:04d})'.format( + self.refTitleLabel, refStartYear, refEndYear) - remappedObsClimatology = xr.open_dataset(remappedFileName) + else: + remappedRefClimatology = None - if depth is not None: - if str(depth) not in remappedObsClimatology.depthSlice.values: + if remappedRefClimatology is not None and depth is not None: + if str(depth) not in remappedRefClimatology.depthSlice.values: raise KeyError('The climatology you are attempting to perform ' 'depth slices of was originally created\n' 'without depth {}. You will need to delete and ' 'regenerate the climatology'.format(depth)) - remappedObsClimatology = remappedObsClimatology.sel( + remappedRefClimatology = remappedRefClimatology.sel( depthSlice=str(depth), drop=True) if self.comparisonGridName == 'latlon': - self._plot_latlon(remappedModelClimatology, remappedObsClimatology) + self._plot_latlon(remappedModelClimatology, remappedRefClimatology) elif self.comparisonGridName == 'antarctic': self._plot_antarctic(remappedModelClimatology, - remappedObsClimatology) + remappedRefClimatology) # }}} - def _plot_latlon(self, remappedModelClimatology, remappedObsClimatology): + def _plot_latlon(self, remappedModelClimatology, remappedRefClimatology): # {{{ """ plotting a global lat-lon data set """ @@ -367,10 +398,14 @@ def _plot_latlon(self, remappedModelClimatology, remappedObsClimatology): lonTarg, latTarg = np.meshgrid(lon, lat) - observations = nans_to_numpy_mask( - remappedObsClimatology[self.obsFieldName].values) + if remappedRefClimatology is None: + refOutput = None + bias = None + else: + refOutput = nans_to_numpy_mask( + remappedRefClimatology[self.refFieldName].values) - bias = modelOutput - observations + bias = modelOutput - refOutput filePrefix = self.filePrefix outFileName = '{}/{}.png'.format(self.plotsDirectory, filePrefix) @@ -381,7 +416,7 @@ def _plot_latlon(self, remappedModelClimatology, remappedObsClimatology): lonTarg, latTarg, modelOutput, - observations, + refOutput, bias, colormapResult, colorbarLevelsResult, @@ -390,7 +425,7 @@ def _plot_latlon(self, remappedModelClimatology, remappedObsClimatology): fileout=outFileName, title=title, modelTitle='{}'.format(mainRunName), - obsTitle=self.observationTitleLabel, + refTitle=self.refTitleLabel, diffTitle=self.diffTitleLabel, cbarlabel=self.unitsLabel) @@ -411,7 +446,7 @@ def _plot_latlon(self, remappedModelClimatology, remappedObsClimatology): # }}} def _plot_antarctic(self, remappedModelClimatology, - remappedObsClimatology): # {{{ + remappedRefClimatology): # {{{ """ plotting an Antarctic data set """ season = self.season @@ -429,10 +464,14 @@ def _plot_antarctic(self, remappedModelClimatology, modelOutput = nans_to_numpy_mask( remappedModelClimatology[self.mpasFieldName].values) - observations = nans_to_numpy_mask( - remappedObsClimatology[self.obsFieldName].values) + if remappedRefClimatology is None: + refOutput = None + bias = None + else: + refOutput = nans_to_numpy_mask( + remappedRefClimatology[self.refFieldName].values) - bias = modelOutput - observations + bias = modelOutput - refOutput x = interp_extrap_corner(remappedModelClimatology['x'].values) y = interp_extrap_corner(remappedModelClimatology['y'].values) @@ -458,14 +497,14 @@ def _plot_antarctic(self, remappedModelClimatology, y, self.landMask, modelOutput, - observations, + refOutput, bias, fileout=outFileName, colorMapSectionName=configSectionName, colorMapType=colorMapType, title=title, modelTitle='{}'.format(mainRunName), - obsTitle=self.observationTitleLabel, + refTitle=self.refTitleLabel, diffTitle=self.diffTitleLabel, cbarlabel=self.unitsLabel) diff --git a/mpas_analysis/ocean/plot_depth_integrated_time_series_subtask.py b/mpas_analysis/ocean/plot_depth_integrated_time_series_subtask.py index 3ff1dc013..10f3b1bf2 100644 --- a/mpas_analysis/ocean/plot_depth_integrated_time_series_subtask.py +++ b/mpas_analysis/ocean/plot_depth_integrated_time_series_subtask.py @@ -72,6 +72,8 @@ class PlotDepthIntegratedTimeSeriesSubtask(AnalysisTask): galleryName : str The name of the gallery in which this plot belongs + refConfig : ``MpasAnalysisConfigParser`` + The configuration options for the reference run (if any) Authors ------- @@ -81,7 +83,8 @@ class PlotDepthIntegratedTimeSeriesSubtask(AnalysisTask): def __init__(self, parentTask, regionName, inFileName, outFileLabel, fieldNameInTitle, mpasFieldName, yAxisLabel, sectionName, thumbnailSuffix, imageCaption, galleryGroup, groupSubtitle, - groupLink, galleryName, subtaskName=None): # {{{ + groupLink, galleryName, subtaskName=None, refConfig=None): + # {{{ """ Construct the analysis task. @@ -133,9 +136,12 @@ def __init__(self, parentTask, regionName, inFileName, outFileLabel, galleryName : str the name of the gallery in which this plot belongs - subtaskName : str + subtaskName : str, optional The name of the subtask (``plotTimeSeries`` by default) + refConfig : ``MpasAnalysisConfigParser``, optional + The configuration options for the reference run (if any) + Authors ------- Xylar Asay-Davis @@ -161,6 +167,8 @@ def __init__(self, parentTask, regionName, inFileName, outFileLabel, self.yAxisLabel = yAxisLabel self.sectionName = sectionName + self.refConfig = refConfig + # xml/html related variables self.thumbnailSuffix = thumbnailSuffix self.imageCaption = imageCaption @@ -188,6 +196,17 @@ def setup_and_check(self): # {{{ config = self.config + if self.refConfig is not None: + # we need to know what file to read from the reference run so + # an absolute path won't work + assert(not os.path.isabs(self.inFileName)) + + baseDirectory = build_config_full_path( + self.refConfig, 'output', 'timeSeriesSubdirectory') + + self.refFileName = '{}/{}'.format(baseDirectory, + self.inFileName) + if not os.path.isabs(self.inFileName): baseDirectory = build_config_full_path( config, 'output', 'timeSeriesSubdirectory') @@ -234,28 +253,13 @@ def run_task(self): # {{{ self.logger.info(' Load ocean data...') ds = open_mpas_dataset(fileName=self.inFileName, calendar=calendar, - variableList=[self.mpasFieldName], + variableList=[self.mpasFieldName, 'depth'], timeVariableNames=None, startDate=startDate, endDate=endDate) ds = ds.isel(nOceanRegionsTmp=regionIndex) - # Note: restart file, not a mesh file because we need refBottomDepth, - # not in a mesh file - try: - restartFile = self.runStreams.readpath('restart')[0] - except ValueError: - raise IOError('No MPAS-O restart file found: need at least one ' - 'restart file for OHC calculation') - - # Define/read in general variables - self.logger.info(' Read in depth...') - with xr.open_dataset(restartFile) as dsRestart: - # reference depth [m] - depths = dsRestart.refBottomDepth.values - - # add depths as a coordinate to the data set - ds.coords['depth'] = (('nVertLevels',), depths) + depths = ds.depth.values divisionDepths = config.getExpression(self.sectionName, 'depths') @@ -374,6 +378,40 @@ def run_task(self): # {{{ maxPoints.append(points[rangeIndex]) legendText.append(None) + if self.refConfig is not None: + + refRunName = self.refConfig.get('runs', 'mainRunName') + + title = '{} \n {} (blue)'.format(title, refRunName) + + self.logger.info(' Load ocean data from reference run...') + refStartYear = self.refConfig.getint('timeSeries', 'startYear') + refEndYear = self.refConfig.getint('timeSeries', 'endYear') + refStartDate = '{:04d}-01-01_00:00:00'.format(refStartYear) + refEndDate = '{:04d}-12-31_23:59:59'.format(refEndYear) + dsRef = open_mpas_dataset(fileName=self.refFileName, + calendar=calendar, + variableList=[self.mpasFieldName, + 'depth'], + timeVariableNames=None, + startDate=refStartDate, + endDate=refEndDate) + dsRef = dsRef.isel(nOceanRegionsTmp=regionIndex) + + color = 'b' + + for rangeIndex in range(len(topDepths)): + top = topDepths[rangeIndex] + bottom = bottomDepths[rangeIndex] + field = dsRef[self.mpasFieldName].where(dsRef.depth > top) + field = field.where(dsRef.depth <= bottom) + timeSeries.append(field.sum('nVertLevels')) + + lineStyles.append('{}{}'.format(color, lines[rangeIndex])) + lineWidths.append(widths[rangeIndex]) + maxPoints.append(points[rangeIndex]) + legendText.append(None) + timeseries_analysis_plot(config=config, dsvalues=timeSeries, N=None, title=title, xlabel=xLabel, ylabel=yLabel, fileout=figureName, lineStyles=lineStyles, diff --git a/mpas_analysis/ocean/remap_depth_slices_subtask.py b/mpas_analysis/ocean/remap_depth_slices_subtask.py index e15fa20fa..6788fe43e 100644 --- a/mpas_analysis/ocean/remap_depth_slices_subtask.py +++ b/mpas_analysis/ocean/remap_depth_slices_subtask.py @@ -88,7 +88,6 @@ def __init__(self, mpasClimatologyTask, parentTask, climatologyName, mpasClimatologyTask, parentTask, climatologyName, variableList, seasons, comparisonGridNames, iselValues) - def run_task(self): # {{{ """ Compute climatologies of T or S from ACME/MPAS output diff --git a/mpas_analysis/ocean/time_series_antarctic_melt.py b/mpas_analysis/ocean/time_series_antarctic_melt.py index a50751eb1..8c668826b 100644 --- a/mpas_analysis/ocean/time_series_antarctic_melt.py +++ b/mpas_analysis/ocean/time_series_antarctic_melt.py @@ -11,7 +11,7 @@ from ..shared.plot.plotting import timeseries_analysis_plot -from ..shared.io import open_mpas_dataset +from ..shared.io import open_mpas_dataset, write_netcdf from ..shared.io.utility import build_config_full_path, make_directories @@ -31,23 +31,30 @@ class TimeSeriesAntarcticMelt(AnalysisTask): mpasTimeSeriesTask : ``MpasTimeSeriesTask`` The task that extracts the time series from MPAS monthly output + refConfig : ``MpasAnalysisConfigParser`` + The configuration options for the reference run (if any) + Authors ------- Xylar Asay-Davis, Stephen Price """ - def __init__(self, config, mpasTimeSeriesTask): # {{{ + def __init__(self, config, mpasTimeSeriesTask, refConfig=None): + # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasTimeSeriesTask : ``MpasTimeSeriesTask`` The task that extracts the time series from MPAS monthly output + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -62,6 +69,7 @@ def __init__(self, config, mpasTimeSeriesTask): # {{{ self.mpasTimeSeriesTask = mpasTimeSeriesTask self.run_after(mpasTimeSeriesTask) + self.refConfig = refConfig # }}} @@ -121,7 +129,7 @@ def setup_and_check(self): # {{{ self.startDate = config.get('timeSeries', 'startDate') self.endDate = config.get('timeSeries', 'endDate') - self.inputFile = self.mpasTimeSeriesTask.outputFile + self.outFileName = 'iceShelfAggregatedFluxes.nc' self.variableList = \ ['timeMonthly_avg_landIceFreshwaterFlux'] @@ -168,7 +176,7 @@ def run_task(self): # {{{ Authors ------- - Xylar Asay-Davis + Xylar Asay-Davis, Stephen Price """ self.logger.info("\nPlotting Antarctic melt rate time series...") @@ -178,12 +186,14 @@ def run_task(self): # {{{ config = self.config calendar = self.calendar - # Load data: - ds = open_mpas_dataset(fileName=self.inputFile, - calendar=calendar, - variableList=self.variableList, - startDate=self.startDate, - endDate=self.endDate) + totalMeltFlux, meltRates = self._compute_ice_shelf_fluxes() + + plotRef = self.refConfig is not None + if plotRef: + refRunName = self.refConfig.get('runs', 'mainRunName') + + refTotalMeltFlux, refMeltRates = \ + self._load_ice_shelf_fluxes(self.refConfig) # Load observations from multiple files and put in dictionary based # on shelf keyname @@ -221,36 +231,11 @@ def run_task(self): # {{{ # If areas from obs file used need to be converted from sq km to sq m - # work on data from simulations - freshwaterFlux = ds.timeMonthly_avg_landIceFreshwaterFlux - mainRunName = config.get('runs', 'mainRunName') movingAverageMonths = config.getint('timeSeriesAntarcticMelt', 'movingAverageMonths') - iceShelvesToPlot = self.iceShelvesToPlot - - dsRestart = xarray.open_dataset(self.restartFileName) - areaCell = dsRestart.landIceFraction.isel(Time=0)*dsRestart.areaCell - - dsRegionMask = xarray.open_dataset(self.regionMaskFileName) - - # select only those regions we want to plot - dsRegionMask = dsRegionMask.isel(nRegions=self.regionIndices) - cellMasks = dsRegionMask.regionCellMasks - nRegions = dsRegionMask.dims['nRegions'] - - # convert from kg/s to kg/yr - totalMeltFlux = constants.sec_per_year * \ - (cellMasks*areaCell*freshwaterFlux).sum(dim='nCells') - - totalArea = (cellMasks*areaCell).sum(dim='nCells') - - # from kg/m^2/yr to m/yr - meltRates = (1./constants.rho_fw) * (totalMeltFlux/totalArea) - - # convert from kg/yr to GT/yr - totalMeltFlux /= constants.kg_per_GT + nRegions = totalMeltFlux.sizes['nRegions'] outputDirectory = build_config_full_path(config, 'output', 'timeseriesSubdirectory') @@ -260,7 +245,7 @@ def run_task(self): # {{{ self.logger.info(' Make plots...') for iRegion in range(nRegions): - regionName = iceShelvesToPlot[iRegion] + regionName = self.iceShelvesToPlot[iRegion] # get obs melt flux and unc. for shelf (similar for rates) obsMeltFlux = [] @@ -298,10 +283,21 @@ def run_task(self): # {{{ filePrefix = 'melt_flux_{}'.format(regionName) figureName = '{}/{}.png'.format(self.plotsDirectory, filePrefix) - timeseries_analysis_plot(config, [timeSeries], movingAverageMonths, + fields = [timeSeries] + lineStyles = ['k-'] + lineWidths = [2.5] + legendText = [mainRunName] + if plotRef: + fields.append(refTotalMeltFlux.isel(nRegions=iRegion)) + lineStyles.append('r-') + lineWidths.append(1.2) + legendText.append(refRunName) + + timeseries_analysis_plot(config, fields, movingAverageMonths, title, xLabel, yLabel, figureName, - lineStyles=['b-'], lineWidths=[1.2], - legendText=[mainRunName], + lineStyles=lineStyles, + lineWidths=lineWidths, + legendText=legendText, calendar=calendar, obsMean=obsMeltFlux, obsUncertainty=obsMeltFluxUnc, obsLegend=list(obsDict.keys())) @@ -328,10 +324,21 @@ def run_task(self): # {{{ filePrefix = 'melt_rate_{}'.format(regionName) figureName = '{}/{}.png'.format(self.plotsDirectory, filePrefix) - timeseries_analysis_plot(config, [timeSeries], movingAverageMonths, + fields = [timeSeries] + lineStyles = ['k-'] + lineWidths = [2.5] + legendText = [mainRunName] + if plotRef: + fields.append(refMeltRates.isel(nRegions=iRegion)) + lineStyles.append('r-') + lineWidths.append(1.2) + legendText.append(refRunName) + + timeseries_analysis_plot(config, fields, movingAverageMonths, title, xLabel, yLabel, figureName, - lineStyles=['b-'], lineWidths=[1.2], - legendText=[mainRunName], + lineStyles=lineStyles, + lineWidths=lineWidths, + legendText=legendText, calendar=calendar, obsMean=obsMeltRate, obsUncertainty=obsMeltRateUnc, obsLegend=list(obsDict.keys())) @@ -351,6 +358,118 @@ def run_task(self): # {{{ imageCaption=caption) # }}} + def _compute_ice_shelf_fluxes(self): # {{{ + """ + Reads melt flux time series and computes regional total melt flux and + mean melt rate. + + Authors + ------- + Xylar Asay-Davis, Stephen Price + """ + + mpasTimeSeriesTask = self.mpasTimeSeriesTask + config = self.config + + baseDirectory = build_config_full_path( + config, 'output', 'timeSeriesSubdirectory') + + outFileName = '{}/{}'.format(baseDirectory, self.outFileName) + + # Load data: + inputFile = mpasTimeSeriesTask.outputFile + dsIn = open_mpas_dataset(fileName=inputFile, + calendar=self.calendar, + variableList=self.variableList, + startDate=self.startDate, + endDate=self.endDate) + try: + if os.path.exists(outFileName): + # The file already exists so load it + dsOut = xarray.open_dataset(outFileName) + if numpy.all(dsOut.Time.values == dsIn.Time.values): + return dsOut.totalMeltFlux, dsOut.meltRates + else: + self.logger.warning('File {} is incomplete. Deleting ' + 'it.'.format(outFileName)) + os.remove(outFileName) + except OSError: + # something is potentailly wrong with the file, so let's delete + # it and try again + self.logger.warning('Problems reading file {}. Deleting ' + 'it.'.format(outFileName)) + os.remove(outFileName) + + # work on data from simulations + freshwaterFlux = dsIn.timeMonthly_avg_landIceFreshwaterFlux + + restartFileName = \ + mpasTimeSeriesTask.runStreams.readpath('restart')[0] + + dsRestart = xarray.open_dataset(restartFileName) + areaCell = dsRestart.landIceFraction.isel(Time=0)*dsRestart.areaCell + + mpasMeshName = config.get('input', 'mpasMeshName') + regionMaskDirectory = config.get('regions', 'regionMaskDirectory') + + regionMaskFileName = '{}/{}_iceShelfMasks.nc'.format( + regionMaskDirectory, mpasMeshName) + + dsRegionMask = xarray.open_dataset(regionMaskFileName) + + # select only those regions we want to plot + dsRegionMask = dsRegionMask.isel(nRegions=self.regionIndices) + cellMasks = dsRegionMask.regionCellMasks + + # convert from kg/s to kg/yr + totalMeltFlux = constants.sec_per_year * \ + (cellMasks*areaCell*freshwaterFlux).sum(dim='nCells') + + totalArea = (cellMasks*areaCell).sum(dim='nCells') + + # from kg/m^2/yr to m/yr + meltRates = (1./constants.rho_fw) * (totalMeltFlux/totalArea) + + # convert from kg/yr to GT/yr + totalMeltFlux /= constants.kg_per_GT + + baseDirectory = build_config_full_path( + config, 'output', 'timeSeriesSubdirectory') + + outFileName = '{}/iceShelfAggregatedFluxes.nc'.format(baseDirectory) + + dsOut = xarray.Dataset() + dsOut['totalMeltFlux'] = totalMeltFlux + dsOut.totalMeltFlux.attrs['units'] = 'GT a$^{-1}$' + dsOut.totalMeltFlux.attrs['description'] = \ + 'Total melt flux summed over each ice shelf or region' + dsOut['meltRates'] = meltRates + dsOut.meltRates.attrs['units'] = 'm a$^{-1}$' + dsOut.meltRates.attrs['description'] = \ + 'Melt rate averaged over each ice shelf or region' + + write_netcdf(dsOut, outFileName) + + return totalMeltFlux, meltRates # }}} + + def _load_ice_shelf_fluxes(self, config): # {{{ + """ + Reads melt flux time series and computes regional total melt flux and + mean melt rate. + + Authors + ------- + Xylar Asay-Davis + """ + baseDirectory = build_config_full_path( + config, 'output', 'timeSeriesSubdirectory') + + outFileName = '{}/{}'.format(baseDirectory, self.outFileName) + + dsOut = xarray.open_dataset(outFileName) + return dsOut.totalMeltFlux, dsOut.meltRates + + # }}} # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/time_series_ohc_anomaly.py b/mpas_analysis/ocean/time_series_ohc_anomaly.py index cab550a5e..0592abaf5 100644 --- a/mpas_analysis/ocean/time_series_ohc_anomaly.py +++ b/mpas_analysis/ocean/time_series_ohc_anomaly.py @@ -3,6 +3,8 @@ from __future__ import absolute_import, division, print_function, \ unicode_literals +import xarray as xr + from ..shared import AnalysisTask from .compute_anomaly_subtask import ComputeAnomalySubtask @@ -20,18 +22,22 @@ class TimeSeriesOHCAnomaly(AnalysisTask): Xylar Asay-Davis, Milena Veneziani, Greg Streletz """ - def __init__(self, config, mpasTimeSeriesTask): # {{{ + def __init__(self, config, mpasTimeSeriesTask, refConfig=None): + # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasTimeSeriesTask : ``MpasTimeSeriesTask`` The task that extracts the time series from MPAS monthly output + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -104,7 +110,8 @@ def __init__(self, config, mpasTimeSeriesTask): # {{{ galleryGroup='Time Series', groupSubtitle=None, groupLink='timeseries', - galleryName=None) + galleryName=None, + refConfig=refConfig) plotTask.run_after(anomalyTask) self.add_subtask(plotTask) @@ -135,6 +142,21 @@ def _compute_ohc(self, ds): # {{{ ds.ohc.attrs['units'] = '$10^{22}$ J' ds.ohc.attrs['description'] = 'Ocean heat content in each region' + # Note: restart file, not a mesh file because we need refBottomDepth, + # not in a mesh file + try: + restartFile = self.runStreams.readpath('restart')[0] + except ValueError: + raise IOError('No MPAS-O restart file found: need at least one ' + 'restart file for OHC calculation') + + # Define/read in general variables + with xr.open_dataset(restartFile) as dsRestart: + # reference depth [m] + # add depths as a coordinate to the data set + ds.coords['depth'] = (('nVertLevels',), + dsRestart.refBottomDepth.values) + return ds # }}} # }}} diff --git a/mpas_analysis/ocean/time_series_sst.py b/mpas_analysis/ocean/time_series_sst.py index bcdeda694..947d2b881 100644 --- a/mpas_analysis/ocean/time_series_sst.py +++ b/mpas_analysis/ocean/time_series_sst.py @@ -27,23 +27,30 @@ class TimeSeriesSST(AnalysisTask): mpasTimeSeriesTask : ``MpasTimeSeriesTask`` The task that extracts the time series from MPAS monthly output + refConfig : ``MpasAnalysisConfigParser`` + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis, Milena Veneziani """ - def __init__(self, config, mpasTimeSeriesTask): # {{{ + def __init__(self, config, mpasTimeSeriesTask, refConfig=None): + # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasTimeSeriesTask : ``MpasTimeSeriesTask`` The task that extracts the time series from MPAS monthly output + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -56,6 +63,7 @@ def __init__(self, config, mpasTimeSeriesTask): # {{{ tags=['timeSeries', 'sst']) self.mpasTimeSeriesTask = mpasTimeSeriesTask + self.refConfig = refConfig self.run_after(mpasTimeSeriesTask) @@ -160,6 +168,27 @@ def run_task(self): # {{{ timeEnd = date_to_days(year=yearEnd, month=12, day=31, calendar=calendar) + if self.refConfig is not None: + baseDirectory = build_config_full_path( + self.refConfig, 'output', 'timeSeriesSubdirectory') + + refFileName = '{}/{}.nc'.format( + baseDirectory, self.mpasTimeSeriesTask.fullTaskName) + + refStartYear = self.refConfig.getint('timeSeries', 'startYear') + refEndYear = self.refConfig.getint('timeSeries', 'endYear') + refStartDate = '{:04d}-01-01_00:00:00'.format(refStartYear) + refEndDate = '{:04d}-12-31_23:59:59'.format(refEndYear) + + dsRefSST = open_mpas_dataset( + fileName=refFileName, + calendar=calendar, + variableList=self.variableList, + startDate=refStartDate, + endDate=refEndDate) + else: + dsRefSST = None + if preprocessedReferenceRunName != 'None': self.logger.info(' Load in SST for a preprocesses reference ' 'run...') @@ -185,8 +214,7 @@ def run_task(self): # {{{ for regionIndex in regionIndicesToPlot: region = regions[regionIndex] - title = plotTitles[regionIndex] - title = 'SST, %s \n %s (black)' % (title, mainRunName) + title = '{} SST'.format(plotTitles[regionIndex]) xLabel = 'Time [years]' yLabel = '[$^\circ$C]' @@ -197,23 +225,32 @@ def run_task(self): # {{{ figureName = '{}/{}.png'.format(self.plotsDirectory, filePrefix) + lineStyles = ['k-'] + lineWidths = [3] + + fields = [SST] + legendText = [mainRunName] + + if dsRefSST is not None: + refSST = dsRefSST[varName].isel(nOceanRegions=regionIndex) + fields.append(refSST) + lineStyles.append('b-') + lineWidths.append(1.5) + refRunName = self.refConfig.get('runs', 'mainRunName') + legendText.append(refRunName) + if preprocessedReferenceRunName != 'None': SST_v0 = dsPreprocessedTimeSlice.SST - - title = '{}\n {} (red)'.format(title, - preprocessedReferenceRunName) - timeseries_analysis_plot(config, [SST, SST_v0], - movingAveragePoints, - title, xLabel, yLabel, figureName, - lineStyles=['k-', 'r-'], - lineWidths=[3, 1.5], - legendText=[None, None], - calendar=calendar) - else: - timeseries_analysis_plot(config, [SST], movingAveragePoints, - title, xLabel, yLabel, figureName, - lineStyles=['k-'], lineWidths=[3], - legendText=[None], calendar=calendar) + fields.append(SST_v0) + lineStyles.append('r-') + lineWidths.append(1.5) + legendText.append(preprocessedReferenceRunName) + + timeseries_analysis_plot(config, fields, movingAveragePoints, + title, xLabel, yLabel, figureName, + lineStyles=lineStyles, + lineWidths=lineWidths, + legendText=legendText, calendar=calendar) caption = 'Running Mean of {} Sea Surface Temperature'.format( region) diff --git a/mpas_analysis/sea_ice/climatology_map_sea_ice_conc.py b/mpas_analysis/sea_ice/climatology_map_sea_ice_conc.py index 9f26dfc70..c749f5c9c 100644 --- a/mpas_analysis/sea_ice/climatology_map_sea_ice_conc.py +++ b/mpas_analysis/sea_ice/climatology_map_sea_ice_conc.py @@ -25,15 +25,15 @@ class ClimatologyMapSeaIceConc(AnalysisTask): # {{{ ------- Luke Van Roekel, Xylar Asay-Davis, Milena Veneziani """ - def __init__(self, config, mpasClimatologyTask, hemisphere): - # {{{ + def __init__(self, config, mpasClimatologyTask, hemisphere, + refConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted @@ -41,18 +41,21 @@ def __init__(self, config, mpasClimatologyTask, hemisphere): hemisphere : {'NH', 'SH'} The hemisphere to plot + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis """ taskName = 'climatologyMapSeaIceConc{}'.format(hemisphere) - self.fieldName = 'seaIceConc' + fieldName = 'seaIceConc' # call the constructor from the base class (AnalysisTask) super(ClimatologyMapSeaIceConc, self).__init__( config=config, taskName=taskName, componentName='seaIce', - tags=['climatology', 'horizontalMap', self.fieldName]) + tags=['climatology', 'horizontalMap', fieldName]) mpasFieldName = 'timeMonthly_avg_iceAreaCell' iselValues = None @@ -64,8 +67,6 @@ def __init__(self, config, mpasClimatologyTask, hemisphere): else: hemisphereLong = 'Southern' - obsFieldName = 'seaIceConc' - # read in what seasons we want to plot seasons = config.getExpression(sectionName, 'seasons') @@ -85,13 +86,30 @@ def __init__(self, config, mpasClimatologyTask, hemisphere): remapClimatologySubtask = RemapMpasClimatologySubtask( mpasClimatologyTask=mpasClimatologyTask, parentTask=self, - climatologyName='{}{}'.format(self.fieldName, hemisphere), + climatologyName='{}{}'.format(fieldName, hemisphere), variableList=[mpasFieldName], comparisonGridNames=comparisonGridNames, seasons=seasons, iselValues=iselValues) - observationPrefixes = config.getExpression(sectionName, + if refConfig is None: + self._add_obs_tasks(seasons, comparisonGridNames, hemisphere, + hemisphereLong, remapClimatologySubtask, + mpasFieldName) + else: + self._add_ref_tasks(seasons, comparisonGridNames, hemisphere, + hemisphereLong, remapClimatologySubtask, + refConfig, mpasFieldName, + fieldName, iselValues) + # }}} + + def _add_obs_tasks(self, seasons, comparisonGridNames, hemisphere, + hemisphereLong, remapClimatologySubtask, + mpasFieldName): # {{{ + config = self.config + obsFieldName = 'seaIceConc' + + observationPrefixes = config.getExpression(self.taskName, 'observationPrefixes') for prefix in observationPrefixes: for season in seasons: @@ -107,46 +125,93 @@ def __init__(self, config, mpasClimatologyTask, hemisphere): remapObservationsSubtask = RemapObservedConcClimatology( parentTask=self, seasons=[season], fileName=obsFileName, - outFilePrefix='{}{}{}_{}'.format(obsFieldName, prefix, - hemisphere, season), + outFilePrefix='{}{}{}_{}'.format( + obsFieldName, prefix, hemisphere, season), comparisonGridNames=comparisonGridNames, - subtaskName='remapObservations_{}{}'.format(prefix, - season)) + subtaskName='remapObservations_{}{}'.format( + prefix, season)) self.add_subtask(remapObservationsSubtask) for comparisonGridName in comparisonGridNames: imageDescription = \ '{} Climatology Map of {}-Hemisphere Sea-Ice ' \ 'Concentration'.format(season, hemisphereLong) - imageCaption = '{}.
Observations: SSM/I {}'.format( - imageDescription, prefix) + imageCaption = \ + '{}.
Observations: SSM/I {}'.format( + imageDescription, prefix) galleryGroup = \ '{}-Hemisphere Sea-Ice Concentration'.format( hemisphereLong) - # make a new subtask for this season and comparison grid + # make a new subtask for this season and comparison + # grid subtask = PlotClimatologyMapSubtask( - self, hemisphere, season, comparisonGridName, - remapClimatologySubtask, remapObservationsSubtask, - subtaskSuffix=prefix) + self, hemisphere, season, comparisonGridName, + remapClimatologySubtask, + remapObsClimatologySubtask=remapObservationsSubtask, + subtaskSuffix=prefix) subtask.set_plot_info( - outFileLabel='iceconc{}{}'.format(prefix, - hemisphere), - fieldNameInTitle='Sea ice concentration', - mpasFieldName=mpasFieldName, - obsFieldName=obsFieldName, - observationTitleLabel=observationTitleLabel, - unitsLabel=r'fraction', - imageDescription=imageDescription, - imageCaption=imageCaption, - galleryGroup=galleryGroup, - groupSubtitle=None, - groupLink='{}_conc'.format(hemisphere.lower()), - galleryName='Observations: SSM/I {}'.format( - prefix)) + outFileLabel='iceconc{}{}'.format(prefix, + hemisphere), + fieldNameInTitle='Sea ice concentration', + mpasFieldName=mpasFieldName, + refFieldName=obsFieldName, + refTitleLabel=observationTitleLabel, + diffTitleLabel='Model - Observations', + unitsLabel=r'fraction', + imageDescription=imageDescription, + imageCaption=imageCaption, + galleryGroup=galleryGroup, + groupSubtitle=None, + groupLink='{}_conc'.format(hemisphere.lower()), + galleryName='Observations: SSM/I {}'.format( + prefix)) self.add_subtask(subtask) # }}} + + def _add_ref_tasks(self, seasons, comparisonGridNames, hemisphere, + hemisphereLong, remapClimatologySubtask, + refConfig, mpasFieldName, fieldName, + iselValues): # {{{ + + refRunName = refConfig.get('runs', 'mainRunName') + galleryName = None + refTitleLabel = 'Ref: {}'.format(refRunName) + + for season in seasons: + for comparisonGridName in comparisonGridNames: + + imageDescription = \ + '{} Climatology Map of {}-Hemisphere Sea-Ice ' \ + 'Concentration'.format(season, hemisphereLong) + imageCaption = imageDescription + galleryGroup = \ + '{}-Hemisphere Sea-Ice Concentration'.format( + hemisphereLong) + # make a new subtask for this season and comparison + # grid + subtask = PlotClimatologyMapSubtask( + self, hemisphere, season, comparisonGridName, + remapClimatologySubtask, refConfig=refConfig) + + subtask.set_plot_info( + outFileLabel='iceconc{}'.format(hemisphere), + fieldNameInTitle='Sea ice concentration', + mpasFieldName=mpasFieldName, + refFieldName=mpasFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel='Main - Reference', + unitsLabel=r'fraction', + imageDescription=imageDescription, + imageCaption=imageCaption, + galleryGroup=galleryGroup, + groupSubtitle=None, + groupLink='{}_conc'.format(hemisphere.lower()), + galleryName=galleryName) + + self.add_subtask(subtask) + # }}} # }}} diff --git a/mpas_analysis/sea_ice/climatology_map_sea_ice_thick.py b/mpas_analysis/sea_ice/climatology_map_sea_ice_thick.py index 503d9b24e..96c3f4f04 100644 --- a/mpas_analysis/sea_ice/climatology_map_sea_ice_thick.py +++ b/mpas_analysis/sea_ice/climatology_map_sea_ice_thick.py @@ -25,15 +25,15 @@ class ClimatologyMapSeaIceThick(AnalysisTask): # {{{ ------- Luke Van Roekel, Xylar Asay-Davis, Milena Veneziani """ - def __init__(self, config, mpasClimatologyTask, hemisphere): - # {{{ + def __init__(self, config, mpasClimatologyTask, hemisphere, + refConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasClimatologyTask : ``MpasClimatologyTask`` The task that produced the climatology to be remapped and plotted @@ -41,18 +41,21 @@ def __init__(self, config, mpasClimatologyTask, hemisphere): hemisphere : {'NH', 'SH'} The hemisphere to plot + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis """ taskName = 'climatologyMapSeaIceThick{}'.format(hemisphere) - self.fieldName = 'seaIceThick' + fieldName = 'seaIceThick' # call the constructor from the base class (AnalysisTask) super(ClimatologyMapSeaIceThick, self).__init__( config=config, taskName=taskName, componentName='seaIce', - tags=['climatology', 'horizontalMap', self.fieldName]) + tags=['climatology', 'horizontalMap', fieldName]) mpasFieldName = 'timeMonthly_avg_iceVolumeCell' iselValues = None @@ -64,8 +67,6 @@ def __init__(self, config, mpasClimatologyTask, hemisphere): else: hemisphereLong = 'Southern' - obsFieldName = 'seaIceThick' - # read in what seasons we want to plot seasons = config.getExpression(sectionName, 'seasons') @@ -85,26 +86,42 @@ def __init__(self, config, mpasClimatologyTask, hemisphere): remapClimatologySubtask = RemapMpasClimatologySubtask( mpasClimatologyTask=mpasClimatologyTask, parentTask=self, - climatologyName='{}{}'.format(self.fieldName, hemisphere), + climatologyName='{}{}'.format(fieldName, hemisphere), variableList=[mpasFieldName], comparisonGridNames=comparisonGridNames, seasons=seasons, iselValues=iselValues) + if refConfig is None: + refTitleLabel = 'Observations (ICESat)' + galleryName = 'Observations: ICESat' + diffTitleLabel = 'Model - Observations' + refFieldName = 'seaIceThick' + else: + refRunName = refConfig.get('runs', 'mainRunName') + galleryName = None + refTitleLabel = 'Ref: {}'.format(refRunName) + refFieldName = mpasFieldName + diffTitleLabel = 'Main - Reference' + + remapObservationsSubtask = None + for season in seasons: + if refConfig is None: + obsFileName = build_config_full_path( + config, 'seaIceObservations', + 'thickness{}_{}'.format(hemisphere, season)) + + remapObservationsSubtask = RemapObservedThickClimatology( + parentTask=self, seasons=[season], + fileName=obsFileName, + outFilePrefix='{}{}_{}'.format(refFieldName, + hemisphere, + season), + comparisonGridNames=comparisonGridNames, + subtaskName='remapObservations{}'.format(season)) + self.add_subtask(remapObservationsSubtask) - obsFileName = build_config_full_path( - config, 'seaIceObservations', - 'thickness{}_{}'.format(hemisphere, season)) - - remapObservationsSubtask = RemapObservedThickClimatology( - parentTask=self, seasons=[season], - fileName=obsFileName, - outFilePrefix='{}{}_{}'.format(obsFieldName, hemisphere, - season), - comparisonGridNames=comparisonGridNames, - subtaskName='remapObservations{}'.format(season)) - self.add_subtask(remapObservationsSubtask) for comparisonGridName in comparisonGridNames: imageDescription = \ @@ -117,25 +134,29 @@ def __init__(self, config, mpasClimatologyTask, hemisphere): # make a new subtask for this season and comparison grid subtask = PlotClimatologyMapSubtask( self, hemisphere, season, comparisonGridName, - remapClimatologySubtask, remapObservationsSubtask) + remapClimatologySubtask, remapObservationsSubtask, + refConfig) subtask.set_plot_info( outFileLabel='icethick{}'.format(hemisphere), fieldNameInTitle='Sea ice thickness', mpasFieldName=mpasFieldName, - obsFieldName=obsFieldName, - observationTitleLabel='Observations (ICESat)', + refFieldName=refFieldName, + refTitleLabel=refTitleLabel, + diffTitleLabel=diffTitleLabel, unitsLabel=r'm', imageDescription=imageDescription, imageCaption=imageCaption, galleryGroup=galleryGroup, groupSubtitle=None, groupLink='{}_thick'.format(hemisphere.lower()), - galleryName='Observations: ICESat', + galleryName=galleryName, maskValue=0) self.add_subtask(subtask) + # }}} + # }}} diff --git a/mpas_analysis/sea_ice/plot_climatology_map_subtask.py b/mpas_analysis/sea_ice/plot_climatology_map_subtask.py index acbf7fdbb..9d5bf9ad2 100644 --- a/mpas_analysis/sea_ice/plot_climatology_map_subtask.py +++ b/mpas_analysis/sea_ice/plot_climatology_map_subtask.py @@ -14,10 +14,13 @@ from ..shared.html import write_image_xml +from ..shared.climatology import get_remapped_mpas_climatology_file_name + class PlotClimatologyMapSubtask(AnalysisTask): # {{{ """ - An analysis task for plotting 2D model fields against observations. + An analysis task for plotting 2D model fields against observations or a + reference run. Attributes ---------- @@ -39,6 +42,9 @@ class PlotClimatologyMapSubtask(AnalysisTask): # {{{ The subtask for remapping the observational climatology that this subtask will plot + refConfig : ``MpasAnalysisConfigParser`` + Configuration options for a reference run (if any) + outFileLabel : str The prefix on each plot and associated XML file @@ -48,11 +54,14 @@ class PlotClimatologyMapSubtask(AnalysisTask): # {{{ mpasFieldName : str The name of the variable in the MPAS timeSeriesStatsMonthly output - obsFieldName : str - The name of the variable to use from the observations file + refFieldName : str + The name of the variable to use from the observations or reference file + + refTitleLabel : str + the title of the observations or reference run subplot - observationTitleLabel : str - the title of the observations subplot + diffTitleLabel : str + the title of the difference (e.g. bias) subplot unitsLabel : str the units of the plotted field, to be displayed on color bars @@ -80,8 +89,9 @@ class PlotClimatologyMapSubtask(AnalysisTask): # {{{ """ def __init__(self, parentTask, hemisphere, season, comparisonGridName, - remapMpasClimatologySubtask, remapObsClimatologySubtask, - subtaskSuffix=None): # {{{ + remapMpasClimatologySubtask, remapObsClimatologySubtask=None, + refConfig=None, subtaskSuffix=None): + # {{{ ''' Construct one analysis subtask for each plot (i.e. each season and comparison grid) and a subtask for computing climatologies. @@ -105,10 +115,13 @@ def __init__(self, parentTask, hemisphere, season, comparisonGridName, The subtask for remapping the MPAS climatology that this subtask will plot - remapObsClimatologySubtask : ``RemapObservedClimatologySubtask`` + remapObsClimatologySubtask : ``RemapObservedClimatologySubtask``, optional The subtask for remapping the observational climatology that this subtask will plot + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + subtaskSuffix : str, optional A suffix on the subtask to ensure that it is unique (e.g. the observations being plotted) @@ -124,6 +137,7 @@ def __init__(self, parentTask, hemisphere, season, comparisonGridName, self.comparisonGridName = comparisonGridName self.remapMpasClimatologySubtask = remapMpasClimatologySubtask self.remapObsClimatologySubtask = remapObsClimatologySubtask + self.refConfig = refConfig subtaskName = 'plot{}'.format(season) if subtaskSuffix is not None: @@ -136,16 +150,17 @@ def __init__(self, parentTask, hemisphere, season, comparisonGridName, # call the constructor from the base class (AnalysisTask) super(PlotClimatologyMapSubtask, self).__init__( config=config, taskName=taskName, subtaskName=subtaskName, - componentName='ocean', tags=tags) + componentName='seaIce', tags=tags) # this task should not run until the remapping subtasks are done, since # it relies on data from those subtasks self.run_after(remapMpasClimatologySubtask) - self.run_after(remapObsClimatologySubtask) + if remapObsClimatologySubtask is not None: + self.run_after(remapObsClimatologySubtask) # }}} def set_plot_info(self, outFileLabel, fieldNameInTitle, mpasFieldName, - obsFieldName, observationTitleLabel, unitsLabel, + refFieldName, refTitleLabel, diffTitleLabel, unitsLabel, imageDescription, imageCaption, galleryGroup, groupSubtitle, groupLink, galleryName, maskValue=None): # {{{ @@ -163,11 +178,15 @@ def set_plot_info(self, outFileLabel, fieldNameInTitle, mpasFieldName, mpasFieldName : str The name of the variable in the MPAS timeSeriesStatsMonthly output - obsFieldName : str - The name of the variable to use from the observations file + refFieldName : str + The name of the variable to use from the observations or reference + run file - observationTitleLabel : str - the title of the observations subplot + refTitleLabel : str + the title of the observations or reference run subplot + + diffTitleLabel : str + the title of the difference (e.g. bias) subplot unitsLabel : str the units of the plotted field, to be displayed on color bars @@ -201,8 +220,9 @@ def set_plot_info(self, outFileLabel, fieldNameInTitle, mpasFieldName, self.outFileLabel = outFileLabel self.fieldNameInTitle = fieldNameInTitle self.mpasFieldName = mpasFieldName - self.obsFieldName = obsFieldName - self.observationTitleLabel = observationTitleLabel + self.refFieldName = refFieldName + self.refTitleLabel = refTitleLabel + self.diffTitleLabel = diffTitleLabel self.unitsLabel = unitsLabel # xml/html related variables @@ -261,9 +281,9 @@ def run_task(self): # {{{ comparisonGridName = self.comparisonGridName self.logger.info("\nPlotting 2-d maps of {} climatologies for " - "{} from {} observations...".format( + "{} against {}...".format( self.fieldNameInTitle, - season, self.observationTitleLabel)) + season, self.refTitleLabel)) mainRunName = config.get('runs', 'mainRunName') startYear = self.startYear @@ -287,9 +307,9 @@ def run_task(self): # {{{ minimumLatitude = config.getfloat(sectionName, 'minimumLatitude') - remappedFileName = self.remapMpasClimatologySubtask.get_file_name( - season=season, stage='remapped', - comparisonGridName=comparisonGridName) + remappedFileName = \ + self.remapMpasClimatologySubtask.get_remapped_file_name( + season=season, comparisonGridName=comparisonGridName) remappedClimatology = xr.open_dataset(remappedFileName) modelOutput = remappedClimatology[self.mpasFieldName].values @@ -300,21 +320,47 @@ def run_task(self): # {{{ lonTarg, latTarg = np.meshgrid(lon, lat) - remappedFileName = self.remapObsClimatologySubtask.get_file_name( - stage='remapped', season=season, - comparisonGridName=comparisonGridName) - - remappedObsClimatology = xr.open_dataset(remappedFileName) + if self.remapObsClimatologySubtask is not None: + remappedFileName = self.remapObsClimatologySubtask.get_file_name( + stage='remapped', season=season, + comparisonGridName=comparisonGridName) + + remappedRefClimatology = xr.open_dataset(remappedFileName) + + refOutput = remappedRefClimatology[self.refFieldName].values + if self.maskValue is not None: + refOutput = ma.masked_values(refOutput, self.maskValue) + + difference = modelOutput - refOutput + elif self.refConfig is not None: + climatologyName = self.remapMpasClimatologySubtask.climatologyName + remappedFileName = \ + get_remapped_mpas_climatology_file_name( + self.refConfig, season=season, + componentName=self.componentName, + climatologyName=climatologyName, + comparisonGridName=comparisonGridName) + remappedRefClimatology = xr.open_dataset(remappedFileName) + refStartYear = self.refConfig.getint('climatology', 'startYear') + refEndYear = self.refConfig.getint('climatology', 'endYear') + if refStartYear != self.startYear or refEndYear != self.endYear: + self.refTitleLabel = '{}\n(years {:04d}-{:04d})'.format( + self.refTitleLabel, refStartYear, refEndYear) + else: + remappedRefClimatology = None - observations = remappedObsClimatology[self.obsFieldName].values - if self.maskValue is not None: - observations = ma.masked_values(observations, self.maskValue) + if remappedRefClimatology is None: + refOutput = None + difference = None + else: + refOutput = remappedRefClimatology[self.refFieldName].values + if self.maskValue is not None: + refOutput = ma.masked_values(refOutput, self.maskValue) - difference = modelOutput - observations + difference = modelOutput - refOutput startYear = self.startYear endYear = self.endYear - observationTitleLabel = self.observationTitleLabel filePrefix = self.filePrefix title = '{} ({}, years {:04d}-{:04d})'.format( self.fieldNameInTitle, season, startYear, endYear) @@ -324,7 +370,7 @@ def run_task(self): # {{{ lonTarg, latTarg, modelOutput, - observations, + refOutput, difference, colormapResult, colorbarLevelsResult, @@ -336,8 +382,8 @@ def run_task(self): # {{{ latmin=minimumLatitude, lon0=referenceLongitude, modelTitle=mainRunName, - obsTitle=observationTitleLabel, - diffTitle='Model-Observations', + refTitle=self.refTitleLabel, + diffTitle=self.diffTitleLabel, cbarlabel=self.unitsLabel, vertical=vertical) diff --git a/mpas_analysis/sea_ice/time_series.py b/mpas_analysis/sea_ice/time_series.py index d045daed1..c890964f3 100644 --- a/mpas_analysis/sea_ice/time_series.py +++ b/mpas_analysis/sea_ice/time_series.py @@ -3,6 +3,7 @@ unicode_literals import xarray as xr +import os from ..shared import AnalysisTask @@ -17,7 +18,7 @@ from ..shared.timekeeping.MpasRelativeDelta import MpasRelativeDelta from ..shared.generalized_reader import open_multifile_dataset -from ..shared.io import open_mpas_dataset +from ..shared.io import open_mpas_dataset, write_netcdf from ..shared.mpas_xarray.mpas_xarray import subset_variables from ..shared.html import write_image_xml @@ -33,23 +34,31 @@ class TimeSeriesSeaIce(AnalysisTask): mpasTimeSeriesTask : ``MpasTimeSeriesTask`` The task that extracts the time series from MPAS monthly output + refConfig : ``MpasAnalysisConfigParser`` + Configuration options for a reference run (if any) + + Authors ------- Xylar Asay-Davis, Milena Veneziani """ - def __init__(self, config, mpasTimeSeriesTask): # {{{ + def __init__(self, config, mpasTimeSeriesTask, + refConfig=None): # {{{ """ Construct the analysis task. Parameters ---------- - config : instance of MpasAnalysisConfigParser - Contains configuration options + config : ``MpasAnalysisConfigParser`` + Configuration options mpasTimeSeriesTask : ``MpasTimeSeriesTask`` The task that extracts the time series from MPAS monthly output + refConfig : ``MpasAnalysisConfigParser``, optional + Configuration options for a reference run (if any) + Authors ------- Xylar Asay-Davis @@ -62,6 +71,7 @@ def __init__(self, config, mpasTimeSeriesTask): # {{{ tags=['timeSeries']) self.mpasTimeSeriesTask = mpasTimeSeriesTask + self.refConfig = refConfig self.run_after(mpasTimeSeriesTask) @@ -227,21 +237,13 @@ def run_task(self): # {{{ self.logger.info(' Load sea-ice data...') # Load mesh - self.dsMesh = xr.open_dataset(self.restartFileName) - self.dsMesh = subset_variables(self.dsMesh, - variableList=['lonCell', 'latCell', - 'areaCell']) - - # Load data - ds = open_mpas_dataset( - fileName=self.inputFile, - calendar=calendar, - variableList=self.variableList, - startDate=self.startDate, - endDate=self.endDate) - - yearStart = days_to_datetime(ds.Time.min(), calendar=calendar).year - yearEnd = days_to_datetime(ds.Time.max(), calendar=calendar).year + + dsTimeSeries = self._compute_area_vol() + + yearStart = days_to_datetime(dsTimeSeries['NH'].Time.min(), + calendar=calendar).year + yearEnd = days_to_datetime(dsTimeSeries['NH'].Time.max(), + calendar=calendar).year timeStart = date_to_days(year=yearStart, month=1, day=1, calendar=calendar) timeEnd = date_to_days(year=yearEnd, month=12, day=31, @@ -266,6 +268,20 @@ def run_task(self): # {{{ 'plotted.') preprocessedReferenceRunName = 'None' + if self.refConfig is not None: + + dsTimeSeriesRef = {} + baseDirectory = build_config_full_path( + self.refConfig, 'output', 'timeSeriesSubdirectory') + + refRunName = self.refConfig.get('runs', 'mainRunName') + + for hemisphere in ['NH', 'SH']: + inFileName = '{}/seaIceAreaVol{}.nc'.format(baseDirectory, + hemisphere) + + dsTimeSeriesRef[hemisphere] = xr.open_dataset(inFileName) + norm = {'iceArea': 1e-6, # m^2 to km^2 'iceVolume': 1e-12, # m^3 to 10^3 km^3 'iceThickness': 1.} @@ -275,18 +291,17 @@ def run_task(self): # {{{ galleryGroup = 'Time Series' groupLink = 'timeseries' - dsTimeSeries = {} obs = {} preprocessed = {} figureNameStd = {} figureNamePolar = {} title = {} plotVars = {} + obsLegend = {} + plotVarsRef = {} for hemisphere in ['NH', 'SH']: - dsTimeSeries[hemisphere] = self._compute_area_vol(ds, hemisphere) - self.logger.info(' Make {} plots...'.format(hemisphere)) for variableName in ['iceArea', 'iceVolume']: @@ -296,6 +311,10 @@ def run_task(self): # {{{ plotVars[key] = (norm[variableName] * dsTimeSeries[hemisphere][variableName]) + if self.refConfig is not None: + plotVarsRef[key] = norm[variableName] * \ + dsTimeSeriesRef[hemisphere][variableName] + prefix = '{}/{}{}_{}'.format(self.plotsDirectory, variableName, hemisphere, @@ -304,23 +323,19 @@ def run_task(self): # {{{ figureNameStd[key] = '{}.png'.format(prefix) figureNamePolar[key] = '{}_polar.png'.format(prefix) - title[key] = '{} ({}) \n {} (black)'.format( - plotTitles[variableName], hemisphere, mainRunName) + title[key] = '{} ({})'.format(plotTitles[variableName], + hemisphere) if compareWithObservations: key = (hemisphere, 'iceArea') - title[key] = '{}\nSSM/I observations, annual cycle ' \ - '(blue)'.format(title[key]) + obsLegend[key] = 'SSM/I observations, annual cycle ' if hemisphere == 'NH': key = (hemisphere, 'iceVolume') - title[key] = \ - '{}\nPIOMAS, annual cycle (blue)'.format(title[key]) + obsLegend[key] = 'PIOMAS, annual cycle (blue)' if preprocessedReferenceRunName != 'None': for variableName in ['iceArea', 'iceVolume']: key = (hemisphere, variableName) - title[key] = '{}\n {} (red)'.format( - title[key], preprocessedReferenceRunName) if compareWithObservations: dsObs = open_multifile_dataset( @@ -376,110 +391,43 @@ def run_task(self): # {{{ for variableName in ['iceArea', 'iceVolume']: key = (hemisphere, variableName) - if compareWithObservations: - if preprocessedReferenceRunName != 'None': - plotVars[key] = [plotVars[key], obs[key], - preprocessed[key]] - lineStyles = ['k-', 'b-', 'r-'] - lineWidths = [3, 1.2, 1.2] - legendText = [None, None, None] - else: - # just v1 model and obs - plotVars[key] = [plotVars[key], obs[key]] - lineStyles = ['k-', 'b-'] - lineWidths = [3, 1.2] - legendText = [None, None] - elif preprocessedReferenceRunName != 'None': - # just v1 and v0 models - plotVars[key] = [plotVars[key], preprocessed[key]] - lineStyles = ['k-', 'r-'] - lineWidths = [3, 1.2] - legendText = [None, None] - - if (compareWithObservations or - preprocessedReferenceRunName != 'None'): - # separate plots for nothern and southern hemispheres - timeseries_analysis_plot(config, plotVars[key], - movingAveragePoints, - title[key], xLabel, - units[variableName], - figureNameStd[key], - lineStyles=lineStyles, - lineWidths=lineWidths, - legendText=legendText, - titleFontSize=titleFontSize, - calendar=calendar) - filePrefix = '{}{}_{}'.format(variableName, - hemisphere, - mainRunName) - thumbnailDescription = '{} {}'.format( - hemisphere, plotTitles[variableName]) - caption = 'Running mean of {}'.format( - thumbnailDescription) - write_image_xml( - config, - filePrefix, - componentName='Sea Ice', - componentSubdirectory='sea_ice', - galleryGroup=galleryGroup, - groupLink=groupLink, - thumbnailDescription=thumbnailDescription, - imageDescription=caption, - imageCaption=caption) - - if (polarPlot): - timeseries_analysis_plot_polar( - config, - plotVars[key], - movingAveragePoints, - title[key], - figureNamePolar[key], - lineStyles=lineStyles, - lineWidths=lineWidths, - legendText=legendText, - titleFontSize=titleFontSize, - calendar=calendar) - - filePrefix = '{}{}_{}_polar'.format(variableName, - hemisphere, - mainRunName) - write_image_xml( - config, - filePrefix, - componentName='Sea Ice', - componentSubdirectory='sea_ice', - galleryGroup=galleryGroup, - groupLink=groupLink, - thumbnailDescription=thumbnailDescription, - imageDescription=caption, - imageCaption=caption) - - if (not compareWithObservations and - preprocessedReferenceRunName == 'None'): - for variableName in ['iceArea', 'iceVolume']: - # we will combine north and south onto a single graph - figureNameStd = '{}/{}.{}.png'.format(self.plotsDirectory, - mainRunName, - variableName) - figureNamePolar = \ - '{}/{}.{}_polar.png'.format(self.plotsDirectory, - mainRunName, - variableName) - title = '{}, NH (black), SH (blue)\n' \ - '{}'.format(plotTitles[variableName], mainRunName) - varList = [plotVars[('NH', variableName)], - plotVars[('SH', variableName)]] - timeseries_analysis_plot(config, varList, + dsvalues = [plotVars[key]] + legendText = [mainRunName] + lineStyles = ['k-'] + lineWidths = [3] + if compareWithObservations and key in obsLegend.keys(): + dsvalues.append(obs[key]) + legendText.append(obsLegend[key]) + lineStyles.append('b-') + lineWidths.append(1.2) + if preprocessedReferenceRunName != 'None': + dsvalues.append(preprocessed[key]) + legendText.append(preprocessedReferenceRunName) + lineStyles.append('r-') + lineWidths.append(1.2) + + if self.refConfig is not None: + dsvalues.append(plotVarsRef[key]) + legendText.append(refRunName) + lineStyles.append('g-') + lineWidths.append(1.2) + + # separate plots for nothern and southern hemispheres + timeseries_analysis_plot(config, dsvalues, movingAveragePoints, - title, xLabel, units[variableName], - figureNameStd, - lineStyles=['k-', 'b-'], - lineWidths=[2, 2], - legendText=[None, None], + title[key], xLabel, + units[variableName], + figureNameStd[key], + lineStyles=lineStyles, + lineWidths=lineWidths, + legendText=legendText, titleFontSize=titleFontSize, calendar=calendar) - filePrefix = '{}.{}'.format(mainRunName, variableName) - thumbnailDescription = plotTitles[variableName] + filePrefix = '{}{}_{}'.format(variableName, + hemisphere, + mainRunName) + thumbnailDescription = '{} {}'.format( + hemisphere, plotTitles[variableName]) caption = 'Running mean of {}'.format( thumbnailDescription) write_image_xml( @@ -492,17 +440,23 @@ def run_task(self): # {{{ thumbnailDescription=thumbnailDescription, imageDescription=caption, imageCaption=caption) - if (polarPlot): - timeseries_analysis_plot_polar(config, varList, - movingAveragePoints, - title, figureNamePolar, - lineStyles=['k-', 'b-'], - lineWidths=[2, 2], - legendText=[None, None], - titleFontSize=titleFontSize, - calendar=calendar) - filePrefix = '{}.{}_polar'.format(mainRunName, - variableName) + + if (polarPlot): + timeseries_analysis_plot_polar( + config, + dsvalues, + movingAveragePoints, + title[key], + figureNamePolar[key], + lineStyles=lineStyles, + lineWidths=lineWidths, + legendText=legendText, + titleFontSize=titleFontSize, + calendar=calendar) + + filePrefix = '{}{}_{}_polar'.format(variableName, + hemisphere, + mainRunName) write_image_xml( config, filePrefix, @@ -513,7 +467,6 @@ def run_task(self): # {{{ thumbnailDescription=thumbnailDescription, imageDescription=caption, imageCaption=caption) - # }}} def _replicate_cycle(self, ds, dsToReplicate, calendar): # {{{ @@ -586,36 +539,72 @@ def _replicate_cycle(self, ds, dsToReplicate, calendar): # {{{ return dsShift # }}} - def _compute_area_vol(self, ds, hemisphere): # {{{ + def _compute_area_vol(self): # {{{ ''' Compute part of the time series of sea ice volume and area, given time indices to process. ''' - if hemisphere == 'NH': - mask = self.dsMesh.latCell > 0 + outFileNames = {} + allExist = True + for hemisphere in ['NH', 'SH']: + baseDirectory = build_config_full_path( + self.config, 'output', 'timeSeriesSubdirectory') + + make_directories(baseDirectory) + + outFileName = '{}/seaIceAreaVol{}.nc'.format(baseDirectory, + hemisphere) + outFileNames[hemisphere] = outFileName + if not os.path.exists(outFileName): + allExist = False + + dsTimeSeries = {} + if allExist: + for hemisphere in ['NH', 'SH']: + dsTimeSeries[hemisphere] = xr.open_dataset( + outFileNames[hemisphere]) else: - mask = self.dsMesh.latCell < 0 - ds = ds.where(mask) - - dsAreaSum = (ds*self.dsMesh.areaCell).sum('nCells') - dsAreaSum = dsAreaSum.rename( - {'timeMonthly_avg_iceAreaCell': 'iceArea', - 'timeMonthly_avg_iceVolumeCell': 'iceVolume'}) - dsAreaSum['iceThickness'] = (dsAreaSum.iceVolume / - self.dsMesh.areaCell.sum('nCells')) - - dsAreaSum['iceArea'].attrs['units'] = 'm$^2$' - dsAreaSum['iceArea'].attrs['description'] = \ - 'Total {} sea ice area'.format(hemisphere) - dsAreaSum['iceVolume'].attrs['units'] = 'm$^3$' - dsAreaSum['iceVolume'].attrs['description'] = \ - 'Total {} sea ice volume'.format(hemisphere) - dsAreaSum['iceThickness'].attrs['units'] = 'm' - dsAreaSum['iceThickness'].attrs['description'] = \ - 'Mean {} sea ice volume'.format(hemisphere) - - return dsAreaSum # }}} + dsMesh = xr.open_dataset(self.restartFileName) + dsMesh = subset_variables(dsMesh, + variableList=['latCell', 'areaCell']) + # Load data + ds = open_mpas_dataset( + fileName=self.inputFile, + calendar=self.calendar, + variableList=self.variableList, + startDate=self.startDate, + endDate=self.endDate) + + for hemisphere in ['NH', 'SH']: + + if hemisphere == 'NH': + mask = dsMesh.latCell > 0 + else: + mask = dsMesh.latCell < 0 + + dsAreaSum = (ds.where(mask)*dsMesh.areaCell).sum('nCells') + dsAreaSum = dsAreaSum.rename( + {'timeMonthly_avg_iceAreaCell': 'iceArea', + 'timeMonthly_avg_iceVolumeCell': 'iceVolume'}) + dsAreaSum['iceThickness'] = (dsAreaSum.iceVolume / + dsMesh.areaCell.sum('nCells')) + + dsAreaSum['iceArea'].attrs['units'] = 'm$^2$' + dsAreaSum['iceArea'].attrs['description'] = \ + 'Total {} sea ice area'.format(hemisphere) + dsAreaSum['iceVolume'].attrs['units'] = 'm$^3$' + dsAreaSum['iceVolume'].attrs['description'] = \ + 'Total {} sea ice volume'.format(hemisphere) + dsAreaSum['iceThickness'].attrs['units'] = 'm' + dsAreaSum['iceThickness'].attrs['description'] = \ + 'Mean {} sea ice volume'.format(hemisphere) + + dsTimeSeries[hemisphere] = dsAreaSum + + write_netcdf(dsAreaSum, outFileNames[hemisphere]) + + return dsTimeSeries # }}} # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/shared/analysis_task.py b/mpas_analysis/shared/analysis_task.py index c60f93a3e..f03625ecd 100644 --- a/mpas_analysis/shared/analysis_task.py +++ b/mpas_analysis/shared/analysis_task.py @@ -293,15 +293,6 @@ def run(self, writeLogFile=True): # {{{ Xylar Asay-Davis ''' # redirect output to a log file - logsDirectory = build_config_full_path(self.config, 'output', - 'logsSubdirectory') - - configFileName = '{}/configs/config.{}'.format(logsDirectory, - self.fullTaskName) - configFile = open(configFileName, 'w') - self.config.write(configFile) - configFile.close() - if writeLogFile: self.logger = logging.getLogger(self.fullTaskName) handler = logging.FileHandler(self._logFileName) diff --git a/mpas_analysis/shared/climatology/__init__.py b/mpas_analysis/shared/climatology/__init__.py index 279616e8d..3be6e13d9 100644 --- a/mpas_analysis/shared/climatology/__init__.py +++ b/mpas_analysis/shared/climatology/__init__.py @@ -1,6 +1,10 @@ from .climatology import get_remapper, \ compute_monthly_climatology, compute_climatology, \ - add_years_months_days_in_month, remap_and_write_climatology + add_years_months_days_in_month, remap_and_write_climatology, \ + get_unmasked_mpas_climatology_directory, \ + get_unmasked_mpas_climatology_file_name, \ + get_masked_mpas_climatology_file_name, \ + get_remapped_mpas_climatology_file_name from .mpas_climatology_task import MpasClimatologyTask from .remap_mpas_climatology_subtask import RemapMpasClimatologySubtask diff --git a/mpas_analysis/shared/climatology/climatology.py b/mpas_analysis/shared/climatology/climatology.py index 9601e3f39..a33469422 100644 --- a/mpas_analysis/shared/climatology/climatology.py +++ b/mpas_analysis/shared/climatology/climatology.py @@ -23,6 +23,7 @@ from ..interpolation import Remapper from ..grid import LatLonGridDescriptor, ProjectionGridDescriptor +from .comparison_descriptors import get_comparison_descriptor def get_remapper(config, sourceDescriptor, comparisonDescriptor, @@ -105,72 +106,6 @@ def get_remapper(config, sourceDescriptor, comparisonDescriptor, return remapper # }}} -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 remapping. - - 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. - - monthNames : str - A string identifying the months in a seasonal climatology (e.g. 'JFM') - - remapper : ``Remapper`` object - A remapper that used to remap files or data sets from the - observation grid to a comparison grid - - Returns - ------- - climatologyFileName : str - The absolute path to a file where the climatology should be stored - before remapping. - - remappedFileName : str - The absolute path to a file where the climatology should be stored - after remapping. - - Authors - ------- - Xylar Asay-Davis - """ - - obsSection = '{}Observations'.format(componentName) - - climatologyDirectory = build_config_full_path( - config=config, section='output', - relativePathOption='climatologySubdirectory', - relativePathSection=obsSection) - - remappedDirectory = build_config_full_path( - config=config, section='output', - relativePathOption='remappedClimSubdirectory', - relativePathSection=obsSection) - - obsGridName = remapper.sourceDescriptor.meshName - comparisonGridName = remapper.destinationDescriptor.meshName - - climatologyFileName = '{}/{}_{}_{}.nc'.format( - climatologyDirectory, fieldName, obsGridName, monthNames) - remappedFileName = '{}/{}_{}_to_{}_{}.nc'.format( - remappedDirectory, fieldName, obsGridName, comparisonGridName, - monthNames) - - make_directories(climatologyDirectory) - - make_directories(remappedDirectory) - - return (climatologyFileName, remappedFileName) # }}} - - def compute_monthly_climatology(ds, calendar=None, maskVaries=True): # {{{ """ Compute monthly climatologies from a data set. The mean is weighted but @@ -419,6 +354,218 @@ def remap_and_write_climatology(config, climatologyDataSet, return remappedClimatology # }}} +def get_unmasked_mpas_climatology_directory(config): # {{{ + """ + Get the directory for an unmasked MPAS climatology produced by ncclimo, + making the directory if it doesn't already exist + + Parameters + ---------- + config : ``MpasAnalysisConfigParser`` + configuration options + + Authors + ------- + Xylar Asay-Davis + """ + + climatologyBaseDirectory = build_config_full_path( + config, 'output', 'mpasClimatologySubdirectory') + + mpasMeshName = config.get('input', 'mpasMeshName') + + directory = '{}/unmasked_{}'.format(climatologyBaseDirectory, + mpasMeshName) + + make_directories(directory) + return directory # }}} + + +def get_unmasked_mpas_climatology_file_name(config, season, componentName): + # {{{ + """ + Get the file name for an unmasked MPAS climatology produced by ncclimo + + Parameters + ---------- + config : ``MpasAnalysisConfigParser`` + configuration options + + season : str + One of the seasons in ``constants.monthDictionary`` + + componentName : {'ocean', 'seaIce'} + The MPAS component for which the climatology is being computed + + Authors + ------- + Xylar Asay-Davis + """ + + startYear = config.getint('climatology', 'startYear') + endYear = config.getint('climatology', 'endYear') + + if componentName == 'ocean': + ncclimoModel = 'mpaso' + elif componentName == 'seaIce': + ncclimoModel = 'mpascice' + else: + raise ValueError('component {} is not supported by ncclimo.\n' + 'Check with Charlie Zender and Xylar Asay-Davis\n' + 'about getting it added'.format(componentName)) + + directory = get_unmasked_mpas_climatology_directory(config) + + make_directories(directory) + monthValues = sorted(constants.monthDictionary[season]) + startMonth = monthValues[0] + endMonth = monthValues[-1] + + suffix = '{:04d}{:02d}_{:04d}{:02d}_climo'.format( + startYear, startMonth, endYear, endMonth) + + if season in constants.abrevMonthNames: + season = '{:02d}'.format(monthValues[0]) + fileName = '{}/{}_{}_{}.nc'.format(directory, ncclimoModel, + season, suffix) + return fileName # }}} + + +def get_masked_mpas_climatology_file_name(config, season, componentName, + climatologyName): # {{{ + """ + Get the file name for a masked MPAS climatology + + Parameters + ---------- + config : ``MpasAnalysisConfigParser`` + Configuration options + + season : str + One of the seasons in ``constants.monthDictionary`` + + componentName : {'ocean', 'seaIce'} + The MPAS component for which the climatology is being computed + + climatologyName : str + The name of the climatology (typically the name of a field to mask + and later remap) + + Authors + ------- + Xylar Asay-Davis + """ + + startYear = config.getint('climatology', 'startYear') + endYear = config.getint('climatology', 'endYear') + mpasMeshName = config.get('input', 'mpasMeshName') + + if componentName == 'ocean': + ncclimoModel = 'mpaso' + elif componentName == 'seaIce': + ncclimoModel = 'mpascice' + else: + raise ValueError('component {} is not supported by ncclimo.\n' + 'Check with Charlie Zender and Xylar Asay-Davis\n' + 'about getting it added'.format(componentName)) + + climatologyBaseDirectory = build_config_full_path( + config, 'output', 'mpasClimatologySubdirectory') + + stageDirectory = '{}/masked'.format(climatologyBaseDirectory) + + directory = '{}/{}_{}'.format( + stageDirectory, climatologyName, + mpasMeshName) + + make_directories(directory) + + monthValues = sorted(constants.monthDictionary[season]) + startMonth = monthValues[0] + endMonth = monthValues[-1] + + suffix = '{:04d}{:02d}_{:04d}{:02d}_climo'.format( + startYear, startMonth, endYear, endMonth) + + if season in constants.abrevMonthNames: + season = '{:02d}'.format(monthValues[0]) + fileName = '{}/{}_{}_{}.nc'.format( + directory, ncclimoModel, season, suffix) + + return fileName # }}} + + +def get_remapped_mpas_climatology_file_name(config, season, componentName, + climatologyName, + comparisonGridName): # {{{ + """ + Get the file name for a masked MPAS climatology + + Parameters + ---------- + config : ``MpasAnalysisConfigParser`` + Configuration options + + season : str + One of the seasons in ``constants.monthDictionary`` + + componentName : {'ocean', 'seaIce'} + The MPAS component for which the climatology is being computed + + climatologyName : str + The name of the climatology (typically the name of a field to mask + and later remap) + + comparisonGridName : {'latlon', 'antarctic'} + The name of the comparison grid to use for remapping + + Authors + ------- + Xylar Asay-Davis + """ + + startYear = config.getint('climatology', 'startYear') + endYear = config.getint('climatology', 'endYear') + mpasMeshName = config.get('input', 'mpasMeshName') + + if componentName == 'ocean': + ncclimoModel = 'mpaso' + elif componentName == 'seaIce': + ncclimoModel = 'mpascice' + else: + raise ValueError('component {} is not supported by ncclimo.\n' + 'Check with Charlie Zender and Xylar Asay-Davis\n' + 'about getting it added'.format(componentName)) + + climatologyBaseDirectory = build_config_full_path( + config, 'output', 'mpasClimatologySubdirectory') + + comparisonDescriptor = get_comparison_descriptor(config, + comparisonGridName) + comparisonFullMeshName = comparisonDescriptor.meshName + + stageDirectory = '{}/remapped'.format(climatologyBaseDirectory) + + directory = '{}/{}_{}_to_{}'.format(stageDirectory, climatologyName, + mpasMeshName, comparisonFullMeshName) + + make_directories(directory) + + monthValues = sorted(constants.monthDictionary[season]) + startMonth = monthValues[0] + endMonth = monthValues[-1] + + suffix = '{:04d}{:02d}_{:04d}{:02d}_climo'.format( + startYear, startMonth, endYear, endMonth) + + if season in constants.abrevMonthNames: + season = '{:02d}'.format(monthValues[0]) + fileName = '{}/{}_{}_{}.nc'.format( + directory, ncclimoModel, season, suffix) + + return fileName # }}} + + def _compute_masked_mean(ds, maskVaries): # {{{ ''' Compute the time average of data set, masked out where the variables in ds diff --git a/mpas_analysis/shared/climatology/mpas_climatology_task.py b/mpas_analysis/shared/climatology/mpas_climatology_task.py index b740cdcad..fa5d05824 100644 --- a/mpas_analysis/shared/climatology/mpas_climatology_task.py +++ b/mpas_analysis/shared/climatology/mpas_climatology_task.py @@ -9,9 +9,8 @@ from ..analysis_task import AnalysisTask -from ..constants import constants - -from ..io.utility import build_config_full_path, make_directories +from .climatology import get_unmasked_mpas_climatology_directory, \ + get_unmasked_mpas_climatology_file_name class MpasClimatologyTask(AnalysisTask): # {{{ @@ -48,7 +47,7 @@ class MpasClimatologyTask(AnalysisTask): # {{{ Xylar Asay-Davis ''' - def __init__(self, config, componentName): # {{{ + def __init__(self, config, componentName, taskName=None): # {{{ ''' Construct the analysis task and adds it as a subtask of the ``parentTask``. @@ -62,6 +61,9 @@ def __init__(self, config, componentName): # {{{ The name of the component (same as the folder where the task resides) + taskName : str, optional + the name of the task, defaults to mpasClimatology + Authors ------- Xylar Asay-Davis @@ -80,8 +82,9 @@ def __init__(self, config, componentName): # {{{ 'Check with Charlie Zender and Xylar Asay-Davis\n' 'about getting it added'.format(componentName)) - suffix = componentName[0].upper() + componentName[1:] - taskName = 'mpasClimatology{}'.format(suffix) + if taskName is None: + suffix = componentName[0].upper() + componentName[1:] + taskName = 'mpasClimatology{}'.format(suffix) # call the constructor from the base class (AnalysisTask) super(MpasClimatologyTask, self).__init__( @@ -122,8 +125,6 @@ def add_variables(self, variableList, seasons=None): # {{{ if season not in self.seasons: self.seasons.append(season) - self._setup_file_names() - # }}} def setup_and_check(self): # {{{ @@ -196,8 +197,9 @@ def run_task(self): # {{{ allExist = True for season in seasonsToCheck: - climatologyFileName, climatologyDirectory = \ - self.get_file_name(season, returnDir=True) + climatologyFileName = self.get_file_name(season) + climatologyDirectory = get_unmasked_mpas_climatology_directory( + self.config) if not os.path.exists(climatologyFileName): allExist = False @@ -205,8 +207,7 @@ def run_task(self): # {{{ if allExist: # make sure all the necessary variables are also present - ds = xarray.open_dataset(self.get_file_name(seasonsToCheck[0], - returnDir=False)) + ds = xarray.open_dataset(self.get_file_name(seasonsToCheck[0])) for variableName in self.variableList: if variableName not in ds.variables: @@ -220,23 +221,16 @@ def run_task(self): # {{{ # }}} - def get_file_name(self, season, returnDir=False): # {{{ + def get_file_name(self, season): # {{{ """ - 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 remapping. + + Returns the full path for MPAS climatology file produced by ncclimo. Parameters ---------- season : str One of the seasons in ``constants.monthDictionary`` - mpasMeshName : str - The name of the MPAS mesh - - returnDir : bool, optional - Return the directory as well - Returns ------- fileName : str @@ -247,13 +241,10 @@ def get_file_name(self, season, returnDir=False): # {{{ Xylar Asay-Davis """ - fileName = self._outputFiles[season] + return get_unmasked_mpas_climatology_file_name(self.config, season, + self.componentName) - if returnDir: - directory = self._outputDirs[season] - return fileName, directory - else: - return fileName # }}} + # }}} def _update_climatology_bounds_from_file_names(self): # {{{ """ @@ -313,46 +304,6 @@ def _update_climatology_bounds_from_file_names(self): # {{{ # }}} - def _setup_file_names(self): # {{{ - """ - Create a dictionary of file names and directories for this climatology - - Authors - ------- - Xylar Asay-Davis - """ - - config = self.config - climatologyBaseDirectory = build_config_full_path( - config, 'output', 'mpasClimatologySubdirectory') - - mpasMeshName = config.get('input', 'mpasMeshName') - - self._outputDirs = {} - self._outputFiles = {} - - directory = '{}/unmasked_{}'.format(climatologyBaseDirectory, - mpasMeshName) - - make_directories(directory) - for season in self.seasons: - monthValues = sorted(constants.monthDictionary[season]) - startMonth = monthValues[0] - endMonth = monthValues[-1] - - suffix = '{:04d}{:02d}_{:04d}{:02d}_climo'.format( - self.startYear, startMonth, self.endYear, endMonth) - - if season in constants.abrevMonthNames: - season = '{:02d}'.format(monthValues[0]) - fileName = '{}/{}_{}_{}.nc'.format(directory, self.ncclimoModel, - season, suffix) - - self._outputDirs[season] = directory - self._outputFiles[season] = fileName - - # }}} - def _compute_climatologies_with_ncclimo(self, inDirectory, outDirectory, remapper=None, remappedDirectory=None): # {{{ diff --git a/mpas_analysis/shared/climatology/remap_mpas_climatology_subtask.py b/mpas_analysis/shared/climatology/remap_mpas_climatology_subtask.py index 70c447e6a..3e9d55067 100644 --- a/mpas_analysis/shared/climatology/remap_mpas_climatology_subtask.py +++ b/mpas_analysis/shared/climatology/remap_mpas_climatology_subtask.py @@ -13,7 +13,8 @@ from ..io.utility import build_config_full_path, make_directories from ..io import write_netcdf -from .climatology import get_remapper +from .climatology import get_remapper, get_masked_mpas_climatology_file_name, \ + get_remapped_mpas_climatology_file_name from .comparison_descriptors import get_comparison_descriptor from ..grid import MpasMeshDescriptor @@ -61,7 +62,7 @@ class RemapMpasClimatologySubtask(AnalysisTask): # {{{ def __init__(self, mpasClimatologyTask, parentTask, climatologyName, variableList, seasons, comparisonGridNames=['latlon'], - iselValues=None): + iselValues=None, subtaskName='remapMpasClimatology'): # {{{ ''' Construct the analysis task and adds it as a subtask of the @@ -97,6 +98,9 @@ def __init__(self, mpasClimatologyTask, parentTask, climatologyName, A dictionary of dimensions and indices (or ``None``) used to extract a slice of the MPAS field(s). + subtaskName : str, optional + The name of the subtask + Authors ------- Xylar Asay-Davis @@ -105,9 +109,9 @@ def __init__(self, mpasClimatologyTask, parentTask, climatologyName, # call the constructor from the base class (AnalysisTask) super(RemapMpasClimatologySubtask, self).__init__( - config=parentTask.config, + config=mpasClimatologyTask.config, taskName=parentTask.taskName, - subtaskName='remapMpasClimatology', + subtaskName=subtaskName, componentName=parentTask.componentName, tags=tags) @@ -166,8 +170,6 @@ def setup_and_check(self): # {{{ # we're sure this subtask is supposed to run self.mpasClimatologyTask.add_variables(self.variableList, self.seasons) - self._setup_file_names() - # make the mapping directory, because doing so within each process # seems to be giving ESMF_RegridWeightGen some trouble mappingSubdirectory = build_config_full_path(self.config, 'output', @@ -204,11 +206,11 @@ def run_task(self): # {{{ for season in self.seasons: - maskedClimatologyFileName = self.get_file_name( - season, 'masked', comparisonGridName) + maskedClimatologyFileName = self.get_masked_file_name( + season) - remappedFileName = self.get_file_name( - season, 'remapped', comparisonGridName) + remappedFileName = self.get_remapped_file_name( + season, comparisonGridName) if not os.path.exists(remappedFileName): self._remap(inFileName=maskedClimatologyFileName, @@ -217,7 +219,7 @@ def run_task(self): # {{{ comparisonGridName=comparisonGridName) # }}} - def get_file_name(self, season, stage, comparisonGridName=None): # {{{ + def get_masked_file_name(self, season): # {{{ """ Given config options, the name of a field and a string identifying the months in a seasonal climatology, returns the full path for MPAS @@ -228,10 +230,35 @@ def get_file_name(self, season, stage, comparisonGridName=None): # {{{ season : str One of the seasons in ``constants.monthDictionary`` - stage : {'masked', 'remapped'} - The stage of the masking and remapping process + Returns + ------- + fileName : str + The path to the climatology file for the specified season. - comparisonGridName : {'latlon', 'antarctic'}, optional + Authors + ------- + Xylar Asay-Davis + """ + + fileName = get_masked_mpas_climatology_file_name(self.config, + season, + self.componentName, + self.climatologyName) + + return fileName # }}} + + def get_remapped_file_name(self, season, comparisonGridName): # {{{ + """ + 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 remapping. + + Parameters + ---------- + season : str + One of the seasons in ``constants.monthDictionary`` + + comparisonGridName : {'latlon', 'antarctic'} The name of the comparison grid to use for remapping. Returns @@ -244,11 +271,9 @@ def get_file_name(self, season, stage, comparisonGridName=None): # {{{ Xylar Asay-Davis """ - if stage == 'remapped': - key = (season, stage, comparisonGridName) - else: - key = (season, stage) - fileName = self._outputFiles[key] + fileName = get_remapped_mpas_climatology_file_name( + self.config, season, self.componentName, self.climatologyName, + comparisonGridName) return fileName # }}} @@ -423,7 +448,7 @@ def _mask_climatologies(self, season, dsMask): # {{{ climatologyFileName = self.mpasClimatologyTask.get_file_name(season) - maskedClimatologyFileName = self.get_file_name(season, 'masked') + maskedClimatologyFileName = self.get_masked_file_name(season) if not os.path.exists(maskedClimatologyFileName): # slice and mask the data set diff --git a/mpas_analysis/shared/html/pages.py b/mpas_analysis/shared/html/pages.py index f01e97c2a..7c3640228 100644 --- a/mpas_analysis/shared/html/pages.py +++ b/mpas_analysis/shared/html/pages.py @@ -11,14 +11,14 @@ from ..io.utility import build_config_full_path -def generate_html(config, analyses): # {{{ +def generate_html(config, analyses, refConfig=None): # {{{ """ Generates webpages for diplaying the plots from each analysis task Parameters ---------- - config : ``MpasAnalysisConfigParser`` object - contains config options + config : ``MpasAnalysisConfigParser`` + Config options analysis : ``OrderedDict`` of ``AnalysisTask`` objects the analysis tasks that generated the plots to include in the webpages. @@ -26,6 +26,10 @@ def generate_html(config, analyses): # {{{ the list of files to include on the webpage for the associated component. + refConfig : ``MpasAnalysisConfigParser``, optional + Config options for a reference run + + Authors ------- Xylar Asay-Davis @@ -36,7 +40,7 @@ def generate_html(config, analyses): # {{{ print("Generating webpage for viewing results...") - page = MainPage(config) + page = MainPage(config, refConfig) components = OrderedDict() @@ -45,7 +49,8 @@ def generate_html(config, analyses): # {{{ for analysisTask in analyses.values(): for fileName in analysisTask.xmlFileNames: try: - ComponentPage.add_image(fileName, config, components) + ComponentPage.add_image(fileName, config, components, + refConfig) except IOError: missingCount += 1 @@ -71,8 +76,11 @@ class MainPage(object): Attributes ---------- - config : ``MpasAnalysisConfigParser`` object - contains config options + config : ``MpasAnalysisConfigParser`` + Config options + + refConfig : ``MpasAnalysisConfigParser`` + Config options for a reference run pageTemplate, componentTemplate : str The contents of templates used to construct the page @@ -85,14 +93,17 @@ class MainPage(object): ------- Xylar Asay-Davis """ - def __init__(self, config): + def __init__(self, config, refConfig=None): """ Create a MainPage object, reading in the templates Parameters ---------- - config : ``MpasAnalysisConfigParser`` object - contains config options + config : ``MpasAnalysisConfigParser`` + Config options + + refConfig : ``MpasAnalysisConfigParser``, optional + Config options for a reference run Authors ------- @@ -100,6 +111,7 @@ def __init__(self, config): """ self.config = config + self.refConfig = refConfig # get template text fileName = \ @@ -155,6 +167,12 @@ def generate(self): """ runName = self.config.get('runs', 'mainRunName') + if self.refConfig is None: + refRunText = '' + else: + refRunText = '
Ref: {}'.format( + self.refConfig.get('runs', 'mainRunName')) + componentsText = '' for componentName, componentDict in self.components.items(): @@ -170,6 +188,7 @@ def generate(self): _replace_tempate_text(self.componentTemplate, replacements) replacements = {'@runName': runName, + '@refRunText': refRunText, '@components': componentsText} pageText = _replace_tempate_text(self.pageTemplate, replacements) @@ -214,8 +233,11 @@ class ComponentPage(object): Attributes ---------- - config : ``MpasAnalysisConfigParser`` object - contains config options + config : ``MpasAnalysisConfigParser`` + Config options + + refConfig : ``MpasAnalysisConfigParser`` + Config options for a reference run name : str The name of the component as it should appear in the list of @@ -236,14 +258,14 @@ class ComponentPage(object): ------- Xylar Asay-Davis """ - def __init__(self, config, name, subdirectory): + def __init__(self, config, name, subdirectory, refConfig=None): """ Create a ComponentPage object, reading in the templates Parameters ---------- - config : ``MpasAnalysisConfigParser`` object - contains config options + config : ``MpasAnalysisConfigParser`` + Config options name : str The name of the component as it should appear in the list of @@ -253,12 +275,16 @@ def __init__(self, config, name, subdirectory): subdirecory : str The subdirectory for the component's webpage + refConfig : ``MpasAnalysisConfigParser``, optional + Config options for a reference run + Authors ------- Xylar Asay-Davis """ self.config = config + self.refConfig = refConfig self.name = name self.subdirectory = subdirectory @@ -283,7 +309,7 @@ def __init__(self, config, name, subdirectory): self.groups = OrderedDict() @staticmethod - def add_image(xmlFileName, config, components): + def add_image(xmlFileName, config, components, refConfig=None): """ Add the image to the appropriate component. Note: this is a static method because we do not know which component to add the image to @@ -303,6 +329,9 @@ def add_image(xmlFileName, config, components): be added. ``components`` should be viewed as an input and output parameter, since it is modified by this function. + refConfig : ``MpasAnalysisConfigParser``, optional + Config options for a reference run + Authors ------- Xylar Asay-Davis @@ -327,7 +356,8 @@ def add_image(xmlFileName, config, components): if componentName not in components: components[componentName] = ComponentPage(config, componentName, - componentSubdirectory) + componentSubdirectory, + refConfig) component = components[componentName] @@ -375,6 +405,12 @@ def generate(self): """ runName = self.config.get('runs', 'mainRunName') + if self.refConfig is None: + refRunText = '' + else: + refRunText = '
Ref: {}'.format( + self.refConfig.get('runs', 'mainRunName')) + quickLinkText = '' galleriesText = '' for groupName, groupDict in self.groups.items(): @@ -385,6 +421,7 @@ def generate(self): self._generate_group_text(groupName, groupDict) replacements = {'@runName': runName, + '@refRunText': refRunText, '@componentName': self.name, '@quickLinks': quickLinkText, '@galleries': galleriesText} diff --git a/mpas_analysis/shared/html/templates/component_page.html b/mpas_analysis/shared/html/templates/component_page.html index 35cc05119..71f3f1929 100644 --- a/mpas_analysis/shared/html/templates/component_page.html +++ b/mpas_analysis/shared/html/templates/component_page.html @@ -42,7 +42,7 @@

MPAS-Analysis Diagnostics: @componentName

-

Run: @runName

+

Run: @runName @refRunText

diff --git a/mpas_analysis/shared/html/templates/main_page.html b/mpas_analysis/shared/html/templates/main_page.html index b426708f0..be3d39190 100644 --- a/mpas_analysis/shared/html/templates/main_page.html +++ b/mpas_analysis/shared/html/templates/main_page.html @@ -17,7 +17,7 @@

MPAS-Analysis Diagnostics

-

Run: @runName

+

Run: @runName @refRunText

diff --git a/mpas_analysis/shared/io/write_netcdf.py b/mpas_analysis/shared/io/write_netcdf.py index b16c57767..03072bae5 100644 --- a/mpas_analysis/shared/io/write_netcdf.py +++ b/mpas_analysis/shared/io/write_netcdf.py @@ -1,17 +1,3 @@ -''' -Functions for writing data sets - -Functions ---------- -write_netcdf - write an xarray data set to a NetCDF file using finite fill - values - -Authors -------- -Xylar Asay-Davis - -''' - from __future__ import absolute_import, division, print_function, \ unicode_literals diff --git a/mpas_analysis/shared/plot/plotting.py b/mpas_analysis/shared/plot/plotting.py index 91dbeea8d..eff883a41 100644 --- a/mpas_analysis/shared/plot/plotting.py +++ b/mpas_analysis/shared/plot/plotting.py @@ -292,10 +292,10 @@ def plot_polar_comparison( Lons, Lats, modelArray, - obsArray, + refArray, diffArray, - cmapModelObs, - clevsModelObs, + cmapModelRef, + clevsModelRef, cmapDiff, clevsDiff, fileout, @@ -304,7 +304,7 @@ def plot_polar_comparison( latmin=50.0, lon0=0, modelTitle='Model', - obsTitle='Observations', + refTitle='Observations', diffTitle='Model-Observations', cbarlabel='units', titleFontSize=None, @@ -324,17 +324,17 @@ def plot_polar_comparison( Lons, Lats : float arrays longitude and latitude arrays - modelArray, obsArray : float arrays - model and observational data sets + modelArray, refArray : float arrays + model and observational or reference run data sets diffArray : float array - difference between modelArray and obsArray + difference between modelArray and refArray - cmapModelObs : str - colormap of model and observations panel + cmapModelRef : str + colormap of model and observations or reference run panel - clevsModelObs : int array - colorbar values for model and observations panel + clevsModelRef : int array + colorbar values for model and observations or reference run panel cmapDiff : str colormap of difference (bias) panel @@ -354,8 +354,8 @@ def plot_polar_comparison( modelTitle : str, optional title of the model panel - obsTitle : str, optional - title of the observations panel + refTitle : str, optional + title of the observations or reference run panel diffTitle : str, optional title of the difference (bias) panel @@ -408,7 +408,11 @@ def do_subplot(ax, field, title, cmap, norm, levels): if dpi is None: dpi = config.getint('plot', 'dpi') - if vertical: + if refArray is None: + if figsize is None: + figsize = (8, 8.5) + subplots = [111] + elif vertical: if figsize is None: figsize = (8, 22) subplots = [311, 312, 313] @@ -428,20 +432,21 @@ def do_subplot(ax, field, title, cmap, norm, levels): fig.suptitle(title, y=0.95, **title_font) axis_font = {'size': config.get('plot', 'axisFontSize')} - normModelObs = cols.BoundaryNorm(clevsModelObs, cmapModelObs.N) + normModelRef = cols.BoundaryNorm(clevsModelRef, cmapModelRef.N) normDiff = cols.BoundaryNorm(clevsDiff, cmapDiff.N) ax = plt.subplot(subplots[0]) - do_subplot(ax=ax, field=modelArray, title=modelTitle, cmap=cmapModelObs, - norm=normModelObs, levels=clevsModelObs) + do_subplot(ax=ax, field=modelArray, title=modelTitle, cmap=cmapModelRef, + norm=normModelRef, levels=clevsModelRef) - ax = plt.subplot(subplots[1]) - do_subplot(ax=ax, field=obsArray, title=obsTitle, cmap=cmapModelObs, - norm=normModelObs, levels=clevsModelObs) + if refArray is not None: + ax = plt.subplot(subplots[1]) + do_subplot(ax=ax, field=refArray, title=refTitle, cmap=cmapModelRef, + norm=normModelRef, levels=clevsModelRef) - ax = plt.subplot(subplots[2]) - do_subplot(ax=ax, field=diffArray, title=diffTitle, cmap=cmapDiff, - norm=normDiff, levels=clevsDiff) + ax = plt.subplot(subplots[2]) + do_subplot(ax=ax, field=diffArray, title=diffTitle, cmap=cmapDiff, + norm=normDiff, levels=clevsDiff) plt.tight_layout(pad=4.) if vertical: @@ -459,16 +464,16 @@ def plot_global_comparison( Lons, Lats, modelArray, - obsArray, + refArray, diffArray, - cmapModelObs, - clevsModelObs, + cmapModelRef, + clevsModelRef, cmapDiff, clevsDiff, fileout, title=None, modelTitle='Model', - obsTitle='Observations', + refTitle='Observations', diffTitle='Model-Observations', cbarlabel='units', titleFontSize=None, @@ -487,17 +492,17 @@ def plot_global_comparison( Lons, Lats : float arrays longitude and latitude arrays - modelArray, obsArray : float arrays - model and observational data sets + modelArray, refArray : float arrays + model and observational or reference run data sets diffArray : float array - difference between modelArray and obsArray + difference between modelArray and refArray - cmapModelObs : str - colormap of model and observations panel + cmapModelRef : str + colormap of model and observations or reference run panel - clevsModelObs : int array - colorbar values for model and observations panel + clevsModelRef : int array + colorbar values for model and observations or reference run panel cmapDiff : str colormap of difference (bias) panel @@ -514,8 +519,8 @@ def plot_global_comparison( modelTitle : str, optional title of the model panel - obsTitle : str, optional - title of the observations panel + refTitle : str, optional + title of the observations or reference run panel diffTitle : str, optional title of the difference (bias) panel @@ -555,10 +560,11 @@ def plot_global_comparison( urcrnrlon=181, resolution='l') x, y = m(Lons, Lats) # compute map proj coordinates - normModelObs = cols.BoundaryNorm(clevsModelObs, cmapModelObs.N) + normModelRef = cols.BoundaryNorm(clevsModelRef, cmapModelRef.N) normDiff = cols.BoundaryNorm(clevsDiff, cmapDiff.N) - plt.subplot(3, 1, 1) + if refArray is not None: + plt.subplot(3, 1, 1) plt.title(modelTitle, y=1.06, **axis_font) m.drawcoastlines() m.fillcontinents(color='grey', lake_color='white') @@ -566,39 +572,40 @@ def plot_global_comparison( labels=[True, False, False, False]) m.drawmeridians(np.arange(-180., 180., 60.), labels=[False, False, False, True]) - cs = m.contourf(x, y, modelArray, cmap=cmapModelObs, norm=normModelObs, - levels=clevsModelObs, extend='both') + cs = m.contourf(x, y, modelArray, cmap=cmapModelRef, norm=normModelRef, + levels=clevsModelRef, extend='both') cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform', - ticks=clevsModelObs, boundaries=clevsModelObs) + ticks=clevsModelRef, boundaries=clevsModelRef) cbar.set_label(cbarlabel) - plt.subplot(3, 1, 2) - plt.title(obsTitle, y=1.06, **axis_font) - m.drawcoastlines() - m.fillcontinents(color='grey', lake_color='white') - m.drawparallels(np.arange(-80., 80., 20.), - labels=[True, False, False, False]) - m.drawmeridians(np.arange(-180., 180., 40.), - labels=[False, False, False, True]) - cs = m.contourf(x, y, obsArray, cmap=cmapModelObs, norm=normModelObs, - levels=clevsModelObs, extend='both') - cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform', - ticks=clevsModelObs, boundaries=clevsModelObs) - cbar.set_label(cbarlabel) + if refArray is not None: + plt.subplot(3, 1, 2) + plt.title(refTitle, y=1.06, **axis_font) + m.drawcoastlines() + m.fillcontinents(color='grey', lake_color='white') + m.drawparallels(np.arange(-80., 80., 20.), + labels=[True, False, False, False]) + m.drawmeridians(np.arange(-180., 180., 40.), + labels=[False, False, False, True]) + cs = m.contourf(x, y, refArray, cmap=cmapModelRef, norm=normModelRef, + levels=clevsModelRef, extend='both') + cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform', + ticks=clevsModelRef, boundaries=clevsModelRef) + cbar.set_label(cbarlabel) - plt.subplot(3, 1, 3) - plt.title(diffTitle, y=1.06, **axis_font) - m.drawcoastlines() - m.fillcontinents(color='grey', lake_color='white') - m.drawparallels(np.arange(-80., 80., 20.), - labels=[True, False, False, False]) - m.drawmeridians(np.arange(-180., 180., 40.), - labels=[False, False, False, True]) - cs = m.contourf(x, y, diffArray, cmap=cmapDiff, norm=normDiff, - levels=clevsDiff, extend='both') - cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform', - ticks=clevsDiff, boundaries=clevsModelObs) - cbar.set_label(cbarlabel) + plt.subplot(3, 1, 3) + plt.title(diffTitle, y=1.06, **axis_font) + m.drawcoastlines() + m.fillcontinents(color='grey', lake_color='white') + m.drawparallels(np.arange(-80., 80., 20.), + labels=[True, False, False, False]) + m.drawmeridians(np.arange(-180., 180., 40.), + labels=[False, False, False, True]) + cs = m.contourf(x, y, diffArray, cmap=cmapDiff, norm=normDiff, + levels=clevsDiff, extend='both') + cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform', + ticks=clevsDiff, boundaries=clevsModelRef) + cbar.set_label(cbarlabel) if (fileout is not None): plt.savefig(fileout, dpi=dpi, bbox_inches='tight', pad_inches=0.1) @@ -613,14 +620,14 @@ def plot_polar_projection_comparison( y, landMask, modelArray, - obsArray, + refArray, diffArray, fileout, colorMapSectionName, colorMapType='norm', title=None, modelTitle='Model', - obsTitle='Observations', + refTitle='Observations', diffTitle='Model-Observations', cbarlabel='units', titleFontSize=None, @@ -641,13 +648,13 @@ def plot_polar_projection_comparison( 1D x and y arrays defining the projection grid landMask : numpy ndarrays - model and observational data sets + model and observational or reference run data sets - modelArray, obsArray : numpy ndarrays - model and observational data sets + modelArray, refArray : numpy ndarrays + model and observational or reference run data sets diffArray : float array - difference between modelArray and obsArray + difference between modelArray and refArray fileout : str the file name to be written @@ -680,8 +687,8 @@ def plot_polar_projection_comparison( modelTitle : str, optional title of the model panel - obsTitle : str, optional - title of the observations panel + refTitle : str, optional + title of the observations or reference run panel diffTitle : str, optional title of the difference (bias) panel @@ -737,7 +744,11 @@ def plot_panel(ax, title, array, cmap, norm, ticks): if dpi is None: dpi = config.getint('plot', 'dpi') - if vertical: + if refArray is None: + if figsize is None: + figsize = (8, 7.5) + subplots = [111] + elif vertical: if figsize is None: figsize = (8, 22) subplots = [311, 312, 313] @@ -747,7 +758,7 @@ def plot_panel(ax, title, array, cmap, norm, ticks): subplots = [131, 132, 133] if colorMapType == 'norm': - (cmapModelObs, normModelObs) = _setup_colormap_and_norm( + (cmapModelRef, normModelRef) = _setup_colormap_and_norm( config, colorMapSectionName, suffix='Result') (cmapDiff, normDiff) = _setup_colormap_and_norm( config, colorMapSectionName, suffix='Difference') @@ -767,12 +778,12 @@ def plot_panel(ax, title, array, cmap, norm, ticks): elif colorMapType == 'indexed': - (cmapModelObs, colorbarTicksResult) = setup_colormap( + (cmapModelRef, colorbarTicksResult) = setup_colormap( config, colorMapSectionName, suffix='Result') (cmapDiff, colorbarTicksDifference) = setup_colormap( config, colorMapSectionName, suffix='Difference') - normModelObs = cols.BoundaryNorm(colorbarTicksResult, cmapModelObs.N) + normModelRef = cols.BoundaryNorm(colorbarTicksResult, cmapModelRef.N) normDiff = cols.BoundaryNorm(colorbarTicksDifference, cmapDiff.N) else: raise ValueError('colorMapType must be one of {norm, indexed}') @@ -797,16 +808,17 @@ def plot_panel(ax, title, array, cmap, norm, ticks): yCenter = 0.5*(y[1:] + y[0:-1]) ax = plt.subplot(subplots[0]) - plot_panel(ax, modelTitle, modelArray, cmapModelObs, normModelObs, + plot_panel(ax, modelTitle, modelArray, cmapModelRef, normModelRef, colorbarTicksResult) - ax = plt.subplot(subplots[1]) - plot_panel(ax, obsTitle, obsArray, cmapModelObs, normModelObs, - colorbarTicksResult) + if refArray is not None: + ax = plt.subplot(subplots[1]) + plot_panel(ax, refTitle, refArray, cmapModelRef, normModelRef, + colorbarTicksResult) - ax = plt.subplot(subplots[2]) - plot_panel(ax, diffTitle, diffArray, cmapDiff, normDiff, - colorbarTicksDifference) + ax = plt.subplot(subplots[2]) + plot_panel(ax, diffTitle, diffArray, cmapDiff, normDiff, + colorbarTicksDifference) if (fileout is not None): plt.savefig(fileout, dpi=dpi, bbox_inches='tight', pad_inches=0.1) diff --git a/mpas_analysis/shared/time_series/__init__.py b/mpas_analysis/shared/time_series/__init__.py index 7106aab94..787036001 100644 --- a/mpas_analysis/shared/time_series/__init__.py +++ b/mpas_analysis/shared/time_series/__init__.py @@ -3,4 +3,3 @@ from .anomaly import compute_moving_avg_anomaly_from_start from .moving_average import compute_moving_avg - diff --git a/mpas_analysis/test/test_climatology.py b/mpas_analysis/test/test_climatology.py index 087ccd7ca..bb3b1f3c8 100644 --- a/mpas_analysis/test/test_climatology.py +++ b/mpas_analysis/test/test_climatology.py @@ -51,8 +51,6 @@ def setup_config(self, autocloseFileLimitFraction=0.5, config.set('output', 'baseDirectory', self.test_dir) config.set('output', 'mappingSubdirectory', '.') config.set('output', 'mpasClimatologySubdirectory', 'clim/mpas') - config.set('output', 'mpasRemappedClimSubdirectory', - 'clim/mpas/remap') config.add_section('climatology') config.set('climatology', 'startYear', '2') diff --git a/mpas_analysis/test/test_mpas_climatology_task.py b/mpas_analysis/test/test_mpas_climatology_task.py index 8ba50800b..e336043cb 100644 --- a/mpas_analysis/test/test_mpas_climatology_task.py +++ b/mpas_analysis/test/test_mpas_climatology_task.py @@ -169,23 +169,23 @@ def test_subtask_run_analysis(self): remapSubtask.run(writeLogFile=False) for season in remapSubtask.seasons: - for stage in ['masked', 'remapped']: - fileName = remapSubtask.get_file_name( - season=season, stage=stage, - comparisonGridName='latlon') - assert(os.path.exists(fileName)) + fileName = remapSubtask.get_masked_file_name(season=season) + assert(os.path.exists(fileName)) + + fileName = remapSubtask.get_remapped_file_name( + season=season, comparisonGridName='latlon') + assert(os.path.exists(fileName)) def test_subtask_get_file_name(self): mpasClimatologyTask = self.setup_task() variableList, seasons = self.add_variables(mpasClimatologyTask) remapSubtask = self.setup_subtask(mpasClimatologyTask) - fileName = remapSubtask.get_file_name(season='JFM', stage='masked', - comparisonGridName='latlon') + fileName = remapSubtask.get_masked_file_name(season='JFM') assert(fileName == '{}/clim/mpas/masked/ssh_oQU240/' 'mpaso_JFM_000201_000203_climo.nc'.format(str(self.test_dir))) - fileName = remapSubtask.get_file_name(season='JFM', stage='remapped', - comparisonGridName='latlon') + fileName = remapSubtask.get_remapped_file_name( + season='JFM', comparisonGridName='latlon') assert(fileName == '{}/clim/mpas/remapped/ssh_oQU240_to_0.5x0.5degree/' 'mpaso_JFM_000201_000203_climo.nc'.format(str(self.test_dir))) diff --git a/run_mpas_analysis b/run_mpas_analysis index d3a4fac02..0472d9a1b 100755 --- a/run_mpas_analysis +++ b/run_mpas_analysis @@ -31,7 +31,7 @@ from mpas_analysis.shared.html import generate_html from mpas_analysis.shared import AnalysisTask -def build_analysis_list(config): # {{{ +def build_analysis_list(config, refConfig): # {{{ """ Build a list of analysis tasks. New tasks should be added here, following the approach used for existing analysis tasks. @@ -41,6 +41,10 @@ def build_analysis_list(config): # {{{ config : ``MpasAnalysisConfigParser`` object contains config options + refConfig : ``MpasAnalysisConfigParser`` object + contains config options for a reference run, or ``None`` if no config + file for a reference run was specified + Returns ------- analyses : list of ``AnalysisTask`` objects @@ -62,7 +66,6 @@ def build_analysis_list(config): # {{{ from mpas_analysis.shared.climatology import MpasClimatologyTask from mpas_analysis.shared.time_series import MpasTimeSeriesTask - # analyses will be a list of analysis classes analyses = [] # Ocean Analyses @@ -73,51 +76,63 @@ def build_analysis_list(config): # {{{ oceanIndexTask = MpasTimeSeriesTask(config=config, componentName='ocean', section='index') + analyses.append(oceanClimatolgyTask) - analyses.append(ocean.ClimatologyMapMLD(config, oceanClimatolgyTask)) - analyses.append(ocean.ClimatologyMapSST(config, oceanClimatolgyTask)) - analyses.append(ocean.ClimatologyMapSSS(config, oceanClimatolgyTask)) - analyses.append(ocean.ClimatologyMapSoseTemperature(config, - oceanClimatolgyTask)) - analyses.append(ocean.ClimatologyMapSoseSalinity(config, - oceanClimatolgyTask)) + + analyses.append(ocean.ClimatologyMapMLD(config, oceanClimatolgyTask, + refConfig)) + analyses.append(ocean.ClimatologyMapSST(config, oceanClimatolgyTask, + refConfig)) + analyses.append(ocean.ClimatologyMapSSS(config, oceanClimatolgyTask, + refConfig)) + analyses.append(ocean.ClimatologyMapSoseTemperature( + config, oceanClimatolgyTask, refConfig)) + analyses.append(ocean.ClimatologyMapSoseSalinity( + config, oceanClimatolgyTask, refConfig)) analyses.append(ocean.ClimatologyMapAntarcticMelt(config, - oceanClimatolgyTask)) + oceanClimatolgyTask, + refConfig)) + + analyses.append(ocean.TimeSeriesAntarcticMelt(config, oceanTimeSeriesTask, + refConfig)) - analyses.append(oceanTimeSeriesTask) - analyses.append(ocean.TimeSeriesAntarcticMelt(config, oceanTimeSeriesTask)) analyses.append(ocean.TimeSeriesTemperatureAnomaly(config, oceanTimeSeriesTask)) analyses.append(ocean.TimeSeriesSalinityAnomaly(config, oceanTimeSeriesTask)) analyses.append(ocean.TimeSeriesOHCAnomaly(config, - oceanTimeSeriesTask)) - - analyses.append(ocean.TimeSeriesSST(config, oceanTimeSeriesTask)) - analyses.append(ocean.MeridionalHeatTransport(config, oceanClimatolgyTask)) + oceanTimeSeriesTask, + refConfig)) + analyses.append(ocean.TimeSeriesSST(config, oceanTimeSeriesTask, + refConfig)) + analyses.append(ocean.MeridionalHeatTransport(config, oceanClimatolgyTask, + refConfig)) analyses.append(ocean.StreamfunctionMOC(config, oceanClimatolgyTask)) - analyses.append(ocean.IndexNino34(config, oceanIndexTask)) + analyses.append(ocean.IndexNino34(config, oceanIndexTask, refConfig)) # Sea Ice Analyses seaIceClimatolgyTask = MpasClimatologyTask(config=config, componentName='seaIce') seaIceTimeSeriesTask = MpasTimeSeriesTask(config=config, componentName='seaIce') + analyses.append(seaIceClimatolgyTask) - analyses.append(sea_ice.ClimatologyMapSeaIceConc(config, - seaIceClimatolgyTask, - hemisphere='NH')) - analyses.append(sea_ice.ClimatologyMapSeaIceThick(config, - seaIceClimatolgyTask, - hemisphere='NH')) - analyses.append(sea_ice.ClimatologyMapSeaIceConc(config, - seaIceClimatolgyTask, - hemisphere='SH')) - analyses.append(sea_ice.ClimatologyMapSeaIceThick(config, - seaIceClimatolgyTask, - hemisphere='SH')) + analyses.append(sea_ice.ClimatologyMapSeaIceConc( + config=config, mpasClimatologyTask=seaIceClimatolgyTask, + hemisphere='NH', refConfig=refConfig)) + analyses.append(sea_ice.ClimatologyMapSeaIceThick( + config=config, mpasClimatologyTask=seaIceClimatolgyTask, + hemisphere='NH', refConfig=refConfig)) + analyses.append(sea_ice.ClimatologyMapSeaIceConc( + config=config, mpasClimatologyTask=seaIceClimatolgyTask, + hemisphere='SH', refConfig=refConfig)) + analyses.append(sea_ice.ClimatologyMapSeaIceThick( + config=config, mpasClimatologyTask=seaIceClimatolgyTask, + hemisphere='SH', refConfig=refConfig)) analyses.append(seaIceTimeSeriesTask) - analyses.append(sea_ice.TimeSeriesSeaIce(config, seaIceTimeSeriesTask)) + + analyses.append(sea_ice.TimeSeriesSeaIce(config, seaIceTimeSeriesTask, + refConfig)) return analyses # }}} @@ -206,7 +221,6 @@ def add_task_and_subtasks(analysisTask, analysesToGenerate, add_task_and_subtasks(prereq, analysesToGenerate, callCheckGenerate=False) if prereq._setupStatus != 'success': - print(prereq._setupStatus) # a prereq failed setup_and_check print("ERROR: prerequisite task {} of analysis task {}" " failed during check,\n" @@ -291,6 +305,18 @@ def run_analysis(config, analyses): # {{{ Xylar Asay-Davis """ + # write the config file the log directory + logsDirectory = build_config_full_path(config, 'output', + 'logsSubdirectory') + + mainRunName = config.get('runs', 'mainRunName') + + configFileName = '{}/config.{}'.format(logsDirectory, mainRunName) + + configFile = open(configFileName, 'w') + config.write(configFile) + configFile.close() + taskCount = config.getWithDefault('execute', 'parallelTaskCount', default=1) @@ -444,19 +470,44 @@ if __name__ == "__main__": print('WARNING: Did not find config.default. Assuming other config ' 'file(s) contain a\n' 'full set of configuration options.') + defaultConfig = None configFiles = args.configFiles config = MpasAnalysisConfigParser() config.read(configFiles) if args.list: - analyses = build_analysis_list(config) + analyses = build_analysis_list(config, refConfig=None) for analysisTask in analyses: print('task: {}'.format(analysisTask.taskName)) print(' component: {}'.format(analysisTask.componentName)), print(' tags: {}'.format(', '.join(analysisTask.tags))) sys.exit(0) + if config.has_option('runs', 'referenceRunConfigFile'): + refConfigFile = config.get('runs', 'referenceRunConfigFile') + if not os.path.exists(refConfigFile): + raise OSError('A reference config file {} was specified but the ' + 'file does not exist'.format(refConfigFile)) + refConfigFiles = [refConfigFile] + if defaultConfig is not None: + refConfigFiles = [defaultConfig] + refConfigFiles + refConfig = MpasAnalysisConfigParser() + refConfig.read(refConfigFiles) + + # replace the log directory so log files get written to this run's + # log directory, not the reference run's + logsDirectory = build_config_full_path(config, 'output', + 'logsSubdirectory') + + refConfig.set('output', 'logsSubdirectory', logsDirectory) + + print('Comparing to reference run {} rather than observations. \n' + 'Make sure that MPAS-Analysis has been run previously with the ' + 'ref config file.'.format(refConfig.get('runs', 'mainRunName'))) + else: + refConfig = None + if args.purge: outputDirectory = config.get('output', 'baseDirectory') if not os.path.exists(outputDirectory): @@ -469,18 +520,22 @@ if __name__ == "__main__": if args.generate: update_generate(config, args.generate) + if refConfig is not None: + # we want to use the "generate" option from the current run, not + # the reference config file + refConfig.set('output', 'generate', config.get('output', 'generate')) + logsDirectory = build_config_full_path(config, 'output', 'logsSubdirectory') make_directories(logsDirectory) - make_directories('{}/configs/'.format(logsDirectory)) - analyses = build_analysis_list(config) + analyses = build_analysis_list(config, refConfig) analyses = determine_analyses_to_generate(analyses) if not args.setup_only and not args.html_only: run_analysis(config, analyses) if not args.setup_only: - generate_html(config, analyses) + generate_html(config, analyses, refConfig) # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python