diff --git a/config.default b/config.default index ca363e284..c8650b848 100644 --- a/config.default +++ b/config.default @@ -61,11 +61,12 @@ logsSubdirectory = logs mpasClimatologySubdirectory = clim/mpas mpasRegriddedClimSubdirectory = clim/mpas/regridded mappingSubdirectory = mapping +timeSeriesSubdirectory = timeseries # a list of analyses to generate. Valid names are: # 'timeSeriesOHC', 'timeSeriesSST', 'regriddedSST', -# 'regriddedSSS', 'regriddedMLD', 'timeSeriesSeaIceAreaVol', -# 'regriddedSeaIceConcThick' +# 'regriddedSSS', 'regriddedMLD', 'streamfunctionMOC', +# 'timeSeriesSeaIceAreaVol', 'regriddedSeaIceConcThick' # the following shortcuts exist: # 'all' -- all analyses will be run # 'all_timeSeries' -- all time-series analyses will be run @@ -85,7 +86,8 @@ generate = ['all'] # alternative examples that would perform all analysis except # 'timeSeriesOHC' -#generate = ['timeSeriesSST', 'all_regriddedHorizontal', 'all_seaIce'] +#generate = ['timeSeriesSST', 'streamfunctionMOC', +# 'all_regriddedHorizontal', 'all_seaIce'] #generate = ['all', 'no_timeSeriesOHC'] # Each subsequent list entry can be used to alter previous list entries. For # example, the following would produce all analyses except regriddedSST, @@ -304,19 +306,47 @@ regionIndicesToPlot = [6] # window) movingAveragePoints = 12 -[timeSeriesMOC] -## options related to plotting time series of meridional overturning +[streamfunctionMOC] +## options related to plotting the streamfunction of the meridional overturning ## circulation (MOC) +# Region names for basin MOC calculation. +# Supported options are Atlantic and IndoPacific +regionNames = ['Atlantic'] + +# Mask file for post-processing regional MOC computation +regionMaskFiles = /path/to/MOCregional/mapping/file + +# Size of latitude bins over which MOC streamfunction is integrated +latBinSizeGlobal = 1. +latBinSizeAtlantic = 0.5 +latBinSizeIndoPacific = 0.5 + +# colormap for model results +colormapNameGlobal = RdYlBu_r +colormapNameAtlantic = RdYlBu_r +colormapNameIndoPacific = RdYlBu_r +# colormap indices for contour color +colormapIndicesGlobal = [0, 40, 80, 110, 140, 170, 200, 230, 255] +colormapIndicesAtlantic = [0, 40, 80, 110, 140, 170, 200, 230, 255] +colormapIndicesIndoPacific = [0, 40, 80, 110, 140, 170, 200, 230, 255] +# colorbar levels/values for contour boundaries +colorbarLevelsGlobal = [-20, -10, -5, -2, 2, 5, 10, 20, 30, 40] +colorbarLevelsAtlantic = [-10, -5, -2, 0, 5, 8, 10, 14, 18, 22] +colorbarLevelsIndoPacific = [-10, -5, -2, 0, 5, 8, 10, 14, 18, 22] +# contour line levels +contourLevelsGlobal = np.arange(-25.1,35.1,10) +contourLevelsAtlantic = np.arange(-8,20.1,2) +contourLevelsIndoPacific = np.arange(-8,20.1,2) + ## compare to output from another model run? #compareWithModel = True # compare to observations? compareWithObservations = True -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] -# Number of points over which to compute moving average (e.g., for monthly -# output, movingAveragePoints=12 corresponds to a 12-month moving average -# window) + +# Number of points over which to compute moving average for +# MOC timeseries (e.g., for monthly output, movingAveragePoints=12 +# corresponds to a 12-month moving average window) movingAveragePoints = 12 [timeSeriesSeaIceAreaVol] @@ -342,18 +372,18 @@ polarPlot = False ## (SST) against reference model results and observations # colormap for model/observations -resultColormap = RdYlBu_r -# indices into resultColormap for contour color -resultColormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255] +colormapNameResult = RdYlBu_r +# color indices into colormapName for filled contours +colormapIndicesResult = [0, 40, 80, 110, 140, 170, 200, 230, 255] # colormap levels/values for contour boundaries -resultContourValues = [-2, 0, 2, 6, 10, 16, 22, 26, 28, 32] +colorbarLevelsResult = [-2, 0, 2, 6, 10, 16, 22, 26, 28, 32] # colormap for differences -differenceColormap = RdBu_r -# indices into differenceColormap for contour color -differenceColormapIndices = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] +colormapNameDifference = RdBu_r +# color indices into colormapName for filled contours +colormapIndicesDifference = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] # colormap levels/values for contour boundaries -differenceContourValues = [-5, -3, -2, -1, 0, 1, 2, 3, 5] +colorbarLevelsDifference = [-5, -3, -2, -1, 0, 1, 2, 3, 5] # Times for comparison times (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, JFM, AMJ, JAS, OND, ANN) comparisonTimes = ['JFM', 'JAS', 'ANN'] @@ -363,18 +393,18 @@ comparisonTimes = ['JFM', 'JAS', 'ANN'] ## (SSS) against reference model results and observations # colormap for model/observations -resultColormap = RdYlBu_r -# indices into resultColormap for contour color -resultColormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255] +colormapNameResult = RdYlBu_r +# color indices into colormapName for filled contours +colormapIndicesResult = [0, 40, 80, 110, 140, 170, 200, 230, 255] # colormap levels/values for contour boundaries -resultContourValues = [28, 29, 30, 31, 32, 33, 34, 35, 36, 38] +colorbarLevelsResult = [28, 29, 30, 31, 32, 33, 34, 35, 36, 38] # colormap for differences -differenceColormap = RdBu_r -# indices into differenceColormap for contour color -differenceColormapIndices = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] +colormapNameDifference = RdBu_r +# color indices into colormapName for filled contours +colormapIndicesDifference = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] # colormap levels/values for contour boundaries -differenceContourValues = [-3, -2, -1, -0.5, 0, 0.5, 1, 2, 3] +colorbarLevelsDifference = [-3, -2, -1, -0.5, 0, 0.5, 1, 2, 3] # Times for comparison times (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, JFM, AMJ, JAS, OND, ANN) comparisonTimes = ['JFM', 'JAS', 'ANN'] @@ -384,18 +414,18 @@ comparisonTimes = ['JFM', 'JAS', 'ANN'] ## (MLD) against reference model results and observations # colormap for model/observations -resultColormap = viridis -# indices into resultColormap for contour color -resultColormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255] +colormapNameResult = viridis +# color indices into colormapName for filled contours +colormapIndicesResult = [0, 40, 80, 110, 140, 170, 200, 230, 255] # colormap levels/values for contour boundaries -resultContourValues = [0, 20, 40, 60, 80, 100, 150, 200, 400, 800] +colorbarLevelsResult = [0, 20, 40, 60, 80, 100, 150, 200, 400, 800] # colormap for differences -differenceColormap = RdBu_r -# indices into differenceColormap for contour color -differenceColormapIndices = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] +colormapNameDifference = RdBu_r +# color indices into colormapName for filled contours +colormapIndicesDifference = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] # colormap levels/values for contour boundaries -differenceContourValues = [-150, -80, -30, -10, 0, 10, 30, 80, 150] +colorbarLevelsDifference = [-150, -80, -30, -10, 0, 10, 30, 80, 150] # Times for comparison times (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, JFM, AMJ, JAS, OND, ANN) comparisonTimes = ['JFM', 'JAS', 'ANN'] @@ -405,28 +435,40 @@ comparisonTimes = ['JFM', 'JAS', 'ANN'] ## and thickness against reference model results and observations # colormap for model/observations -resultColormap = inferno -# indices into resultColormap for contour color -resultColormapIndices = [20, 80, 110, 140, 170, 200, 230, 255] +colormapNameConcResultWinter = inferno +colormapNameConcResultSummer = inferno +colormapNameThickResultNH = inferno +colormapNameThickResultSH = inferno +# color indices into colormapName for filled contours +colormapIndicesConcResultWinter = [20, 80, 110, 140, 170, 200, 230, 255] +colormapIndicesConcResultSummer = [20, 80, 110, 140, 170, 200, 230, 255] +colormapIndicesThickResultNH = [20, 80, 110, 140, 170, 200, 230, 255] +colormapIndicesThickResultSH = [20, 80, 110, 140, 170, 200, 230, 255] # colormap levels/values for contour boundaries for: # concentration in winter and summer -resultConcWinterContourValues = [0.15, 0.4, 0.7, 0.9, 0.94, 0.96, 0.98, 0.99, 1] -resultConcSummerContourValues = [0.15, 0.3, 0.5, 0.7, 0.8, 0.85, 0.9, 0.95, 1] +colorbarLevelsConcResultWinter = [0.15, 0.4, 0.7, 0.9, 0.94, 0.96, 0.98, 0.99, 1] +colorbarLevelsConcResultSummer = [0.15, 0.3, 0.5, 0.7, 0.8, 0.85, 0.9, 0.95, 1] # thickness in the northern and southern hemispheres -resultThickNHContourValues = [0, 0.25, 0.5, 1, 1.5, 2, 2.5, 3, 3.5] -resultThickSHContourValues = [0, 0.2, 0.4, 0.6, 0.8, 1, 1.5, 2, 2.5] +colorbarLevelsThickResultNH = [0, 0.25, 0.5, 1, 1.5, 2, 2.5, 3, 3.5] +colorbarLevelsThickResultSH = [0, 0.2, 0.4, 0.6, 0.8, 1, 1.5, 2, 2.5] # colormap for differences -differenceColormap = RdBu_r -# indices into differenceColormap for contour color -differenceColormapIndices = [0, 40, 80, 127, 127, 170, 210, 255] +colormapNameConcDifferenceWinter = RdBu_r +colormapNameConcDifferenceSummer = RdBu_r +colormapNameThickDifferenceNH = RdBu_r +colormapNameThickDifferenceSH = RdBu_r +# color indices into colormapName for filled contours +colormapIndicesConcDifferenceWinter = [0, 40, 80, 127, 127, 170, 210, 255] +colormapIndicesConcDifferenceSummer = [0, 40, 80, 127, 127, 170, 210, 255] +colormapIndicesThickDifferenceNH = [0, 40, 80, 127, 127, 170, 210, 255] +colormapIndicesThickDifferenceSH = [0, 40, 80, 127, 127, 170, 210, 255] # colormap levels/values for contour boundaries for: # concentration in winter and summer -differenceConcWinterContourValues = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8] -differenceConcSummerContourValues = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8] +colorbarLevelsConcDifferenceWinter = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8] +colorbarLevelsConcDifferenceSummer = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8] # thickness in the northern and southern hemispheres -differenceThickNHContourValues = [-2.5, -2, -0.5, -0.1, 0, 0.1, 0.5, 2, 2.5] -differenceThickSHContourValues = [-2.5, -2, -0.5, -0.1, 0, 0.1, 0.5, 2, 2.5] +colorbarLevelsThickDifferenceNH = [-2.5, -2, -0.5, -0.1, 0, 0.1, 0.5, 2, 2.5] +colorbarLevelsThickDifferenceSH = [-2.5, -2, -0.5, -0.1, 0, 0.1, 0.5, 2, 2.5] # reference lat/lon for sea ice plots in the northern hemisphere minimumLatitudeNH = 50 diff --git a/configs/edison/config.20161006bugfix.alpha8.A_WCYCL1850S.ne30_oEC_ICG.edison b/configs/edison/config.20161006bugfix.alpha8.A_WCYCL1850S.ne30_oEC_ICG.edison index 4804d0e22..a225448da 100644 --- a/configs/edison/config.20161006bugfix.alpha8.A_WCYCL1850S.ne30_oEC_ICG.edison +++ b/configs/edison/config.20161006bugfix.alpha8.A_WCYCL1850S.ne30_oEC_ICG.edison @@ -4,10 +4,6 @@ # mainRunName is a name that identifies the simulation being analyzed. mainRunName = 20161006bugfix.alpha8.A_WCYCL1850S.ne30_oEC_ICG.edison -# 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 @@ -20,11 +16,6 @@ preprocessedReferenceRunName = B1850C5_ne30_v0.4 # directory containing model results baseDirectory = /global/cscratch1/sd/jonbob/ACME_simulations/20161006bugfix.alpha8.A_WCYCL1850S.ne30_oEC_ICG.edison/run -# names of namelist and streams files. If not in baseDirectory, give full path -oceanNamelistFileName = mpas-o_in -oceanStreamsFileName = streams.ocean -seaIceNamelistFileName = mpas-cice_in -seaIceStreamsFileName = streams.cice # names of ocean and sea ice meshes (e.g. EC60to30, QU240, RRS30to10, etc.) mpasMeshName = EC60to30 @@ -121,3 +112,15 @@ baseDirectory = /global/project/projectdirs/acme/observations/SeaIce # directory where ocean reference simulation results are stored baseDirectory = /global/project/projectdirs/acme/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing +[timeSeriesSeaIceAreaVol] +## options related to plotting time series of sea ice area and volume + +# plot on polar plot +polarPlot = False + +[streamfunctionMOC] +## options related to plotting the streamfunction of the meridional overturning +## circulation (MOC) + +# Mask file for ocean basin regional computation +regionMaskFiles = /global/project/projectdirs/acme/mapping/grids/EC60to30v1_SingleRegionAtlanticWTransportTransects_masks.nc diff --git a/configs/edison/config.20161117.beta0.A_WCYCL1850.ne30_oEC.edison b/configs/edison/config.20161117.beta0.A_WCYCL1850.ne30_oEC.edison index 3185763a7..466330840 100644 --- a/configs/edison/config.20161117.beta0.A_WCYCL1850.ne30_oEC.edison +++ b/configs/edison/config.20161117.beta0.A_WCYCL1850.ne30_oEC.edison @@ -4,10 +4,6 @@ # mainRunName is a name that identifies the simulation being analyzed. mainRunName = 20161117.beta0.A_WCYCL1850.ne30_oEC.edison -# 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 @@ -20,11 +16,6 @@ preprocessedReferenceRunName = B1850C5_ne30_v0.4 # directory containing model results baseDirectory = /scratch2/scratchdirs/golaz/ACME_simulations/20161117.beta0.A_WCYCL1850.ne30_oEC.edison/run -# names of namelist and streams files. If not in baseDirectory, give full path -oceanNamelistFileName = mpas-o_in -oceanStreamsFileName = streams.ocean -seaIceNamelistFileName = mpas-cice_in -seaIceStreamsFileName = streams.cice # names of ocean and sea ice meshes (e.g. EC60to30, QU240, RRS30to10, etc.) mpasMeshName = EC60to30 @@ -122,3 +113,15 @@ baseDirectory = /global/project/projectdirs/acme/observations/SeaIce # directory where ocean reference simulation results are stored baseDirectory = /global/project/projectdirs/acme/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing +[timeSeriesSeaIceAreaVol] +## options related to plotting time series of sea ice area and volume + +# plot on polar plot +polarPlot = False + +[streamfunctionMOC] +## options related to plotting the streamfunction of the meridional overturning +## circulation (MOC) + +# Mask file for ocean basin regional computation +regionMaskFiles = /global/project/projectdirs/acme/mapping/grids/EC60to30v1_SingleRegionAtlanticWTransportTransects_masks.nc diff --git a/configs/edison/config.20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison b/configs/edison/config.20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison new file mode 100644 index 000000000..e63e6b4a4 --- /dev/null +++ b/configs/edison/config.20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison @@ -0,0 +1,127 @@ +[runs] +## options related to the run to be analyzed and reference runs to be +## compared against + +# mainRunName is a name that identifies the simulation being analyzed. +mainRunName = 20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison +# 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 +# performed with MPAS components (so they cannot be easily ingested by +# MPAS-Analysis) +preprocessedReferenceRunName = B1850C5_ne30_v0.4 + +[input] +## options related to reading in the results to be analyzed + +# directory containing model results +baseDirectory = /scratch2/scratchdirs/golaz/ACME_simulations/20170313.beta1.A_WCYCL1850S.ne30_oECv3_ICG.edison/run + +# names of ocean and sea ice meshes (e.g. EC60to30, QU240, RRS30to10, etc.) +mpasMeshName = EC60to30v3 + +[output] +## options related to writing out plots, intermediate cached data sets, logs, +## etc. + +# directory where analysis should be written +baseDirectory = /dir/to/analysis/output + +# a list of analyses to generate. Valid names are: +# 'timeSeriesOHC', 'timeSeriesSST', 'regriddedSST', +# 'regriddedSSS', 'regriddedMLD', 'timeSeriesSeaIceAreaVol', +# 'regriddedSeaIceConcThick' +# the following shortcuts exist: +# 'all' -- all analyses will be run +# 'all_timeSeries' -- all time-series analyses will be run +# 'all_regriddedHorizontal' -- all analyses involving regridded horizontal +# fields will be run +# 'all_ocean' -- all ocean analyses will be run +# 'all_seaIce' -- all sea-ice analyses will be run +# 'no_timeSeriesOHC' -- skip 'timeSeriesOHC' (and similarly with the +# other analyses). +# 'no_ocean', 'no_timeSeries', etc. -- in analogy to 'all_*', skip the +# given category of analysis +# an equivalent syntax can be used on the command line to override this +# option: +# ./run_analysis.py config.analysis --generate \ +# all,no_ocean,all_timeSeries +generate = ['all'] + +# alternative examples that would perform all analysis except +# 'timeSeriesOHC' +#generate = ['timeSeriesSST', 'all_regriddedHorizontal', 'all_seaIce'] +#generate = ['all', 'no_timeSeriesOHC'] +# Each subsequent list entry can be used to alter previous list entries. For +# example, the following would produce all analyses except regriddedSST, +# regriddedSSS and regriddedMLD (albeit not in a very intuitive way): +#generate = ['all', 'no_ocean', 'all_timeSeries'] + +[climatology] +## options related to producing climatologies, typically to compare against +## observations and previous runs + +# the first year over which to average climatalogies +startYear = 11 +# the last year over which to average climatalogies +endYear = 20 + +# The names of the mapping file used for interpolation. If a mapping file has +# already been generated, supplying the absolute path can save the time of +# generating a new one. If nothing is supplied, the file name is automatically +# generated based on the MPAS mesh name, the comparison grid resolution, and +# the interpolation method +mpasMappingFile = /global/project/projectdirs/acme/mapping/maps/map_oEC60to30v3_TO_0.5x0.5degree_blin.nc + +[timeSeries] +## options related to producing time series plots, often to compare against +## observations and previous runs + +# start and end years for timeseries analysis. Using out-of-bounds values +# like start_year = 1 and end_year = 9999 will be clipped to the valid range +# of years, and is a good way of insuring that all values are used. +startYear = 1 +endYear = 22 + +[oceanObservations] +## options related to ocean observations with which the results will be compared + +# directory where ocean observations are stored +baseDirectory = /global/project/projectdirs/acme/observations/Ocean/ +sstSubdirectory = SST +sssSubdirectory = SSS +mldSubdirectory = MLD + +[oceanPreprocessedReference] +## options related to preprocessed ocean reference run with which the results +## will be compared (e.g. a POP, CESM or ACME v0 run) + +# directory where ocean reference simulation results are stored +baseDirectory = /global/project/projectdirs/acme/ACMEv0_lowres/B1850C5_ne30_v0.4/ocn/postprocessing + +[seaIceObservations] +## options related to sea ice observations with which the results will be +## compared + +# directory where sea ice observations are stored +baseDirectory = /global/project/projectdirs/acme/observations/SeaIce + +[seaIcePreprocessedReference] +## options related to preprocessed sea ice reference run with which the results +## will be compared (e.g. a CICE, CESM or ACME v0 run) + +# directory where ocean reference simulation results are stored +baseDirectory = /global/project/projectdirs/acme/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing + +[timeSeriesSeaIceAreaVol] +## options related to plotting time series of sea ice area and volume + +# plot on polar plot +polarPlot = False + +[streamfunctionMOC] +## options related to plotting the streamfunction of the meridional overturning +## circulation (MOC) + +# Mask file for ocean basin regional computation +regionMaskFiles = /global/project/projectdirs/acme/mapping/grids/EC60to30v3_SingleRegionAtlanticWTransportTransects_masks.nc diff --git a/configs/lanl/config.20161117.beta0.A_WCYCL1850S.ne30_oEC_ICG.edison b/configs/lanl/config.20161117.beta0.A_WCYCL1850S.ne30_oEC_ICG.edison index 4ad368482..035ccd206 100644 --- a/configs/lanl/config.20161117.beta0.A_WCYCL1850S.ne30_oEC_ICG.edison +++ b/configs/lanl/config.20161117.beta0.A_WCYCL1850S.ne30_oEC_ICG.edison @@ -4,10 +4,6 @@ # mainRunName is a name that identifies the simulation being analyzed. mainRunName = 20161117.beta0.A_WCYCL1850S.ne30_oEC_ICG.edison -# 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 @@ -107,3 +103,15 @@ baseDirectory = /usr/projects/climate/SHARED_CLIMATE/observations/SeaIce # directory where ocean reference simulation results are stored baseDirectory = /usr/projects/climate/SHARED_CLIMATE/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing +[timeSeriesSeaIceAreaVol] +## options related to plotting time series of sea ice area and volume + +# plot on polar plot +polarPlot = False + +[streamfunctionMOC] +## options related to plotting the streamfunction of the meridional overturning +## circulation (MOC) + +# Mask file for ocean basin regional computation +regionMaskFiles = /usr/projects/climate/mpeterse/analysis_input_files/EC60to30v1/EC60to30v1_SingleRegionAtlanticWTransportTransects_masks.nc diff --git a/configs/lanl/config.20170120.beta0.GMPAS-QU240.wolf b/configs/lanl/config.20170120.beta0.GMPAS-QU240.wolf index f63a107fe..98d17b047 100644 --- a/configs/lanl/config.20170120.beta0.GMPAS-QU240.wolf +++ b/configs/lanl/config.20170120.beta0.GMPAS-QU240.wolf @@ -4,10 +4,6 @@ # mainRunName is a name that identifies the simulation being analyzed. mainRunName = 20170120.beta0.GMPAS-QU240.wolf -# 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 @@ -107,3 +103,8 @@ baseDirectory = /usr/projects/climate/SHARED_CLIMATE/observations/SeaIce # directory where ocean reference simulation results are stored baseDirectory = /usr/projects/climate/SHARED_CLIMATE/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing +[timeSeriesSeaIceAreaVol] +## options related to plotting time series of sea ice area and volume + +# plot on polar plot +polarPlot = False diff --git a/configs/lanl/config.20170207.MPAS-SeaIce.QU60km_polar.wolf b/configs/lanl/config.20170207.MPAS-SeaIce.QU60km_polar.wolf index 850dd7f11..9af6b5f2b 100644 --- a/configs/lanl/config.20170207.MPAS-SeaIce.QU60km_polar.wolf +++ b/configs/lanl/config.20170207.MPAS-SeaIce.QU60km_polar.wolf @@ -4,10 +4,6 @@ # mainRunName is a name that identifies the simulation being analyzed. mainRunName = MPAS-SeaIce.QU60km_polar -# 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 @@ -20,11 +16,6 @@ preprocessedReferenceRunName = B1850C5_ne30_v0.4 # directory containing model results baseDirectory = /net/scratch2/akt/MPAS/rundirs/rundir_QU60km_polar -# names of namelist and streams files. If not in baseDirectory, give full path -oceanNamelistFileName = mpas-o_in -oceanStreamsFileName = streams.ocean -seaIceNamelistFileName = namelist.cice -seaIceStreamsFileName = streams.cice # names of ocean and sea ice meshes (e.g. EC60to30, QU240, RRS30to10, etc.) mpasMeshName = QU60 @@ -112,3 +103,8 @@ baseDirectory = /usr/projects/climate/SHARED_CLIMATE/observations/SeaIce # directory where ocean reference simulation results are stored baseDirectory = /usr/projects/climate/SHARED_CLIMATE/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing +[timeSeriesSeaIceAreaVol] +## options related to plotting time series of sea ice area and volume + +# plot on polar plot +polarPlot = False diff --git a/configs/olcf/config.20161117.beta0.A_WCYCL1850S.ne30_oEC_ICG.titan b/configs/olcf/config.20161117.beta0.A_WCYCL1850S.ne30_oEC_ICG.titan index 38b080d22..4ac3b63ba 100644 --- a/configs/olcf/config.20161117.beta0.A_WCYCL1850S.ne30_oEC_ICG.titan +++ b/configs/olcf/config.20161117.beta0.A_WCYCL1850S.ne30_oEC_ICG.titan @@ -64,11 +64,7 @@ generate = ['all'] # the first year over which to average climatalogies startYear = 131 # the last year over which to average climatalogies -endYear = 20 - -# The comparison grid resolution in degrees -comparisonLatResolution = 0.5 -comparisonLonResolution = 0.5 +endYear = 140 # The names of the mapping file used for interpolation. If a mapping file has # already been generated, supplying the absolute path can save the time of @@ -77,14 +73,6 @@ comparisonLonResolution = 0.5 # the interpolation method # mpasMappingFile = /path/to/mapping/file -# overwrite files when building climatologies? -overwriteMapping = False -overwriteMpasClimatology = False - -# interpolation order for model and observation results. Likely values are -# 'bilinear', 'neareststod' (nearest neighbor) or 'conserve' -mpasInterpolationMethod = bilinear - [timeSeries] ## options related to producing time series plots, often to compare against ## observations and previous runs @@ -100,6 +88,9 @@ endYear = 140 # directory where ocean observations are stored baseDirectory = /lustre/atlas/proj-shared/cli115/observations +sstSubdirectory = SST +sssSubdirectory = SSS +mldSubdirectory = MLD [oceanPreprocessedReference] ## options related to preprocessed ocean reference run with which the results @@ -114,52 +105,6 @@ baseDirectory = /lustre/atlas/proj-shared/cli115/milena/ACMEv0_lowres/B1850C5_ne # directory where sea ice observations are stored baseDirectory = /lustre/atlas/proj-shared/cli115/observations/SeaIce -areaNH = IceArea_timeseries/iceAreaNH_climo.nc -areaSH = IceArea_timeseries/iceAreaSH_climo.nc -volNH = PIOMAS/PIOMASvolume_monthly_climo.nc -volSH = none -concentrationNASATeamNH_JFM = SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_NH_jfm.interp0.5x0.5.nc -concentrationNASATeamNH_JAS = SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_NH_jas.interp0.5x0.5.nc -concentrationNASATeamSH_DJF = SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_SH_djf.interp0.5x0.5.nc -concentrationNASATeamSH_JJA = SSMI/NASATeam_NSIDC0051/SSMI_NASATeam_gridded_concentration_SH_jja.interp0.5x0.5.nc -concentrationBootstrapNH_JFM = SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_NH_jfm.interp0.5x0.5.nc -concentrationBootstrapNH_JAS = SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_NH_jas.interp0.5x0.5.nc -concentrationBootstrapSH_DJF = SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_SH_djf.interp0.5x0.5.nc -concentrationBootstrapSH_JJA = SSMI/Bootstrap_NSIDC0079/SSMI_Bootstrap_gridded_concentration_SH_jja.interp0.5x0.5.nc -thicknessNH_ON = ICESat/ICESat_gridded_mean_thickness_NH_on.interp0.5x0.5.nc -thicknessNH_FM = ICESat/ICESat_gridded_mean_thickness_NH_fm.interp0.5x0.5.nc -thicknessSH_ON = ICESat/ICESat_gridded_mean_thickness_SH_on.interp0.5x0.5.nc -thicknessSH_FM = ICESat/ICESat_gridded_mean_thickness_SH_fm.interp0.5x0.5.nc - -# The name of mapping files used for interpolating observations to the -# comparison grid. Interpolation is only performed if the observation grid has -# a different resolution from the comparison grid. If nothing is supplied, the -# file name is automatically generated based on the MPAS mesh name, the -# comparison grid resolution, and the interpolation method -# seaIceClimatologyMappingFile = /path/to/mapping/file - -# interpolation order for observations. Likely values are -# 'bilinear', 'neareststod' (nearest neighbor) or 'conserve' -interpolationMethod = bilinear - -# The directories where observation climatologies will be stored if they need -# to be computed. If a relative path is supplied, it is relative to the output -# base directory. If an absolute path is supplied, this should point to -# cached climatology files on the desired comparison grid, in which case -# overwriteObsClimatology should be False. If cached regridded files are -# supplied, there is no need to provide cached files before regridding. -climatologySubdirectory = clim/obs -regriddedClimSubdirectory = clim/obs/regridded - -# overwrite files when building climatologies? -overwriteObsClimatology = False - -[seaIceReference] -## options related to sea ice reference run with which the results will be -## compared - -# directory where sea ice reference simulation results are stored -baseDirectory = /lustre/atlas/proj-shared/cli115/milena/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing [seaIcePreprocessedReference] ## options related to preprocessed sea ice reference run with which the results @@ -168,221 +113,15 @@ baseDirectory = /lustre/atlas/proj-shared/cli115/milena/ACMEv0_lowres/B1850C5_ne # directory where ocean reference simulation results are stored baseDirectory = /lustre/atlas/proj-shared/cli115/milena/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing -[timeSeriesOHC] -## options related to plotting time series of ocean heat content (OHC) - -## compare to output from another model run? -#compareWithModel = True -# compare to observations? -compareWithObservations = False -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] -# Number of points over which to compute moving average (e.g., for monthly -# output, movingAveragePoints=12 corresponds to a 12-month moving average -# window) -movingAveragePoints = 12 - -[timeSeriesSST] -## options related to plotting time series of sea surface temperature (SST) - -## compare to output from another model run? -#compareWithModel = True -# compare to observations? -compareWithObservations = True -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] -# Number of points over which to compute moving average (e.g., for monthly -# output, movingAveragePoints=12 corresponds to a 12-month moving average -# window) -movingAveragePoints = 12 - -[timeSeriesNino34] -## options related to plotting time series of the El Nino 3.4 index - -## compare to output from another model run? -#compareWithModel = True -# compare to observations? -compareWithObservations = True -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] -# Number of points over which to compute moving average (e.g., for monthly -# output, movingAveragePoints=12 corresponds to a 12-month moving average -# window) -movingAveragePoints = 12 - -[timeSeriesMHT] -## options related to plotting time series of meridional heat transport (MHT) - -## compare to output from another model run? -#compareWithModel = True -# compare to observations? -compareWithObservations = True -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] -# Number of points over which to compute moving average (e.g., for monthly -# output, movingAveragePoints=12 corresponds to a 12-month moving average -# window) -movingAveragePoints = 12 - -[timeSeriesMOC] -## options related to plotting time series of meridional overturning -## circulation (MOC) - -## compare to output from another model run? -#compareWithModel = True -# compare to observations? -compareWithObservations = True -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] -# Number of points over which to compute moving average (e.g., for monthly -# output, movingAveragePoints=12 corresponds to a 12-month moving average -# window) -movingAveragePoints = 12 - [timeSeriesSeaIceAreaVol] ## options related to plotting time series of sea ice area and volume -## compare to output from another model run? -#compareWithModel = True -# compare to observations? -compareWithObservations = True -# list of region indices to plot from the region list in [regions] below -regionIndicesToPlot = [6] -# Number of points over which to compute moving average (e.g., for monthly -# output, movingAveragePoints=12 corresponds to a 12-month moving average -# window) -movingAveragePoints = 1 -# title font properties -titleFontSize = 18 +# plot on polar plot +polarPlot = False -[regriddedSST] -## options related to plotting horizontally regridded sea surface temperature -## (SST) against reference model results and observations - -# colormap for model/observations -resultColormap = RdYlBu_r -# indices into resultColormap for contour color -resultColormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255] -# colormap levels/values for contour boundaries -resultContourValues = [-2, 0, 2, 6, 10, 16, 22, 26, 28, 32] - -# colormap for differences -differenceColormap = RdBu_r -# indices into differenceColormap for contour color -differenceColormapIndices = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] -# colormap levels/values for contour boundaries -differenceContourValues = [-5, -3, -2, -1, 0, 1, 2, 3, 5] - -# alternative colormap/index suggestions -#resultColormap = viridis -#differenceColormap = RdBu_r -#differenceColormapIndices = [0, 40, 80, 127, 170, 210, 255] -#differenceContourValues = [-3, -2, -1, -0.5, 0.5, 1, 2, 3] - -# Times for comparison times (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, JFM, AMJ, JAS, OND, ANN) -comparisonTimes = ['JFM', 'JAS', 'ANN'] - -[regriddedSSS] -## options related to plotting horizontally regridded sea surface salinity -## (SSS) against reference model results and observations - -# colormap for model/observations -resultColormap = RdYlBu_r -# indices into resultColormap for contour color -resultColormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255] -# colormap levels/values for contour boundaries -resultContourValues = [28, 29, 30, 31, 32, 33, 34, 35, 36, 38] - -# colormap for differences -differenceColormap = RdBu_r -# indices into differenceColormap for contour color -differenceColormapIndices = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] -# colormap levels/values for contour boundaries -differenceContourValues = [-3, -2, -1, -0.5, 0, 0.5, 1, 2, 3] - -# Times for comparison times (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, JFM, AMJ, JAS, OND, ANN) -comparisonTimes = ['JFM', 'JAS', 'ANN'] - -[regriddedMLD] -## options related to plotting horizontally regridded mixed layer depth -## (MLD) against reference model results and observations - -# colormap for model/observations -resultColormap = viridis -# indices into resultColormap for contour color -resultColormapIndices = [0, 40, 80, 110, 140, 170, 200, 230, 255] -# colormap levels/values for contour boundaries -resultContourValues = [0, 20, 40, 60, 80, 100, 150, 200, 400, 800] - -###### Xylar's note: I think something is off here. I think there should -###### be one more contour value than there are indices but there are 2 in MLD -# colormap for differences -differenceColormap = RdBu_r -# indices into differenceColormap for contour color -differenceColormapIndices = [0, 28, 57, 85, 113, 142, 170, 198, 227, 255] -# colormap levels/values for contour boundaries -differenceContourValues = [-150, -80, -30, -10, 0, 10, 30, 80, 150] - -# Times for comparison times (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, JFM, AMJ, JAS, OND, ANN) -comparisonTimes = ['JFM', 'JAS', 'ANN'] - -[regriddedSeaIceConcThick] -## options related to plotting horizontally regridded sea ice concentration -## and thickness against reference model results and observations - -# colormap for model/observations -resultColormap = inferno -# indices into resultColormap for contour color -resultColormapIndices = [20, 80, 110, 140, 170, 200, 230, 255] -# colormap levels/values for contour boundaries for: -# concentration in winter and summer -resultConcWinterContourValues = [0.15, 0.4, 0.7, 0.9, 0.94, 0.96, 0.98, 0.99, 1] -resultConcSummerContourValues = [0.15, 0.3, 0.5, 0.7, 0.8, 0.85, 0.9, 0.95, 1] -# thickness in the northern and southern hemispheres -resultThickNHContourValues = [0, 0.25, 0.5, 1, 1.5, 2, 2.5, 3, 3.5] -resultThickSHContourValues = [0, 0.2, 0.4, 0.6, 0.8, 1, 1.5, 2, 2.5] - -# colormap for differences -differenceColormap = RdBu_r -# indices into differenceColormap for contour color -differenceColormapIndices = [0, 40, 80, 127, 127, 170, 210, 255] -# colormap levels/values for contour boundaries for: -# concentration in winter and summer -differenceConcWinterContourValues = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8] -differenceConcSummerContourValues = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8] -# thickness in the northern and southern hemispheres -differenceThickNHContourValues = [-2.5, -2, -0.5, -0.1, 0, 0.1, 0.5, 2, 2.5] -differenceThickSHContourValues = [-2.5, -2, -0.5, -0.1, 0, 0.1, 0.5, 2, 2.5] - -# reference lat/lon for sea ice plots in the northern hemisphere -minimumLatitudeNH = 50 -referenceLongitudeNH = 0 -# reference lat/lon for sea ice plots in the southern hemisphere -minimumLatitudeSH = -50 -referenceLongitudeSH = 180 - -[regions] -## options related to ocean regions used in several analysis modules - -# list of region names (needs to be in the same order as region indices in -# time-series stats) -regions = ['arctic', 'equatorial', 'so', 'nino3', 'nino4', 'nino3.4', 'global'] -# list of plot titles (needs to be in the same order as region indices in -# time-series stats) -plotTitles = ['Arctic', 'Equatorial (15S-15N)', 'Southern Ocean', 'Nino 3', 'Nino 4', 'Nino 3.4', 'Global Ocean'] - -[plot] -## options related to plotting that are the defaults across all analysis modules - -# set to true if you want plots to be displayed (one by one) to the screen in -# addition to being written out to png files -# Note: displayToScreen = True seems to hang on Edison on large data sets, -# so suggested use is just for debugging either locally or with small data sets -displayToScreen = False +[streamfunctionMOC] +## options related to plotting the streamfunction of the meridional overturning +## circulation (MOC) -# font size on axes -axisFontSize = 16 -# title font properties -titleFontSize = 20 -titleFontColor = black -titleFontWeight = normal +# Mask file for ocean basin regional computation +regionMaskFiles = /lustre/atlas1/cli115/proj-shared/mapping/grids/EC60to30v1_SingleRegionAtlanticWTransportTransects_masks.nc diff --git a/mpas_analysis/ocean/meridional_overturning_circulation.py b/mpas_analysis/ocean/meridional_overturning_circulation.py new file mode 100644 index 000000000..09d98330b --- /dev/null +++ b/mpas_analysis/ocean/meridional_overturning_circulation.py @@ -0,0 +1,505 @@ +""" +Computation and plotting of model meridional overturning circulation. +Will eventually support: + * MOC streamfunction, post-processed (currently supported) + * MOC streamfunction, from MOC analysis member + * MOC time series (max value at 24.5N), post-processed + * MOC time series (max value at 24.5N), from MOC analysis member + +Authors +------- +Milena Veneziani, Mark Petersen, Phillip Wolfram + +Last Modified +------------- +03/14/2017 +""" + +import xarray as xr +import numpy as np +import netCDF4 +import os + +from ..shared.constants.constants import m3ps_to_Sv, rad_to_deg +from ..shared.plot.plotting import plot_vertical_section,\ + timeseries_analysis_plot, setup_colormap + +from ..shared.io import NameList, StreamsFile +from ..shared.io.utility import buildConfigFullPath + +from ..shared.generalized_reader.generalized_reader \ + import open_multifile_dataset + +from ..shared.timekeeping.utility import get_simulation_start_time, \ + days_to_datetime + +from ..shared.climatology import climatology + + +def moc_streamfunction(config, streamMap=None, variableMap=None): # {{{ + """ + Process MOC analysis member data if available, or compute MOC at + post-processing if not. Plots streamfunction climatolgoical sections + as well as time series of max Atlantic MOC at 26.5N (latitude of + RAPID MOC Array). + + Parameters + ---------- + config : instance of MpasAnalysisConfigParser + configuration options used to customize the analysis task + + streamMap : dict, optional + a dictionary of MPAS-O stream names that map to + their mpas_analysis counterparts. + + variableMap : dict, optional + a dictionary of MPAS-O variable names that map + to their mpas_analysis counterparts. + + Authors + ------- + Milena Veneziani, Mark Petersen, Phillip J. Wolfram + + Last Modified + ------------- + 03/14/2017 + """ + + # **** Initial settings **** + # read parameters from config file + inDirectory = config.get('input', 'baseDirectory') + + namelistFileName = config.get('input', 'oceanNamelistFileName') + namelist = NameList(namelistFileName, path=inDirectory) + + streamsFileName = config.get('input', 'oceanStreamsFileName') + streams = StreamsFile(streamsFileName, streamsdir=inDirectory) + + timeSeriesStatsMonthlyAnalysisMemberFlag = namelist.get( + 'config_am_timeseriesstatsmonthly_enable') + if timeSeriesStatsMonthlyAnalysisMemberFlag == '.false.': + raise RuntimeError('*** MPAS-Analysis relies on the existence of monthly\n' + '*** averaged files: make sure to enable the\n' + '*** timeSeriesStatsMonthly analysis member!') + + # Get a list of timeSeriesStats output files from the streams file, + # reading only those that are between the start and end dates + # First a list necessary for the streamfunctionMOC climatology + streamName = streams.find_stream(streamMap['timeSeriesStats']) + startDateClimo = config.get('climatology', 'startDate') + endDateClimo = config.get('climatology', 'endDate') + calendar = namelist.get('config_calendar_type') + inputFilesClimo = streams.readpath(streamName, + startDate=startDateClimo, + endDate=endDateClimo, + calendar=calendar) + print '\n List of files for climatologies:\n{} through {}'.format( + inputFilesClimo[0], inputFilesClimo[-1]) + startYearClimo = config.getint('climatology', 'startYear') + endYearClimo = config.getint('climatology', 'endYear') + # Create dictionary to store Climo related variables + dictClimo = {'inputFilesClimo': inputFilesClimo, + 'startDateClimo': startDateClimo, + 'endDateClimo': endDateClimo, + 'startYearClimo': startYearClimo, + 'endYearClimo': endYearClimo} + + # Then a list necessary for the streamfunctionMOC Atlantic timeseries + startDateTseries = config.get('timeSeries', 'startDate') + endDateTseries = config.get('timeSeries', 'endDate') + inputFilesTseries = streams.readpath(streamName, + startDate=startDateTseries, + endDate=endDateTseries, + calendar=calendar) + print '\n List of files for timeSeries:\n{} through {}'.format( + inputFilesTseries[0], inputFilesTseries[-1]) + startYearTseries = config.getint('timeSeries', 'startYear') + endYearTseries = config.getint('timeSeries', 'endYear') + # Create dictionary to store Tseries related variables + dictTseries = {'inputFilesTseries': inputFilesTseries, + 'startDateTseries': startDateTseries, + 'endDateTseries': endDateTseries, + 'startYearTseries': startYearTseries, + 'endYearTseries': endYearTseries} + + sectionName = 'streamfunctionMOC' + regionNames = config.getExpression(sectionName, 'regionNames') + + # **** Compute MOC **** + mocAnalysisMemberFlag = namelist.get('config_am_mocstreamfunction_enable') + # Check whether MOC Analysis Member is enabled + if mocAnalysisMemberFlag == '.true.': + # Add a moc_analisysMember_processing + print '*** MOC Analysis Member is on ***' + # (mocDictClimo, mocDictTseries) = _compute_moc_analysismember(config, + # streams, calendar, sectionName, dictClimo, dictTseries) + else: + (mocDictClimo, mocDictTseries) = _compute_moc_postprocess( + config, streams, variableMap, calendar, sectionName, regionNames, + dictClimo, dictTseries) + + # **** Plot MOC **** + # Define plotting variables + plotsDirectory = buildConfigFullPath(config, 'output', 'plotsSubdirectory') + mainRunName = config.get('runs', 'mainRunName') + movingAveragePoints = config.getint(sectionName, 'movingAveragePoints') + colorbarLabel = '[Sv]' + xLabel = 'latitude [deg]' + yLabel = 'depth [m]' + + for region in regionNames: + print ' Plot climatological {} MOC...'.format(region) + title = '{} MOC (ANN, years {:04d}-{:04d})\n {}'.format( + region, startYearClimo, endYearClimo, mainRunName) + figureName = '{}/moc{}_{}_years{:04d}-{:04d}.png'.format( + plotsDirectory, region, mainRunName, + startYearClimo, endYearClimo) + contourLevels = config.getExpression(sectionName, + 'contourLevels{}'.format(region), + usenumpyfunc=True) + (colormapName, colorbarLevels) = setup_colormap(config, sectionName, + suffix=region) + + x = mocDictClimo['lat{}'.format(region)]['data'] + y = mocDictClimo['depth']['data'] + z = mocDictClimo['moc{}'.format(region)]['data'] + plot_vertical_section(config, x, y, z, colormapName, colorbarLevels, + contourLevels, colorbarLabel, title, + xLabel, yLabel, figureName) + + # Plot time series + print ' Plot time series of max Atlantic MOC at 26.5N...' + xLabel = 'Time [years]' + yLabel = '[Sv]' + title = 'Max Atlantic MOC at $26.5^\circ$N\n {}'.format(mainRunName) + figureName = '{}/mocTimeseries_{}.png'.format(plotsDirectory, + mainRunName) + + dsmoc = xr.Dataset.from_dict(mocDictTseries) + timeseries_analysis_plot(config, [dsmoc.mocRapidArrayMax], + movingAveragePoints, title, + xLabel, yLabel, figureName, + lineStyles=['k-'], lineWidths=[1.5], + calendar=calendar) + + return # }}} + + +def _compute_moc_postprocess(config, streams, variableMap, calendar, + sectionName, regionNames, dictClimo, + dictTseries): + + '''compute MOC streamfunction at postprocessing''' + + # Load mesh related variables + try: + restartFile = streams.readpath('restart')[0] + except ValueError: + raise IOError('No MPAS-O restart file found: need at least one ' + 'restart file for MOC calculation') + ncFile = netCDF4.Dataset(restartFile, mode='r') + dvEdge = ncFile.variables['dvEdge'][:] + areaCell = ncFile.variables['areaCell'][:] + refBottomDepth = ncFile.variables['refBottomDepth'][:] + latCell = ncFile.variables['latCell'][:] + latCell = latCell * rad_to_deg # convert to degree + ncFile.close() + + nVertLevels = len(refBottomDepth) + refTopDepth = np.zeros(nVertLevels+1) + refTopDepth[1:nVertLevels+1] = refBottomDepth[0:nVertLevels] + refLayerThickness = np.zeros(nVertLevels) + refLayerThickness[0] = refBottomDepth[0] + refLayerThickness[1:nVertLevels] = (refBottomDepth[1:nVertLevels] - + refBottomDepth[0:nVertLevels-1]) + + variableList = ['avgNormalVelocity', + 'avgVertVelocityTop'] + + # Load basin region related variables and save them to dictionary + # NB: The following will need to change with new regional mapping files + regionMaskFiles = config.get(sectionName, 'regionMaskFiles') + if not os.path.exists(regionMaskFiles): + raise IOError('Regional masking file for MOC calculation ' + 'does not exist') + iRegion = 0 + for region in regionNames: + print '\n Reading region and transect mask for {}...'.format(region) + ncFileRegional = netCDF4.Dataset(regionMaskFiles, mode='r') + maxEdgesInTransect = ncFileRegional.dimensions['maxEdgesInTransect'].size + transectEdgeMaskSigns = ncFileRegional.variables['transectEdgeMaskSigns'][:, iRegion] + transectEdgeGlobalIDs = ncFileRegional.variables['transectEdgeGlobalIDs'][iRegion, :] + regionCellMask = ncFileRegional.variables['regionCellMasks'][:, iRegion] + ncFileRegional.close() + iRegion += 1 + + indRegion = np.where(regionCellMask == 1) + dictRegion = { + 'ind{}'.format(region): indRegion, + '{}CellMask'.format(region): regionCellMask, + 'maxEdgesInTransect{}'.format(region): maxEdgesInTransect, + 'transectEdgeMaskSigns{}'.format(region): transectEdgeMaskSigns, + 'transectEdgeGlobalIDs{}'.format(region): transectEdgeGlobalIDs} + # Add Global regionCellMask=1 everywhere to make the algorithm + # for the global moc similar to that of the regional moc + dictRegion['GlobalCellMask'] = np.ones(np.size(latCell)) + regionNames[:0] = ['Global'] + + # Compute and plot annual climatology of MOC streamfunction + print '\n Compute and/or plot post-processed MOC climatological '\ + 'streamfunction...' + simulationStartTime = get_simulation_start_time(streams) + outputDirectory = buildConfigFullPath(config, 'output', + 'mpasClimatologySubdirectory') + try: + os.makedirs(outputDirectory) + except OSError: + pass + outputFileClimo = '{}/mocStreamfunction_years{:04d}-{:04d}.nc'.format( + outputDirectory, dictClimo['startYearClimo'], + dictClimo['endYearClimo']) + if not os.path.exists(outputFileClimo): + print ' Load data...' + ds = open_multifile_dataset( + fileNames=dictClimo['inputFilesClimo'], + calendar=calendar, + simulationStartTime=simulationStartTime, + timeVariableName='Time', + variableList=variableList, + variableMap=variableMap, + startDate=dictClimo['startDateClimo'], + endDate=dictClimo['endDateClimo']) + + # Compute annual climatology + annualClimatology = climatology.compute_annual_climatology( + ds, calendar) + + horizontalVel = annualClimatology.avgNormalVelocity + verticalVel = annualClimatology.avgVertVelocityTop + # Convert to numpy + horizontalVel = horizontalVel.values + verticalVel = verticalVel.values + velArea = verticalVel * areaCell[:, np.newaxis] + + # Create dictionary for MOC climatology (NB: need this form + # in order to convert it to xarray dataset later in the script) + mocDictClimo = {'depth': {'dims': ('nz'), 'data': refTopDepth}} + for region in regionNames: + print ' Compute {} MOC...'.format(region) + print ' Compute transport through region southern transect...' + if region == 'Global': + transportZ = np.zeros(nVertLevels) + else: + maxEdgesInTransect = dictRegion['maxEdgesInTransect{}'.format(region)] + transectEdgeGlobalIDs = dictRegion['transectEdgeGlobalIDs{}'.format(region)] + transectEdgeMaskSigns = dictRegion['transectEdgeMaskSigns{}'.format(region)] + transportZ = _compute_transport(maxEdgesInTransect, + transectEdgeGlobalIDs, + transectEdgeMaskSigns, + nVertLevels, dvEdge, + refLayerThickness, + horizontalVel) + + regionCellMask = dictRegion['{}CellMask'.format(region)] + latBinSize = config.getExpression(sectionName, + 'latBinSize{}'.format(region)) + if region == 'Global': + latBins = np.arange(-90.0, 90.1, latBinSize) + else: + indRegion = dictRegion['ind{}'.format(region)] + latBins = latCell[indRegion] + latBins = np.arange(np.amin(latBins), + np.amax(latBins)+latBinSize, + latBinSize) + mocTop = _compute_moc(latBins, nVertLevels, latCell, + regionCellMask, transportZ, velArea) + + # Store computed MOC to dictionary + mocDictClimo['lat{}'.format(region)] = { + 'dims': ('nx{}'.format(region)), 'data': latBins} + mocDictClimo['moc{}'.format(region)] = { + 'dims': ('nz', 'nx{}'.format(region)), 'data': mocTop} + + # Save to file + print ' Save global and regional MOC to file...' + ncFile = netCDF4.Dataset(outputFileClimo, mode='w') + # create dimensions + ncFile.createDimension('nz', len(refTopDepth)) + for region in regionNames: + latBins = mocDictClimo['lat{}'.format(region)]['data'] + mocTop = mocDictClimo['moc{}'.format(region)]['data'] + ncFile.createDimension('nx{}'.format(region), len(latBins)) + # create variables + x = ncFile.createVariable('lat{}'.format(region), 'f4', + ('nx{}'.format(region),)) + x.description = 'latitude bins for MOC {}'\ + ' streamfunction'.format(region) + x.units = 'degrees (-90 to 90)' + y = ncFile.createVariable('moc{}'.format(region), 'f4', + ('nz', 'nx{}'.format(region))) + y.description = 'MOC {} streamfunction, annual'\ + ' climatology'.format(region) + y.units = 'Sv (10^6 m^3/s)' + # save variables + x[:] = latBins + y[:, :] = mocTop + depth = ncFile.createVariable('depth', 'f4', ('nz',)) + depth.description = 'depth' + depth.units = 'meters' + depth[:] = refTopDepth + ncFile.close() + else: + # Read from file + print ' Read previously computed MOC streamfunction from file...' + ncFile = netCDF4.Dataset(outputFileClimo, mode='r') + refTopDepth = ncFile.variables['depth'][:] + mocDictClimo = {'depth': {'dims': ('nz'), 'data': refTopDepth}} + for region in regionNames: + latBins = ncFile.variables['lat{}'.format(region)][:] + mocTop = ncFile.variables['moc{}'.format(region)][:, :] + mocDictClimo['lat{}'.format(region)] = { + 'dims': ('nx{}'.format(region)), 'data': latBins} + mocDictClimo['moc{}'.format(region)] = { + 'dims': ('nz', 'nx{}'.format(region)), 'data': mocTop} + ncFile.close() + + # Compute and plot time series of Atlantic MOC at 26.5N (RAPID array) + print '\n Compute and/or plot post-processed Atlantic MOC '\ + 'time series...' + print ' Load data...' + ds = open_multifile_dataset(fileNames=dictTseries['inputFilesTseries'], + calendar=calendar, + simulationStartTime=simulationStartTime, + timeVariableName='Time', + variableList=variableList, + variableMap=variableMap, + startDate=dictTseries['startDateTseries'], + endDate=dictTseries['endDateTseries']) + latAtlantic = mocDictClimo['latAtlantic']['data'] + dLat = latAtlantic - 26.5 + indlat26 = np.where(dLat == np.amin(np.abs(dLat))) + + startYear = days_to_datetime(ds.Time.min()).year + endDate = days_to_datetime(ds.Time.max()) + endYear = endDate.year + if endDate.month < 12: + endYear -= 1 + nyears = (endYear - startYear + 1) + + Time = np.zeros(nyears*12) + mocAtlantic26 = np.zeros(nyears*12) + maxEdgesInTransect = dictRegion['maxEdgesInTransectAtlantic'] + transectEdgeGlobalIDs = dictRegion['transectEdgeGlobalIDsAtlantic'] + transectEdgeMaskSigns = dictRegion['transectEdgeMaskSignsAtlantic'] + regionCellMask = dictRegion['AtlanticCellMask'] + + outputDirectory = buildConfigFullPath(config, 'output', + 'timeseriesSubdirectory') + try: + os.makedirs(outputDirectory) + except OSError: + pass + + ncount = 0 + for ny in range(nyears): + year = dictTseries['startYearTseries'] + ny + print ' year = ', year + outputFileTseries = '{}/mocTimeSeries_year{:04d}.nc'.format( + outputDirectory, year) + if not os.path.exists(outputFileTseries): + for nm in range(12): + print ' month = ', nm + 1 + horizontalVel = ds.avgNormalVelocity[ncount, :, :].values + verticalVel = ds.avgVertVelocityTop[ncount, :, :].values + velArea = verticalVel * areaCell[:, np.newaxis] + transportZ = _compute_transport(maxEdgesInTransect, + transectEdgeGlobalIDs, + transectEdgeMaskSigns, + nVertLevels, dvEdge, + refLayerThickness, + horizontalVel) + mocTop = _compute_moc(latAtlantic, nVertLevels, latCell, + regionCellMask, transportZ, velArea) + mocAtlantic26[ncount] = np.amax(mocTop[:, indlat26]) + ncount += 1 + + Time[ncount-12:ncount] = ds.Time[ncount-12:ncount].values + + # Save to file + print ' save to yearly file...' + ncFile = netCDF4.Dataset(outputFileTseries, mode='w') + # create dimension + ncFile.createDimension('Time', 0) + # create variables + Timevar = ncFile.createVariable('Time', 'f4', ('Time',)) + Timevar.description = 'days since start of simulation' + mocAtlantic26var = ncFile.createVariable('mocRapidArrayMax', + 'f4', ('Time',)) + mocAtlantic26var.description = 'Max MOC Atlantic streamfunction'\ + ' nearest to RAPID Array latitude'\ + ' (26.5N)' + mocAtlantic26var.units = 'Sv (10^6 m^3/s)' + # save variables + Timevar[:] = Time[ncount-12:ncount] + mocAtlantic26var[:] = mocAtlantic26[ncount-12:ncount] + ncFile.close() + else: + ncount += 12 + # Read from file + print ' Read previously computed Atlantic MOC at 26.5N'\ + ' from file...' + ncFile = netCDF4.Dataset(outputFileTseries, mode='r') + Time[ncount-12:ncount] = ncFile.variables['Time'][:] + mocAtlantic26[ncount-12:ncount] = ncFile.variables['mocRapidArrayMax'][:] + ncFile.close() + + # Store max Atlantic MOC at 26.5 to dictionary + mocDictTseries = {'Time': {'dims': ('Time'), 'data': Time}, + 'mocRapidArrayMax': {'dims': ('Time'), + 'data': mocAtlantic26}} + + return (mocDictClimo, mocDictTseries) + + +# def _compute_moc_analysismember(config): +# +# return (mocDictClimo, mocDictTseries) + + +def _compute_transport(maxEdgesInTransect, transectEdgeGlobalIDs, + transectEdgeMaskSigns, nz, dvEdge, refLayerThickness, + horizontalVel): + + '''compute mass transport across southern transect of ocean basin''' + + transportZEdge = np.zeros([nz, maxEdgesInTransect]) + for i in range(maxEdgesInTransect): + if transectEdgeGlobalIDs[i] == 0: + break + # subtract 1 because of python 0-indexing + iEdge = transectEdgeGlobalIDs[i] - 1 + transportZEdge[:, i] = horizontalVel[iEdge, :] * \ + transectEdgeMaskSigns[iEdge, np.newaxis] * \ + dvEdge[iEdge, np.newaxis] * \ + refLayerThickness[np.newaxis, :] + transportZ = transportZEdge.sum(axis=1) + return transportZ + + +def _compute_moc(latBins, nz, latCell, regionCellMask, transportZ, velArea): + + '''compute meridionally integrated MOC streamfunction''' + + mocTop = np.zeros([np.size(latBins), nz+1]) + mocTop[0, range(1, nz+1)] = transportZ.cumsum() + for iLat in range(1, np.size(latBins)): + indlat = np.logical_and(np.logical_and( + regionCellMask == 1, latCell >= latBins[iLat-1]), + latCell < latBins[iLat]) + mocTop[iLat, :] = mocTop[iLat-1, :] + velArea[indlat, :].sum(axis=0) + # convert m^3/s to Sverdrup + mocTop = mocTop * m3ps_to_Sv + mocTop = mocTop.T + return mocTop + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/ocean_modelvsobs.py b/mpas_analysis/ocean/ocean_modelvsobs.py index 968fe2bbe..96d67b90c 100644 --- a/mpas_analysis/ocean/ocean_modelvsobs.py +++ b/mpas_analysis/ocean/ocean_modelvsobs.py @@ -23,7 +23,8 @@ from ..shared.interpolation import interpolate -from ..shared.plot.plotting import plot_global_comparison +from ..shared.plot.plotting import plot_global_comparison, \ + setup_colormap from ..shared.constants import constants from ..shared.io import NameList, StreamsFile @@ -114,7 +115,6 @@ def ocn_modelvsobs(config, field, streamMap=None, variableMap=None): # get a list of regridded observations files and check if they exist. If # they are all there, we don't have to do anything else with the # observations - obsFileNames = \ {'mld': "{}/holtetalley_mld_climatology.nc".format( observationsDirectory), @@ -253,10 +253,10 @@ def ocn_modelvsobs(config, field, streamMap=None, variableMap=None): else: obsMappingFileName = None - (resultColormap, resultContourValues) = _setup_colormap( - config, sectionName, prefix='result') - (differenceColormap, differenceContourValues) = _setup_colormap( - config, sectionName, prefix='difference') + (colormapResult, colorbarLevelsResult) = setup_colormap( + config, sectionName, suffix='Result') + (colormapDifference, colorbarLevelsDifference) = setup_colormap( + config, sectionName, suffix='Difference') # Interpolate and compute biases for months in outputTimes: @@ -333,10 +333,10 @@ def ocn_modelvsobs(config, field, streamMap=None, variableMap=None): modelOutput, observations, bias, - resultColormap, - resultContourValues, - differenceColormap, - differenceContourValues, + colormapResult, + colorbarLevelsResult, + colormapDifference, + colorbarLevelsDifference, fileout=outFileName, title=title, modelTitle="{}".format(mainRunName), @@ -345,29 +345,4 @@ def ocn_modelvsobs(config, field, streamMap=None, variableMap=None): cbarlabel=unitsLabel) -def _setup_colormap(config, sectionName, prefix): - - '''set up a colormap from the registry''' - - contourValues = config.getExpression(sectionName, - '{}ContourValues'.format(prefix)) - colormap = plt.get_cmap(config.get(sectionName, - '{}Colormap'.format(prefix))) - indices = config.getExpression(sectionName, - '{}ColormapIndices'.format(prefix)) - - # set under/over values based on the first/last indices in the colormap - underColor = colormap(indices[0]) - overColor = colormap(indices[-1]) - if len(contourValues)+1 == len(indices): - # we have 2 extra values for the under/over so make the colormap - # without these values - indices = indices[1:-1] - colormap = cols.ListedColormap(colormap(indices), - '{}Colormap'.format(prefix)) - colormap.set_under(underColor) - colormap.set_over(overColor) - return (colormap, contourValues) - - # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/ocean/variable_stream_map.py b/mpas_analysis/ocean/variable_stream_map.py index 1bccabb7a..26d3180ec 100644 --- a/mpas_analysis/ocean/variable_stream_map.py +++ b/mpas_analysis/ocean/variable_stream_map.py @@ -34,6 +34,16 @@ 'time_avg_avgValueWithinOceanLayerRegion_avgLayerThickness_1', 'timeMonthly_avg_avgValueWithinOceanLayerRegion_avgLayerThickness'] +# MOC +oceanVariableMap['avgNormalVelocity'] = \ + ['time_avg_normalVelocity', + 'time_avg_normalVelocity_1', + 'timeMonthly_avg_normalVelocity'] +oceanVariableMap['avgVertVelocityTop'] = \ + ['time_avg_vertVelocityTop', + 'time_avg_vertVelocityTop_1', + 'timeMonthly_avg_vertVelocityTop'] + # model vs. obs. oceanVariableMap['mld'] = \ ['time_avg_dThreshMLD', diff --git a/mpas_analysis/sea_ice/modelvsobs.py b/mpas_analysis/sea_ice/modelvsobs.py index b8bf42bd4..94bc1cb30 100644 --- a/mpas_analysis/sea_ice/modelvsobs.py +++ b/mpas_analysis/sea_ice/modelvsobs.py @@ -27,7 +27,8 @@ from ..shared.climatology import climatology -from ..shared.plot.plotting import plot_polar_comparison +from ..shared.plot.plotting import plot_polar_comparison, \ + setup_colormap from ..shared.io import NameList, StreamsFile from ..shared.io.utility import buildConfigFullPath @@ -247,26 +248,12 @@ def _compute_and_plot_concentration(config, monthlyClimatology, else: plotProjection = 'spstere' - resultConcContourValues = config.getExpression( + (colormapResult, colorbarLevelsResult) = setup_colormap(config, 'regriddedSeaIceConcThick', - 'resultConc{}ContourValues'.format(season)) - resultColormap = plt.get_cmap(config.get( - 'regriddedSeaIceConcThick', 'resultColormap')) - resultColormapIndices = config.getExpression( - 'regriddedSeaIceConcThick', 'resultColormapIndices') - resultColormap = cols.ListedColormap( - resultColormap(resultColormapIndices), "resultColormap") - - differenceConcContourValues = config.getExpression( + suffix='ConcResult{}'.format(season)) + (colormapDifference, colorbarLevelsDifference) = setup_colormap(config, 'regriddedSeaIceConcThick', - 'differenceConc{}ContourValues'.format(season)) - differenceColormap = plt.get_cmap(config.get( - 'regriddedSeaIceConcThick', 'differenceColormap')) - differenceColormapIndices = config.getExpression( - 'regriddedSeaIceConcThick', 'differenceColormapIndices') - differenceColormap = cols.ListedColormap( - differenceColormap(differenceColormapIndices), - "differenceColormap") + suffix='ConcDifference{}'.format(season)) referenceLongitude = config.getfloat( 'regriddedSeaIceConcThick', @@ -319,10 +306,10 @@ def _compute_and_plot_concentration(config, monthlyClimatology, iceConcentration, obsIceConcentration, difference, - resultColormap, - resultConcContourValues, - differenceColormap, - differenceConcContourValues, + colormapResult, + colorbarLevelsResult, + colormapDifference, + colorbarLevelsDifference, title=title, fileout=fileout, plotProjection=plotProjection, @@ -437,26 +424,12 @@ def _compute_and_plot_thickness(config, monthlyClimatology, for hemisphere in ['NH', 'SH']: - resultThickContourValues = config.getExpression( + (colormapResult, colorbarLevelsResult) = setup_colormap(config, 'regriddedSeaIceConcThick', - 'resultThick{}ContourValues'.format(hemisphere)) - resultColormap = plt.get_cmap(config.get( - 'regriddedSeaIceConcThick', 'resultColormap')) - resultColormapIndices = config.getExpression( - 'regriddedSeaIceConcThick', 'resultColormapIndices') - resultColormap = cols.ListedColormap( - resultColormap(resultColormapIndices), "resultColormap") - - differenceThickContourValues = config.getExpression( + suffix='ThickResult{}'.format(hemisphere)) + (colormapDifference, colorbarLevelsDifference) = setup_colormap(config, 'regriddedSeaIceConcThick', - 'differenceThick{}ContourValues'.format(hemisphere)) - differenceColormap = plt.get_cmap(config.get( - 'regriddedSeaIceConcThick', 'differenceColormap')) - differenceColormapIndices = config.getExpression( - 'regriddedSeaIceConcThick', 'differenceColormapIndices') - differenceColormap = cols.ListedColormap( - differenceColormap(differenceColormapIndices), - "differenceColormap") + suffix='ThickDifference{}'.format(hemisphere)) referenceLongitude = config.getfloat( 'regriddedSeaIceConcThick', @@ -517,10 +490,10 @@ def _compute_and_plot_thickness(config, monthlyClimatology, iceThickness, obsIceThickness, difference, - resultColormap, - resultThickContourValues, - differenceColormap, - differenceThickContourValues, + colormapResult, + colorbarLevelsResult, + colormapDifference, + colorbarLevelsDifference, title=title, fileout=fileout, plotProjection=plotProjection, diff --git a/mpas_analysis/sea_ice/timeseries.py b/mpas_analysis/sea_ice/timeseries.py index f079bdf11..4bb73ffc9 100644 --- a/mpas_analysis/sea_ice/timeseries.py +++ b/mpas_analysis/sea_ice/timeseries.py @@ -1,6 +1,7 @@ import xarray as xr -from ..shared.plot.plotting import timeseries_analysis_plot, timeseries_analysis_plot_polar +from ..shared.plot.plotting import timeseries_analysis_plot, \ + timeseries_analysis_plot_polar from ..shared.io import NameList, StreamsFile from ..shared.io.utility import buildConfigFullPath diff --git a/mpas_analysis/shared/climatology/climatology.py b/mpas_analysis/shared/climatology/climatology.py index 1a3e0f13a..8aba57cfd 100644 --- a/mpas_analysis/shared/climatology/climatology.py +++ b/mpas_analysis/shared/climatology/climatology.py @@ -475,6 +475,37 @@ def compute_seasonal_climatology(monthlyClimatology, monthValues, return seasonalClimatology +def compute_annual_climatology(ds, calendar): + """ + Compute an annual climatology data set from a data set with Time expressed + as days since 0001-01-01 with the given calendar. + + Parameters + ---------- + ds : an xarray data set with a 'Time' coordinate expressed as days since + 0001-01-01 + + calendar: {'gregorian', 'gregorian_noleap'} + The name of one of the calendars supported by MPAS cores + + Returns + ------- + annualClimatology : an xarray data set containing annual climatologies + of all variables in ds (time coordinate is squashed) + + Authors + ------- + Milena Veneziani + + Last Modified + ------------- + 02/28/2017 + """ + monthlyClimatology = compute_monthly_climatology(ds, calendar) + annualClimatology = monthlyClimatology.mean('month') + return annualClimatology + + def _get_comparison_lat_lon(comparisonLatRes, comparisonLonRes): ''' Returns the lat and lon arrays defining the corners of the comparison diff --git a/mpas_analysis/shared/constants/constants.py b/mpas_analysis/shared/constants/constants.py index f32e4412e..ad608fdfc 100644 --- a/mpas_analysis/shared/constants/constants.py +++ b/mpas_analysis/shared/constants/constants.py @@ -3,8 +3,13 @@ """ Constants that are common to all analysis tasks -Luke Van Roekel, Xylar Asay-Davis -02/26/2017 +Authors +------- +Luke Van Roekel, Xylar Asay-Davis, Milena Veneziani + +Last modified +------------- +03/15/2017 """ # set parameters for default climatology comparison grid @@ -27,4 +32,13 @@ abrevMonthNames = ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"] +# conversion factor from m^3/s to Sverdrups +m3ps_to_Sv = 1e-6 + +# conversion factor from radians to degrees +rad_to_deg = 180./np.pi + +# conversion factor from degrees to radians +deg_to_rad = np.pi/180. + # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/mpas_analysis/shared/plot/plotting.py b/mpas_analysis/shared/plot/plotting.py index 03602acd0..8a534b356 100644 --- a/mpas_analysis/shared/plot/plotting.py +++ b/mpas_analysis/shared/plot/plotting.py @@ -1,11 +1,12 @@ """ -Plotting utilities, including routine for plotting: +Plotting utilities, including routines for plotting: * time series (and comparing with reference data sets) * regridded horizontal fields (and comparing with reference data sets) + * vertical sections on native grid Authors ------- -Xylar Asay-Davis +Xylar Asay-Davis, Milena Veneziani Last Modified ------------- @@ -13,8 +14,9 @@ """ import matplotlib.pyplot as plt -import pandas as pd +import matplotlib.colors as cols import xarray as xr +import pandas as pd from mpl_toolkits.basemap import Basemap import matplotlib.colors as cols from matplotlib.ticker import FuncFormatter, FixedLocator @@ -74,7 +76,7 @@ def timeseries_analysis_plot(config, dsvalues, N, title, xlabel, ylabel, Authors ------- - Xylar Asay-Davis + Xylar Asay-Davis, Milena Veneziani Last Modified ------------- @@ -96,6 +98,15 @@ def timeseries_analysis_plot(config, dsvalues, N, title, xlabel, ylabel, lineStyles[dsIndex], linewidth=lineWidths[dsIndex]) + ax = plt.gca() + + # Add a y=0 line if y ranges between positive and negative values + yaxLimits = ax.get_ylim() + if yaxLimits[0]*yaxLimits[1] < 0: + indgood = np.where(np.logical_not(np.isnan(mean))) + x = mean['Time'][indgood] + plt.plot(x, np.zeros(np.size(x)), 'k-', linewidth=1.2) + start = days_to_datetime(np.amin(minDays), calendar=calendar) end = days_to_datetime(np.amax(maxDays), calendar=calendar) @@ -114,7 +125,6 @@ def timeseries_analysis_plot(config, dsvalues, N, title, xlabel, ylabel, formatterFun = partial(_date_tick, calendar=calendar, includeMonth=True) - ax = plt.gca() ax.xaxis.set_major_locator(FixedLocator(major, maxXTicks)) ax.xaxis.set_major_formatter(FuncFormatter(formatterFun)) @@ -142,10 +152,10 @@ def timeseries_analysis_plot(config, dsvalues, N, title, xlabel, ylabel, plt.close() - def timeseries_analysis_plot_polar(config, dsvalues, N, title, - fileout, lineStyles, lineWidths, calendar, - titleFontSize=None, figsize=(15, 6), dpi=300): + fileout, lineStyles, lineWidths, + calendar, titleFontSize=None, + figsize=(15, 6), dpi=300): """ Plots the list of time series data sets on a polar plot and stores @@ -254,60 +264,85 @@ def plot_polar_comparison( cmapDiff, clevsDiff, fileout, - title = None, - plotProjection = "npstere", - latmin = 50.0, - lon0 = 0, - modelTitle = "Model", - obsTitle = "Observations", - diffTitle = "Model-Observations", - cbarlabel = "units", - titleFontSize = None, - figsize = (8,22), - dpi = 300): + title=None, + plotProjection='npstere', + latmin=50.0, + lon0=0, + modelTitle='Model', + obsTitle='Observations', + diffTitle='Model-Observations', + cbarlabel='units', + titleFontSize=None, + figsize=(8, 22), + dpi=300): """ Plots a data set around either the north or south pole. - config is an instance of ConfigParser + Parameters + ---------- + config : instance of ConfigParser + the configuration, containing a [plot] section with options that + control plotting - Lons and Lats are the longitude and latitude arrays + Lons, Lats : float arrays + longitude and latitude arrays - modelArray and obsArray are the model and observational data sets + modelArray, obsArray : float arrays + model and observational data sets - diffArray is the difference between modelArray and obsArray + diffArray : float array + difference between modelArray and obsArray - cmapModelObs is the colormap of model and observations + cmapModelObs : str + colormap of model and observations panel - clevsModleObs are contour values for model and observations + clevsModelObs : int array + colorbar values for model and observations panel - cmapDiff is the colormap of difference + cmapDiff : str + colormap of difference (bias) panel - clevsDiff are contour values fordifference + clevsDiff : int array + colorbar values for difference (bias) panel - fileout is the file name to be written + fileout : str + the file name to be written - title is a string with the subtitle of the plot + title : str, optional + the subtitle of the plot - plotProjection is a Basemap projection for the plot + plotProjection : str, optional + Basemap projection for the plot - modelTitle is the title of the model plot + modelTitle : str, optional + title of the model panel - obsTitle is the title of the observations plot + obsTitle : str, optional + title of the observations panel - diffTitle is the title of the difference plot + diffTitle : str, optional + title of the difference (bias) panel - cbarlabel is the label on the colorbar + cbarlabel : str, optional + label on the colorbar - titleFontSize is the size of the title font + titleFontSize : int, optional + size of the title font - figsize is the size of the figure in inches + figsize : tuple of float, optional + the size of the figure in inches - dpi is the number of dots per inch of the figure + dpi : int, optional + the number of dots per inch of the figure - Author: Xylar Asay-Davis + Authors + ------- + Xylar Asay-Davis, Milena Veneziani - Last Modified: 02/02/2017 + Last Modified + ------------- + 03/17/2017 """ # set up figure @@ -316,53 +351,64 @@ def plot_polar_comparison( if titleFontSize is None: titleFontSize = config.get('plot', 'titleFontSize') title_font = {'size': titleFontSize, - 'color':config.get('plot', 'titleFontColor'), - 'weight':config.get('plot', 'titleFontWeight')} + 'color': config.get('plot', 'titleFontColor'), + 'weight': config.get('plot', 'titleFontWeight')} fig.suptitle(title, y=0.95, **title_font) - axis_font = {'size':config.get('plot', 'axisFontSize')} + axis_font = {'size': config.get('plot', 'axisFontSize')} - m = Basemap(projection=plotProjection,boundinglat=latmin,lon_0=lon0,resolution='l') - x, y = m(Lons, Lats) # compute map proj coordinates + m = Basemap(projection=plotProjection, boundinglat=latmin, + lon_0=lon0, resolution='l') + x, y = m(Lons, Lats) # compute map proj coordinates normModelObs = cols.BoundaryNorm(clevsModelObs, cmapModelObs.N) normDiff = cols.BoundaryNorm(clevsDiff, cmapDiff.N) - plt.subplot(3,1,1) + plt.subplot(3, 1, 1) plt.title(modelTitle, y=1.06, **axis_font) m.drawcoastlines() - m.fillcontinents(color='grey',lake_color='white') - m.drawparallels(np.arange(-80.,81.,10.)) - m.drawmeridians(np.arange(-180.,181.,20.),labels=[True,True,True,True]) - cs = m.contourf(x,y,modelArray,cmap=cmapModelObs,norm=normModelObs,spacing='uniform',levels=clevsModelObs) - cbar = m.colorbar(cs,location='right',pad="15%",spacing='uniform',ticks=clevsModelObs,boundaries=clevsModelObs) + m.fillcontinents(color='grey', lake_color='white') + m.drawparallels(np.arange(-80., 81., 10.)) + m.drawmeridians(np.arange(-180., 181., 20.), + labels=[True, True, True, True]) + cs = m.contourf(x, y, modelArray, cmap=cmapModelObs, norm=normModelObs, + spacing='uniform', levels=clevsModelObs) + cbar = m.colorbar(cs, location='right', pad="15%", spacing='uniform', + ticks=clevsModelObs, boundaries=clevsModelObs) cbar.set_label(cbarlabel) - plt.subplot(3,1,2) + 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.,81.,10.)) - m.drawmeridians(np.arange(-180.,181.,20.),labels=[True,True,True,True]) - cs = m.contourf(x,y,obsArray,cmap=cmapModelObs,norm=normModelObs,spacing='uniform',levels=clevsModelObs) - cbar = m.colorbar(cs,location='right',pad="15%",spacing='uniform',ticks=clevsModelObs,boundaries=clevsModelObs) + m.fillcontinents(color='grey', lake_color='white') + m.drawparallels(np.arange(-80., 81., 10.)) + m.drawmeridians(np.arange(-180., 181., 20.), + labels=[True, True, True, True]) + cs = m.contourf(x, y, obsArray, cmap=cmapModelObs, norm=normModelObs, + spacing='uniform', levels=clevsModelObs) + cbar = m.colorbar(cs, location='right', pad="15%", spacing='uniform', + ticks=clevsModelObs, boundaries=clevsModelObs) cbar.set_label(cbarlabel) - plt.subplot(3,1,3) + 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.,81.,10.)) - m.drawmeridians(np.arange(-180.,181.,20.),labels=[True,True,True,True]) - cs = m.contourf(x,y,diffArray,cmap=cmapDiff,norm=normDiff,spacing='uniform',levels=clevsDiff) - cbar = m.colorbar(cs,location='right',pad="15%",spacing='uniform',ticks=clevsDiff,boundaries=clevsModelObs) + m.fillcontinents(color='grey', lake_color='white') + m.drawparallels(np.arange(-80., 81., 10.)) + m.drawmeridians(np.arange(-180., 181., 20.), + labels=[True, True, True, True]) + cs = m.contourf(x, y, diffArray, cmap=cmapDiff, norm=normDiff, + spacing='uniform', levels=clevsDiff) + cbar = m.colorbar(cs, location='right', pad="15%", spacing='uniform', + ticks=clevsDiff, boundaries=clevsModelObs) cbar.set_label(cbarlabel) if (fileout is not None): - plt.savefig(fileout,dpi=dpi,bbox_inches='tight',pad_inches=0.1) + plt.savefig(fileout, dpi=dpi, bbox_inches='tight', pad_inches=0.1) - if not config.getboolean('plot','displayToScreen'): + if not config.getboolean('plot', 'displayToScreen'): plt.close() + def plot_global_comparison( config, Lons, @@ -375,107 +421,145 @@ def plot_global_comparison( cmapDiff, clevsDiff, fileout, - title = None, - modelTitle = "Model", - obsTitle = "Observations", - diffTitle = "Model-Observations", - cbarlabel = "units", - titleFontSize = None, - figsize = (8,12), - dpi = 300): + title=None, + modelTitle='Model', + obsTitle='Observations', + diffTitle='Model-Observations', + cbarlabel='units', + titleFontSize=None, + figsize=(8, 12), + dpi=300): """ Plots a data set as a longitude/latitude map. - config is an instance of ConfigParser + Parameters + ---------- + config : instance of ConfigParser + the configuration, containing a [plot] section with options that + control plotting - Lons and Lats are the longitude and latitude arrays + Lons, Lats : float arrays + longitude and latitude arrays - modelArray and obsArray are the model and observational data sets + modelArray, obsArray : float arrays + model and observational data sets - diffArray is the difference between modelArray and obsArray + diffArray : float array + difference between modelArray and obsArray - cmapModelObs is the colormap of model and observations + cmapModelObs : str + colormap of model and observations panel - clevsModleObs are contour values for model and observations + clevsModelObs : int array + colorbar values for model and observations panel - cmapDiff is the colormap of difference + cmapDiff : str + colormap of difference (bias) panel - clevsDiff are contour values fordifference + clevsDiff : int array + colorbar values for difference (bias) panel - fileout is the file name to be written + fileout : str + the file name to be written - title is a string with the subtitle of the plot + title : str, optional + the subtitle of the plot - modelTitle is the title of the model plot + modelTitle : str, optional + title of the model panel - obsTitle is the title of the observations plot + obsTitle : str, optional + title of the observations panel - diffTitle is the title of the difference plot + diffTitle : str, optional + title of the difference (bias) panel - cbarlabel is the label on the colorbar + cbarlabel : str, optional + label on the colorbar - titleFontSize is the size of the title font + titleFontSize : int, optional + size of the title font - figsize is the size of the figure in inches + figsize : tuple of float, optional + the size of the figure in inches - dpi is the number of dots per inch of the figure + dpi : int, optional + the number of dots per inch of the figure - Author: Xylar Asay-Davis + Authors + ------- + Xylar Asay-Davis, Milena Veneziani - Last Modified: 03/09/2017 + Last Modified + ------------- + 03/13/2017 """ + # set up figure fig = plt.figure(figsize=figsize, dpi=dpi) if (title is not None): if titleFontSize is None: titleFontSize = config.get('plot', 'titleFontSize') title_font = {'size': titleFontSize, - 'color':config.get('plot', 'titleFontColor'), - 'weight':config.get('plot', 'titleFontWeight')} + 'color': config.get('plot', 'titleFontColor'), + 'weight': config.get('plot', 'titleFontWeight')} fig.suptitle(title, y=0.95, **title_font) - axis_font = {'size':config.get('plot', 'axisFontSize')} + axis_font = {'size': config.get('plot', 'axisFontSize')} - m = Basemap(projection='cyl',llcrnrlat=-85,urcrnrlat=86,llcrnrlon=-180,urcrnrlon=181,resolution='l') - x, y = m(Lons, Lats) # compute map proj coordinates + m = Basemap(projection='cyl', llcrnrlat=-85, urcrnrlat=86, llcrnrlon=-180, + urcrnrlon=181, resolution='l') + x, y = m(Lons, Lats) # compute map proj coordinates normModelObs = cols.BoundaryNorm(clevsModelObs, cmapModelObs.N) normDiff = cols.BoundaryNorm(clevsDiff, cmapDiff.N) - plt.subplot(3,1,1) + plt.subplot(3, 1, 1) plt.title(modelTitle, 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.,60.),labels=[False,False,False,True]) - cs = m.contourf(x,y,modelArray,cmap=cmapModelObs,norm=normModelObs,spacing='uniform',levels=clevsModelObs,extend='both') - cbar = m.colorbar(cs,location='right',pad="5%",spacing='uniform',ticks=clevsModelObs,boundaries=clevsModelObs) + 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., 60.), + labels=[False, False, False, True]) + cs = m.contourf(x, y, modelArray, cmap=cmapModelObs, norm=normModelObs, + spacing='uniform', levels=clevsModelObs, extend='both') + cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform', + ticks=clevsModelObs, boundaries=clevsModelObs) cbar.set_label(cbarlabel) - plt.subplot(3,1,2) + 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,spacing='uniform',levels=clevsModelObs,extend='both') - cbar = m.colorbar(cs,location='right',pad="5%",spacing='uniform',ticks=clevsModelObs,boundaries=clevsModelObs) + 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, + spacing='uniform', levels=clevsModelObs, extend='both') + cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform', + ticks=clevsModelObs, boundaries=clevsModelObs) cbar.set_label(cbarlabel) - plt.subplot(3,1,3) + 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,spacing='uniform',levels=clevsDiff,extend='both') - cbar = m.colorbar(cs,location='right',pad="5%",spacing='uniform',ticks=clevsDiff,boundaries=clevsDiff) + 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, + spacing='uniform', levels=clevsDiff, extend='both') + cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform', + ticks=clevsDiff, boundaries=clevsModelObs) cbar.set_label(cbarlabel) if (fileout is not None): - plt.savefig(fileout,dpi=dpi,bbox_inches='tight',pad_inches=0.1) + plt.savefig(fileout, dpi=dpi, bbox_inches='tight', pad_inches=0.1) - if not config.getboolean('plot','displayToScreen'): + if not config.getboolean('plot', 'displayToScreen'): plt.close() @@ -486,3 +570,169 @@ def _date_tick(days, pos, calendar='gregorian', includeMonth=True): return '{:04d}-{:02d}'.format(date.year, date.month) else: return '{:04d}'.format(date.year) + + +def plot_vertical_section( + config, + xArray, + depthArray, + fieldArray, + colormapName, + colorbarLevels, + contourLevels, + colorbarLabel=None, + title=None, + xlabel=None, + ylabel=None, + fileout='moc.png', + figsize=(10, 4), + dpi=300): # {{{ + + """ + Plots a data set as a x distance (latitude, longitude, + or spherical distance) vs depth map (vertical section). + + Parameters + ---------- + config : instance of ConfigParser + the configuration, containing a [plot] section with options that + control plotting + + xArray : float array + x array (latitude, longitude, or spherical distance) + + depthArray : float array + depth array [m] + + fieldArray : float array + field array to plot + + colormapName : str + colormap of plot + + colorbarLevels : int array + colorbar levels of plot + + contourLevels : int levels + levels of contours to be drawn + + colorbarLabel : str, optional + label of the colorbar + + title : str, optional + title of plot + + xlabel, ylabel : str, optional + label of x- and y-axis + + fileout : str, optional + the file name to be written + + figsize : tuple of float, optional + size of the figure in inches + + dpi : int, optional + number of dots per inch of the figure + + Authors + ------- + Milena Veneziani, Mark Petersen + + Last Modified + ------------- + 03/13/2017 + """ + + # set up figure + fig = plt.figure(figsize=figsize, dpi=dpi) + + x, y = np.meshgrid(xArray, depthArray) # change to zMid + + normModelObs = cols.BoundaryNorm(colorbarLevels, colormapName.N) + + cs = plt.contourf(x, y, fieldArray, cmap=colormapName, norm=normModelObs, + spacing='uniform', levels=colorbarLevels, extend='both') + plt.contour(x, y, fieldArray, levels=contourLevels[::2], colors='k') + + cbar = plt.colorbar(cs, orientation='vertical', spacing='uniform', + ticks=colorbarLevels, boundaries=colorbarLevels) + if colorbarLabel is not None: + cbar.set_label(colorbarLabel) + + axis_font = {'size': config.get('plot', 'axisFontSize')} + title_font = {'size': config.get('plot', 'titleFontSize'), + 'color': config.get('plot', 'titleFontColor'), + 'weight': config.get('plot', 'titleFontWeight')} + if title is not None: + plt.title(title, **title_font) + if xlabel is not None: + plt.xlabel(xlabel, **axis_font) + if ylabel is not None: + plt.ylabel(ylabel, **axis_font) + + plt.gca().invert_yaxis() + + if (fileout is not None): + plt.savefig(fileout, dpi=dpi, bbox_inches='tight', pad_inches=0.1) + + if not config.getboolean('plot', 'displayToScreen'): + plt.close() + + return # }}} + + +def setup_colormap(config, configSectionName, suffix=''): + + ''' + Set up a colormap from the registry + + Parameters + ---------- + config : instance of ConfigParser + the configuration, containing a [plot] section with options that + control plotting + + configSectionName : str + name of config section + + suffix: str, optional + suffix of colormap related options + + Returns + ------- + colormap : srt + new colormap + + colorbarLevels : int array + colorbar levels + + Authors + ------- + Xylar Asay-Davis, Milena Veneziani + + Last modified + ------------- + 03/17/2017 + ''' + + colormap = plt.get_cmap(config.get(configSectionName, + 'colormapName{}'.format(suffix))) + indices = config.getExpression(configSectionName, + 'colormapIndices{}'.format(suffix)) + colorbarLevels = config.getExpression(configSectionName, + 'colorbarLevels{}'.format(suffix)) + + # set under/over values based on the first/last indices in the colormap + underColor = colormap(indices[0]) + overColor = colormap(indices[-1]) + if len(colorbarLevels)+1 == len(indices): + # we have 2 extra values for the under/over so make the colormap + # without these values + indices = indices[1:-1] + colormap = cols.ListedColormap(colormap(indices), + 'colormapName{}'.format(suffix)) + colormap.set_under(underColor) + colormap.set_over(overColor) + return (colormap, colorbarLevels) + +# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python diff --git a/run_analysis.py b/run_analysis.py index 541942705..244a67e21 100755 --- a/run_analysis.py +++ b/run_analysis.py @@ -163,13 +163,6 @@ def analysis(config): # {{{ # from mpas_analysis.ocean.mht_timeseries import mht_timeseries # mht_timeseries(config) -# if checkGenerate(config, analysisName='timeSeriesMOC', mpasCore='ocean', -# analysisCategory='timeSeries'): -# print "" -# print "Plotting Meridional Overturning Circulation (MOC)..." -# from mpas_analysis.ocean.moc_timeseries import moc_timeseries -# moc_timeseries(config) - if checkGenerate(config, analysisName='regriddedSST', mpasCore='ocean', analysisCategory='regriddedHorizontal'): print "" @@ -194,6 +187,14 @@ def analysis(config): # {{{ ocn_modelvsobs(config, 'sss', streamMap=oceanStreamMap, variableMap=oceanVariableMap) + if checkGenerate(config, analysisName='streamfunctionMOC', mpasCore='ocean', + analysisCategory='streamfunctionMOC'): + print "" + print "Plotting streamfunction of Meridional Overturning Circulation (MOC)..." + from mpas_analysis.ocean.meridional_overturning_circulation import moc_streamfunction + moc_streamfunction(config, streamMap=oceanStreamMap, + variableMap=oceanVariableMap) + # GENERATE SEA-ICE DIAGNOSTICS if checkGenerate(config, analysisName='timeSeriesSeaIceAreaVol', mpasCore='seaIce', analysisCategory='timeSeries'):