From 9f0b90f1150b1642889364844501783ea38befda Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Fri, 20 Jun 2014 17:33:08 -0400 Subject: [PATCH] Adding a new example analysis script. - EddyKineticEnergy.py calculates the surface eddy kinetic energy annual mean - Required adding a logarithmic color plot option to 6plot.xyplot - Serves as an example analysis script that can be called inside the refineDiag script to be done annually - Motivated by the plot in http://www.gfdl.noaa.gov/cms-filesystem-action/user_files/kd/pdf/c37_dixon_th85b.pdf - Unaware of the scientific value, just an example. --- tools/analysis/EddyKineticEnergy.py | 56 +++++++++++++++++++++++++++++ tools/analysis/m6plot.py | 7 ++-- 2 files changed, 61 insertions(+), 2 deletions(-) create mode 100755 tools/analysis/EddyKineticEnergy.py diff --git a/tools/analysis/EddyKineticEnergy.py b/tools/analysis/EddyKineticEnergy.py new file mode 100755 index 0000000000..decf07d0cc --- /dev/null +++ b/tools/analysis/EddyKineticEnergy.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python + +import netCDF4 +import matplotlib.pyplot as plt +import numpy +import m6plot +import math + +try: import argparse +except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.') + +parser = argparse.ArgumentParser(description='''Script for plotting annual-average SST bias.''') +parser.add_argument('annual_file', type=str, help='''Annually-averaged file contained 3D temp.''') +parser.add_argument('-l','--label', type=str, default='', help='''Label to add to the plot.''') +parser.add_argument('-o','--outdir', type=str, default='.', help='''Directory in which to place plots.''') +parser.add_argument('-g','--gridspec', type=str, + default='/archive/gold/datasets/MOM6z_SIS_025/siena/mosaic.unpacked', + help='''Directory containing mosaic/grid-spec files (ocean_hgrid.nc and ocean_mask.nc).''') +cmdLineArgs = parser.parse_args() + +rootGroup = netCDF4.Dataset( cmdLineArgs.annual_file ) +if 'ssu' not in rootGroup.variables: raise Exception('Could not find "ssu" in file "%s"'%(cmdLineArgs.annual_file)) +if 'ssv' not in rootGroup.variables: raise Exception('Could not find "ssv" in file "%s"'%(cmdLineArgs.annual_file)) + +x = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['x'][::2,::2] +y = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['y'][::2,::2] +msk = netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_mask.nc').variables['mask'][:] +area = msk*netCDF4.Dataset(cmdLineArgs.gridspec+'/ocean_hgrid.nc').variables['area'][:,:].reshape([msk.shape[0], 2, msk.shape[1], 2]).sum(axis=-3).sum(axis=-1) +msk = numpy.ma.array(msk, mask=(msk==0)) + +#Tobs = netCDF4.Dataset( '/archive/gold/datasets/MOM6z_SIS_025/siena/obs/WOA05_ptemp_salt_annual_v1.nc' ).variables['temp'][1] + +ssu = rootGroup.variables['ssu'] +ssu_mean = ssu[:].mean(axis=0) +ssv = rootGroup.variables['ssv'] +ssv_mean = ssv[:].mean(axis=0) + +eke = ((ssu[:] - ssu_mean)**2 + (ssv[:] - ssv_mean)**2)/2 + +eke_mean = 10000 * eke[:].mean(axis=0) +#eke_mean = numpy.log10(eke_mean[:]) + +ci=m6plot.pmCI(0.0,0.5,0.1) + +m6plot.xyplot( eke_mean , x, y, area=area, + suptitle=rootGroup.title+' '+cmdLineArgs.label, title='Eddy Kinetic Energy annual mean [(cm/s)^2]', + clim=ci, logscale=True, + save=cmdLineArgs.outdir+'/EKE_mean.png') + +#m6plot.xycompare( temp, Tobs , x, y, area=area, +# suptitle=rootGroup.title+' '+cmdLineArgs.label, +# title1='SST [$\degree$C]', +# title2='WOA\'05 SST [$\degree$C]', +# clim=m6plot.linCI(-2,29,.5), colormap='dunneRainbow', extend='max', +# dlim=ci, dcolormap='dunnePM', dextend='both', centerDCB=True, +# save=cmdLineArgs.outdir+'/SST_bias_WOA05.3_panel.png') diff --git a/tools/analysis/m6plot.py b/tools/analysis/m6plot.py index aea2c8a0e2..63855c2b85 100644 --- a/tools/analysis/m6plot.py +++ b/tools/analysis/m6plot.py @@ -4,7 +4,7 @@ import numpy, numpy.matlib import matplotlib.pyplot as plt -from matplotlib.colors import BoundaryNorm, ListedColormap +from matplotlib.colors import BoundaryNorm, ListedColormap, LogNorm from matplotlib.ticker import MaxNLocator import math import m6toolbox @@ -15,7 +15,7 @@ def xyplot(field, x=None, y=None, area=None, clim=None, colormap=None, extend=None, centerlabels=False, nbins=None, landcolor=[.5,.5,.5], aspect=[16,9], resolution=576, - ignore=None, save=None, debug=False, show=False, interactive=False): + ignore=None, save=None, debug=False, show=False, interactive=False, logscale=False): """ Renders plot of scalar field, field(x,y). @@ -46,6 +46,7 @@ def xyplot(field, x=None, y=None, area=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. + logscale If true, use logaritmic coloring scheme. Default False. """ # Create coordinates if not provided @@ -65,6 +66,8 @@ def xyplot(field, x=None, y=None, area=None, if colormap==None: colormap = chooseColorMap(sMin, sMax) cmap, norm, extend = chooseColorLevels(sMin, sMax, colormap, clim=clim, nbins=nbins, extend=extend) + if logscale: norm = LogNorm() + if axis==None: setFigureSize(aspect, resolution, debug=debug) #plt.gcf().subplots_adjust(left=.08, right=.99, wspace=0, bottom=.09, top=.9, hspace=0)