diff --git a/tools/analysis/SSS_bias_WOA05.py b/tools/analysis/SSS_bias_WOA05.py index 3ef9cf6230..d9976396e2 100755 --- a/tools/analysis/SSS_bias_WOA05.py +++ b/tools/analysis/SSS_bias_WOA05.py @@ -13,6 +13,7 @@ def run(): parser = argparse.ArgumentParser(description='''Script for plotting annual-average SSS bias.''') parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'salt'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-s','--suptitle', type=str, default='', help='''Super-title for experiment. Default is to read from netCDF file.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') @@ -56,15 +57,18 @@ def main(cmdLineArgs,stream=None): if rootGroup.variables[varName].shape[0]>1: Smod = rootGroup.variables[varName][:,0].mean(axis=0) else: Smod = rootGroup.variables[varName][0,0] + if cmdLineArgs.suptitle != '': suptitle = cmdLineArgs.suptitle + ' ' + cmdLineArgs.label + else: suptitle = rootGroup.title + ' ' + cmdLineArgs.label + ci=m6plot.pmCI(0.125,2.25,.25) if stream == None: stream = cmdLineArgs.outdir+'/SSS_bias_WOA05.png' m6plot.xyplot( Smod - Sobs , x, y, area=area, - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='SSS bias (w.r.t. WOA\'05) [ppt]', + suptitle=suptitle, title='SSS bias (w.r.t. WOA\'05) [ppt]', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=stream) m6plot.xycompare( Smod, Sobs , x, y, area=area, - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1='SSS [ppt]', title2='WOA\'05 SSS [ppt]', clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', diff --git a/tools/analysis/SST_bias_WOA05.py b/tools/analysis/SST_bias_WOA05.py index 6a4b18ee04..1bc6c63541 100755 --- a/tools/analysis/SST_bias_WOA05.py +++ b/tools/analysis/SST_bias_WOA05.py @@ -13,6 +13,7 @@ def run(): parser = argparse.ArgumentParser(description='''Script for plotting annual-average SST bias.''') parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'temp'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-s','--suptitle', type=str, default='', help='''Super-title for experiment. Default is to read from netCDF file.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') @@ -59,15 +60,18 @@ def main(cmdLineArgs,stream=None): if rootGroup.variables[varName].shape[0]>1: Tmod = rootGroup.variables[varName][:,0].mean(axis=0) else: Tmod = rootGroup.variables[varName][0,0] + if cmdLineArgs.suptitle != '': suptitle = cmdLineArgs.suptitle + ' ' + cmdLineArgs.label + else: suptitle = rootGroup.title + ' ' + cmdLineArgs.label + ci=m6plot.pmCI(0.25,4.5,.5) if stream == None: stream = cmdLineArgs.outdir+'/SST_bias_WOA05.png' m6plot.xyplot( Tmod - Tobs , x, y, area=area, - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='SST bias (w.r.t. WOA\'05) [$\degree$C]', + suptitle=suptitle, title='SST bias (w.r.t. WOA\'05) [$\degree$C]', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=stream) m6plot.xycompare( Tmod, Tobs , x, y, area=area, - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1='SST [$\degree$C]', title2='WOA\'05 SST [$\degree$C]', clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', diff --git a/tools/analysis/TS_drift.py b/tools/analysis/TS_drift.py new file mode 100755 index 0000000000..8570b7e6a2 --- /dev/null +++ b/tools/analysis/TS_drift.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python + +import netCDF4 +import numpy +import m6plot +import m6toolbox +import matplotlib.pyplot as plt +import os + +try: import argparse +except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.') + +def run(): + parser = argparse.ArgumentParser(description='''Script for plotting depth vs. time plots of temperature and salinity drift''') + parser.add_argument('annual_directory', type=str, help='''Directory containing annual time series thetao and so xyave files''') + parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-s','--suptitle', type=str, default='', help='''Super-title for experiment. Default is to read from netCDF file.''') + parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') + parser.add_argument('-t','--trange', type=str, default=None, help='''Tuple containing start and end years to plot''') + cmdLineArgs = parser.parse_args() + main(cmdLineArgs) + +def main(cmdLineArgs,stream=None): + rootGroupT = netCDF4.MFDataset( cmdLineArgs.annual_directory + '/*.thetao_xyave.nc' ) + rootGroupS = netCDF4.MFDataset( cmdLineArgs.annual_directory + '/*.so_xyave.nc' ) + if 'thetao_xyave' not in rootGroupT.variables: raise Exception('Could not find "thetao_xyave" files "%s"'%(cmdLineArgs.annual_directory)) + if 'so_xyave' not in rootGroupS.variables: raise Exception('Could not find "so_xyave" files "%s"'%(cmdLineArgs.annual_directory)) + + zt = rootGroupT.variables['zt'][::-1] * -1 + timeT = rootGroupT.variables['time'] + timeS = rootGroupS.variables['time'] + timeT = numpy.array([int(x.year) for x in netCDF4.num2date(timeT[:],timeT.units,calendar=timeT.calendar)]) + timeS = numpy.array([int(x.year) for x in netCDF4.num2date(timeS[:],timeS.units,calendar=timeS.calendar)]) + + if cmdLineArgs.trange != None: + start = list(timeT).index(cmdLineArgs.trange[0]) + end = list(timeT).index(cmdLineArgs.trange[1]) + else: + start = 0 + end = -1 + + variable = rootGroupT.variables['thetao_xyave'] + T = variable[start:end] - variable[start] + T = T[:,::-1] + timeT = timeT[start:end] + + variable = rootGroupS.variables['so_xyave'] + S = variable[start:end] - variable[start] + S = S[:,::-1] + timeS = timeS[start:end] + + if cmdLineArgs.suptitle != '': suptitle = cmdLineArgs.suptitle + ' ' + cmdLineArgs.label + else: suptitle = rootGroupT.title + ' ' + cmdLineArgs.label + + if stream != None: objOut = stream[0] + else: objOut = cmdLineArgs.outdir+'/T_drift.png' + m6plot.ztplot( T, timeT, zt, splitscale=[0., -1000., -6500.], + suptitle=suptitle, title='Potential Temperature [C]', + extend='both', colormap='dunnePM', centered=True, + save=objOut) + + if stream != None: objOut = stream[1] + else: objOut = cmdLineArgs.outdir+'/S_drift.png' + m6plot.ztplot( S, timeS, zt, splitscale=[0., -1000., -6500.], + suptitle=suptitle, title='Salinity [psu]', + extend='both', colormap='dunnePM', centered=True, + save=objOut) + +if __name__ == '__main__': + run() + diff --git a/tools/analysis/m6plot.py b/tools/analysis/m6plot.py index 481c7dfa94..f913b7020f 100644 --- a/tools/analysis/m6plot.py +++ b/tools/analysis/m6plot.py @@ -71,7 +71,10 @@ def xyplot(field, x=None, y=None, area=None, # Choose colormap if nbins is None and (clim is None or len(clim)==2): nbins=35 if colormap is None: colormap = chooseColorMap(sMin, sMax) - cmap, norm, extend = chooseColorLevels(sMin, sMax, colormap, clim=clim, nbins=nbins, extend=extend, logscale=logscale) + if clim is None and sStd>0: + cmap, norm, extend = chooseColorLevels(sMean-2.*sStd, sMean+2.*sStd, colormap, clim=clim, nbins=nbins, extend=extend, logscale=logscale) + else: + cmap, norm, extend = chooseColorLevels(sMin, sMax, colormap, clim=clim, nbins=nbins, extend=extend, logscale=logscale) if axis is None: setFigureSize(aspect, resolution, debug=debug) @@ -81,6 +84,7 @@ def xyplot(field, x=None, y=None, area=None, if interactive: addStatusBar(xCoord, yCoord, maskedField) cb = plt.colorbar(fraction=.08, pad=0.02, extend=extend) if centerlabels and len(clim)>2: cb.set_ticks( 0.5*(clim[:-1]+clim[1:]) ) + elif len(clim)>2: cb.set_ticks( clim ) axis.set_axis_bgcolor(landcolor) plt.xlim( xLims ) plt.ylim( yLims ) @@ -218,7 +222,10 @@ def annotateStats(axis, sMin, sMax, sMean, sStd, sRMS): if npanels in [1,3]: axis = plt.subplot(npanels,1,npanels) if dcolormap is None: dcolormap = chooseColorMap(dMin, dMax) - cmap, norm, extend = chooseColorLevels(dMin, dMax, dcolormap, clim=dlim, nbins=nbins, extend=dextend) + if dlim is None and dStd>0: + cmap, norm, dextend = chooseColorLevels(dMean-2.*dStd, dMean+2.*dStd, dcolormap, clim=dlim, nbins=nbins, extend='both', autocenter=True) + else: + cmap, norm, dextend = chooseColorLevels(dMin, dMax, dcolormap, clim=dlim, nbins=nbins, extend=dextend, autocenter=True) plt.pcolormesh(xCoord, yCoord, maskedField1 - maskedField2, cmap=cmap, norm=norm) if interactive: addStatusBar(xCoord, yCoord, maskedField1 - maskedField2) if dextend is None: dextend = extend @@ -230,7 +237,7 @@ def annotateStats(axis, sMin, sMax, sMean, sStd, sRMS): if len(ylabel+yunits)>0: plt.ylabel(label(ylabel, yunits)) if len(title3)>0: plt.title(title3) - if dRxy is not None: axis.annotate(' r(A,B)=%.5g\n'%(dRxy), xy=(1.0,-0.20), xycoords='axes fraction', verticalalignment='bottom', horizontalalignment='center', fontsize=10) + if dRxy is not None: axis.annotate(' r(A,B)=%.5g\n'%(dRxy), xy=(1.0,-1.1), xycoords='axes fraction', verticalalignment='top', horizontalalignment='center', fontsize=10) if len(xlabel+xunits)>0: plt.xlabel(label(xlabel, xunits)) if len(suptitle)>0: plt.suptitle(suptitle) @@ -458,7 +465,10 @@ def annotateStats(axis, sMin, sMax, sMean, sStd, sRMS): if npanels in [1, 3]: axis = plt.subplot(npanels,1,npanels) if dcolormap is None: dcolormap = chooseColorMap(dMin, dMax) - cmap, norm, extend = chooseColorLevels(dMin, dMax, dcolormap, clim=dlim, nbins=nbins, extend=dextend) + if dlim is None and dStd>0: + cmap, norm, dextend = chooseColorLevels(dMean-2.*dStd, dMean+2.*dStd, dcolormap, clim=dlim, nbins=nbins, extend='both', autocenter=True) + else: + cmap, norm, dextend = chooseColorLevels(dMin, dMax, dcolormap, clim=dlim, nbins=nbins, extend=dextend, autocenter=True) plt.pcolormesh(yCoord, zCoord, field1 - field2, cmap=cmap, norm=norm) if interactive: addStatusBar(yCoord, zCoord, field1 - field2) cb3 = plt.colorbar(fraction=.08, pad=0.02, extend=dextend) @@ -471,7 +481,7 @@ def annotateStats(axis, sMin, sMax, sMean, sStd, sRMS): annotateStats(axis, dMin, dMax, dMean, dStd, dRMS) if len(zlabel+zunits)>0: plt.ylabel(label(zlabel, zunits)) - axis.annotate(' r(A,B)=%.5g\n'%(dRxy), xy=(1.0,-0.20), xycoords='axes fraction', verticalalignment='bottom', horizontalalignment='center', fontsize=10) + axis.annotate(' r(A,B)=%.5g\n'%(dRxy), xy=(1.0,-1.07), xycoords='axes fraction', verticalalignment='top', horizontalalignment='center', fontsize=10) if len(ylabel+yunits)>0: plt.xlabel(label(ylabel, yunits)) if len(title3)>0: plt.title(title3) if len(suptitle)>0: plt.suptitle(suptitle) @@ -480,6 +490,89 @@ def annotateStats(axis, sMin, sMax, sMean, sStd, sRMS): if interactive: addInteractiveCallbacks() if show: plt.show(block=False) +def ztplot(field, t=None, z=None, + tlabel=None, tunits=None, zlabel=None, zunits=None, + splitscale=None, + title='', suptitle='', autocenter=False, + clim=None, colormap=None, extend=None, centerlabels=False, + nbins=None, landcolor=[.5,.5,.5], + aspect=[16,9], resolution=576, axis=None, + ignore=None, save=None, debug=False, show=False, interactive=False): + """ + Renders section plot of scalar field, field(x,z). + + Arguments: + field Scalar 2D array to be plotted. + t t (or time) coordinate (1D array). + z z coordinate (1D array). + tlabel The label for the t axis. Default 'Time'. + tunits The units for the t axis. Default 'Years'. + zlabel The label for the z axis. Default 'Elevation'. + zunits The units for the z axis. Default 'm'. + splitscale A list of depths to define equal regions of projection in the vertical, e.g. [0.,-1000,-6500] + title The title to place at the top of the panel. Default ''. + suptitle The super-title to place at the top of the figure. Default ''. + autocenter If clim generated by script, set to be centered on zero. Default False. + clim A tuple of (min,max) color range OR a list of contour levels. Default None. + colormap The name of the colormap to use. Default None. + extend Can be one of 'both', 'neither', 'max', 'min'. Default None. + centerlabels If True, will move the colorbar labels to the middle of the interval. Default False. + nbins The number of colors levels (used is clim is missing or only specifies the color range). + landcolor An rgb tuple to use for the color of land (no data). Default [.5,.5,.5]. + aspect The aspect ratio of the figure, given as a tuple (W,H). Default [16,9]. + resolution The vertical resolution of the figure given in pixels. Default 720. + axis The axis handle to plot to. Default None. + ignore A value to use as no-data (NaN). Default None. + save Name of file to save figure in. Default None. + debug If true, report stuff for debugging. Default False. + show If true, causes the figure to appear on screen. Used for testing. Default False. + interactive If true, adds interactive features such as zoom, close and cursor. Default False. + """ + + # Create coordinates if not provided + tlabel, tunits, zlabel, zunits = createTZlabels(t, z, tlabel, tunits, zlabel, zunits) + if debug: print('t,z label/units=',tlabel,tunits,zlabel,zunits) + if ignore is not None: maskedField = numpy.ma.masked_array(field, mask=[field==ignore]) + else: maskedField = field.copy() + field2 = maskedField.T + tCoord = t; zCoord = z + + # Diagnose statistics + sMin, sMax, sMean, sStd, sRMS = myStats(maskedField, None, debug=debug) + tLims = numpy.amin(tCoord), numpy.amax(tCoord) + zLims = numpy.amin(zCoord), numpy.amax(zCoord) + #zLims = boundaryStats(zCoord) + + # Choose colormap + if nbins is None and (clim is None or len(clim)==2): nbins=35 + if colormap is None: colormap = chooseColorMap(sMin, sMax) + cmap, norm, extend = chooseColorLevels(sMin, sMax, colormap, clim=clim, nbins=nbins, extend=extend, autocenter=autocenter) + + if axis is None: + setFigureSize(aspect, resolution, debug=debug) + #plt.gcf().subplots_adjust(left=.10, right=.99, wspace=0, bottom=.09, top=.9, hspace=0) + axis = plt.gca() + plt.pcolormesh(tCoord, zCoord, field2, cmap=cmap, norm=norm) + if interactive: addStatusBar(tCoord, zCoord, field2) + cb = plt.colorbar(fraction=.08, pad=0.02, extend=extend) + if centerlabels and len(clim)>2: cb.set_ticks( 0.5*(clim[:-1]+clim[1:]) ) + axis.set_axis_bgcolor(landcolor) + if splitscale is not None: + for zzz in splitscale[1:-1]: plt.axhline(zzz,color='k',linestyle='--') + axis.set_yscale('splitscale', zval=splitscale) + plt.xlim( tLims ); plt.ylim( zLims ) + axis.annotate('max=%.5g\nmin=%.5g'%(sMax,sMin), xy=(0.0,1.01), xycoords='axes fraction', verticalalignment='bottom', fontsize=10) + if sMean is not None: + axis.annotate('mean=%.5g\nrms=%.5g'%(sMean,sRMS), xy=(1.0,1.01), xycoords='axes fraction', verticalalignment='bottom', horizontalalignment='right', fontsize=10) + axis.annotate(' sd=%.5g\n'%(sStd), xy=(1.0,1.01), xycoords='axes fraction', verticalalignment='bottom', horizontalalignment='left', fontsize=10) + if len(tlabel+tunits)>0: plt.xlabel(label(tlabel, tunits)) + if len(zlabel+zunits)>0: plt.ylabel(label(zlabel, zunits)) + if len(title)>0: plt.title(title) + if len(suptitle)>0: plt.suptitle(suptitle) + + if save is not None: plt.savefig(save) + if interactive: addInteractiveCallbacks() + if show: plt.show(block=False) def chooseColorMap(sMin, sMax): """ @@ -491,13 +584,14 @@ def chooseColorMap(sMin, sMax): else: return 'dunneRainbow' -def chooseColorLevels(sMin, sMax, colorMapName, clim=None, nbins=None, steps=[1,2,2.5,5,10], extend=None, logscale=False): +def chooseColorLevels(sMin, sMax, colorMapName, clim=None, nbins=None, steps=[1,2,2.5,5,10], extend=None, logscale=False, autocenter=False): """ If nbins is a positive integer, choose sensible color levels with nbins colors. If clim is a 2-element tuple, create color levels within the clim range or if clim is a vector, use clim as contour levels. If clim provides more than 2 color interfaces, nbins must be absent. If clim is absent, the sMin,sMax are used as the color range bounds. + If autocenter is True and clim is None then the automatic color levels are centered. Returns cmap, norm and extend. """ @@ -506,7 +600,11 @@ def chooseColorLevels(sMin, sMax, colorMapName, clim=None, nbins=None, steps=[1, if len(clim)<2: raise Exception('clim must be at least 2 values long.') if nbins is None and len(clim)==2: raise Exception('nbins must be provided when clims specifies a color range.') if nbins is not None and len(clim)>2: raise Exception('nbins cannot be provided when clims specifies color levels.') - if clim is None: levels = MaxNLocator(nbins=nbins, steps=steps).tick_values(sMin, sMax) + if clim is None: + if autocenter: + levels = MaxNLocator(nbins=nbins, steps=steps).tick_values(min(sMin,-sMax), max(sMax,-sMin)) + else: + levels = MaxNLocator(nbins=nbins, steps=steps).tick_values(sMin, sMax) elif len(clim)==2: levels = MaxNLocator(nbins=nbins, steps=steps).tick_values(clim[0], clim[1]) else: levels = clim @@ -832,6 +930,24 @@ def createYZlabels(y, z, ylabel, yunits, zlabel, zunits): if zunits is None: zunits='m' return ylabel, yunits, zlabel, zunits +def createTZlabels(t, z, tlabel, tunits, zlabel, zunits): + """ + Checks that y and z labels are appropriate and tries to make some if they are not. + """ + if t is None: + if tlabel is None: tlabel='t' + if tunits is None: tunits='' + else: + if tlabel is None: tlabel='Time' + if tunits is None: tunits='' + if z is None: + if zlabel is None: zlabel='k' + if zunits is None: zunits='' + else: + if zlabel is None: zlabel='Elevation' + if zunits is None: zunits='m' + return tlabel, tunits, zlabel, zunits + def yzWeight(y, z): """ diff --git a/tools/analysis/meridional_overturning.py b/tools/analysis/meridional_overturning.py index cc069df0da..a9666bacf7 100755 --- a/tools/analysis/meridional_overturning.py +++ b/tools/analysis/meridional_overturning.py @@ -14,6 +14,7 @@ def run(): parser = argparse.ArgumentParser(description='''Script for plotting meridional overturning.''') parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'vh' and 'e'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-s','--suptitle', type=str, default='', help='''Super-title for experiment. Default is to read from netCDF file.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') @@ -101,13 +102,16 @@ def findExtrema(y, z, psi, min_lat=-90., max_lat=90., min_depth=0., mult=1.): m6plot.setFigureSize(npanels=1) cmap = plt.get_cmap('dunnePM') + if cmdLineArgs.suptitle != '': suptitle = cmdLineArgs.suptitle + ' ' + cmdLineArgs.label + else: suptitle = rootGroup.title + ' ' + cmdLineArgs.label + # Global MOC z = Zmod.min(axis=-1); psiPlot = MOCpsi(VHmod)*conversion_factor yy = y[1:,:].max(axis=-1)+0*z ci=m6plot.pmCI(0.,40.,5.) plotPsi(yy, z, psiPlot, ci, 'Global MOC [Sv]') plt.xlabel(r'Latitude [$\degree$N]') - plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + plt.suptitle(suptitle) findExtrema(yy, z, psiPlot, max_lat=-30.) findExtrema(yy, z, psiPlot, min_lat=25.) findExtrema(yy, z, psiPlot, min_depth=2000., mult=-1.) @@ -124,7 +128,7 @@ def findExtrema(y, z, psi, min_lat=-90., max_lat=90., min_depth=0., mult=1.): yy = y[1:,:].max(axis=-1)+0*z plotPsi(yy, z, psiPlot, ci, 'Atlantic MOC [Sv]') plt.xlabel(r'Latitude [$\degree$N]') - plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + plt.suptitle(suptitle) findExtrema(yy, z, psiPlot, min_lat=26.5, max_lat=27.) # RAPID findExtrema(yy, z, psiPlot, max_lat=-33.) findExtrema(yy, z, psiPlot) diff --git a/tools/analysis/poleward_heat_transport.py b/tools/analysis/poleward_heat_transport.py index 9bb97ddc38..1c82e0bad8 100755 --- a/tools/analysis/poleward_heat_transport.py +++ b/tools/analysis/poleward_heat_transport.py @@ -16,6 +16,7 @@ def run(): parser = argparse.ArgumentParser(description='''Script for plotting plotting poleward heat transport.''') parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'T_ady_2d' and 'T_diffy_2d'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-s','--suptitle', type=str, default='', help='''Super-title for experiment. Default is to read from netCDF file.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory or tarfile containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') @@ -134,6 +135,9 @@ def plotGandW(lat,trans,err): plt.plot([lat[n],lat[n]], [low[n],high[n]], 'c', linewidth=2.0) plt.scatter(lat,trans,marker='s',facecolor='cyan') + if cmdLineArgs.suptitle != '': suptitle = cmdLineArgs.suptitle + ' ' + cmdLineArgs.label + else: suptitle = rootGroup.title + ' ' + cmdLineArgs.label + # Global Heat Transport HTplot = heatTrans(advective,diffusive) yy = y[1:,:].max(axis=-1) @@ -142,7 +146,7 @@ def plotGandW(lat,trans,err): plt.plot(yobs,ECMWF['Global'],'k.',linewidth=0.5,label='ECMWF') plotGandW(GandW['Global']['lat'],GandW['Global']['trans'],GandW['Global']['err']) plt.xlabel(r'Latitude [$\degree$N]') - plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + plt.suptitle(suptitle) plt.legend(loc=0,fontsize=10) annotateObs() if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') @@ -162,7 +166,7 @@ def plotGandW(lat,trans,err): plt.plot(yobs,ECMWF['Atlantic'],'k.',linewidth=0.5,label='ECMWF') plotGandW(GandW['Atlantic']['lat'],GandW['Atlantic']['trans'],GandW['Atlantic']['err']) plt.xlabel(r'Latitude [$\degree$N]') - plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + plt.suptitle(suptitle) plt.legend(loc=0,fontsize=10) annotateObs() if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') @@ -184,7 +188,7 @@ def plotGandW(lat,trans,err): plt.xlabel(r'Latitude [$\degree$N]') annotateObs() if diffusive == None: annotatePlot('Warning: Diffusive component of transport is missing.') - plt.suptitle(rootGroup.title+' '+cmdLineArgs.label) + plt.suptitle(suptitle) plt.legend(loc=0,fontsize=10) if stream != None: plt.savefig(stream[2]) diff --git a/tools/analysis/zonal_S_bias_WOA05.py b/tools/analysis/zonal_S_bias_WOA05.py index 0ac7f7f804..30d32f1c0b 100755 --- a/tools/analysis/zonal_S_bias_WOA05.py +++ b/tools/analysis/zonal_S_bias_WOA05.py @@ -13,6 +13,7 @@ def run(): parser = argparse.ArgumentParser(description='''Script for plotting annual-average zonal salinity bias.''') parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'salt' and 'e'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-s','--suptitle', type=str, default='', help='''Super-title for experiment. Default is to read from netCDF file.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') @@ -67,6 +68,9 @@ def zonalAverage(S, eta, area, mask=1.): return numpy.sum( vols * S, axis=-1 ) / numpy.sum( vols, axis=-1 ), (mask*eta).min(axis=-1) ci=m6plot.pmCI(0.125,2.25,.25) + + if cmdLineArgs.suptitle != '': suptitle = cmdLineArgs.suptitle + ' ' + cmdLineArgs.label + else: suptitle = rootGroup.title + ' ' + cmdLineArgs.label # Global sPlot, z = zonalAverage(Smod, Zmod, area) @@ -74,13 +78,13 @@ def zonalAverage(S, eta, area, mask=1.): if stream != None: objOut = stream[0] else: objOut = cmdLineArgs.outdir+'/S_global_xave_bias_WOA05.png' m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Global zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', + suptitle=suptitle, title='''Global zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=objOut) if stream == None: m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1='Global zonal-average salinity [ppt]', title2='''WOA'05 salinity [ppt]''', clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', @@ -94,13 +98,13 @@ def zonalAverage(S, eta, area, mask=1.): if stream != None: objOut = stream[1] else: objOut = cmdLineArgs.outdir+'/S_Atlantic_xave_bias_WOA05.png' m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Atlantic zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', + suptitle=suptitle, title='''Atlantic zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=objOut) if stream == None: m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1='Atlantic zonal-average salinity [ppt]', title2='''WOA'05 salinity [ppt]''', clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', @@ -114,13 +118,13 @@ def zonalAverage(S, eta, area, mask=1.): if stream != None: objOut = stream[2] else: objOut = cmdLineArgs.outdir+'/S_Pacific_xave_bias_WOA05.png' m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Pacific zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', + suptitle=suptitle, title='''Pacific zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=objOut) if stream == None: m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1='Pacific zonal-average salinity [ppt]', title2='''WOA'05 salinity [ppt]''', clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', @@ -134,13 +138,13 @@ def zonalAverage(S, eta, area, mask=1.): if stream != None: objOut = stream[3] else: objOut = cmdLineArgs.outdir+'/S_Indian_xave_bias_WOA05.png' m6plot.yzplot( sPlot - sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title='''Indian zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', + suptitle=suptitle, title='''Indian zonal-average salinity bias (w.r.t. WOA'05) [ppt]''', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=objOut) if stream == None: m6plot.yzcompare( sPlot, sObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1='Indian zonal-average salinity [ppt]', title2='''WOA'05 salinity [ppt]''', clim=m6plot.linCI(20,30,10, 31,39,.5), colormap='dunneRainbow', extend='both', diff --git a/tools/analysis/zonal_T_bias_WOA05.py b/tools/analysis/zonal_T_bias_WOA05.py index 66c36540ea..f354e206fb 100755 --- a/tools/analysis/zonal_T_bias_WOA05.py +++ b/tools/analysis/zonal_T_bias_WOA05.py @@ -13,6 +13,7 @@ def run(): parser = argparse.ArgumentParser(description='''Script for plotting annual-average zonal temperature bias.''') parser.add_argument('annual_file', type=str, help='''Annually-averaged file containing 3D 'temp' and 'e'.''') parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') + parser.add_argument('-s','--suptitle', type=str, default='', help='''Super-title for experiment. Default is to read from netCDF file.''') parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') parser.add_argument('-g','--gridspec', type=str, required=True, help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') @@ -70,20 +71,23 @@ def zonalAverage(T, eta, area, mask=1.): return numpy.sum( vols * T, axis=-1 ) / numpy.sum( vols, axis=-1 ), (mask*eta).min(axis=-1) ci=m6plot.pmCI(0.25,4.5,.5) - + + if cmdLineArgs.suptitle != '': suptitle = cmdLineArgs.suptitle + ' ' + cmdLineArgs.label + else: suptitle = rootGroup.title + ' ' + cmdLineArgs.label + # Global tPlot, z = zonalAverage(Tmod, Zmod, area) tObsPlot, _ = zonalAverage(Tobs, Zobs, area) if stream != None: objOut = stream[0] else: objOut = cmdLineArgs.outdir+'/T_global_xave_bias_WOA05.png' m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Global zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', + suptitle=suptitle, title=r'''Global zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=objOut) if stream == None: m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1=r'Global zonal-average $\theta$ [$\degree$C]', title2=r'''WOA'05 $\theta$ [$\degree$C]''', clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', @@ -97,13 +101,13 @@ def zonalAverage(T, eta, area, mask=1.): if stream != None: objOut = stream[1] else: objOut = cmdLineArgs.outdir+'/T_Atlantic_xave_bias_WOA05.png' m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Atlantic zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', + suptitle=suptitle, title=r'''Atlantic zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=objOut) if stream == None: m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1=r'Atlantic zonal-average $\theta$ [$\degree$C]', title2=r'''WOA'05 $\theta$ [$\degree$C]''', clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', @@ -117,13 +121,13 @@ def zonalAverage(T, eta, area, mask=1.): if stream != None: objOut = stream[2] else: objOut = cmdLineArgs.outdir+'/T_Pacific_xave_bias_WOA05.png' m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Pacific zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', + suptitle=suptitle, title=r'''Pacific zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=objOut) if stream == None: m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1=r'Pacific zonal-average $\theta$ [$\degree$C]', title2=r'''WOA'05 $\theta$ [$\degree$C]''', clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', @@ -137,13 +141,13 @@ def zonalAverage(T, eta, area, mask=1.): if stream != None: objOut = stream[3] else: objOut = cmdLineArgs.outdir+'/T_Indian_xave_bias_WOA05.png' m6plot.yzplot( tPlot - tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, title=r'''Indian zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', + suptitle=suptitle, title=r'''Indian zonal-average $\theta$ bias (w.r.t. WOA'05) [$\degree$C]''', clim=ci, colormap='dunnePM', centerlabels=True, extend='both', save=objOut) if stream == None: m6plot.yzcompare( tPlot, tObsPlot , y, z, splitscale=[0., -1000., -6500.], - suptitle=rootGroup.title+' '+cmdLineArgs.label, + suptitle=suptitle, title1=r'Indian zonal-average $\theta$ [$\degree$C]', title2=r'''WOA'05 $\theta$ [$\degree$C]''', clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max',