Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding a new example analysis script. #26

Merged
merged 2 commits into from
Jul 22, 2014
Merged

Adding a new example analysis script. #26

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9f0b90f
Adding a new example analysis script.
nikizadehgfdl Jun 20, 2014
70d47e6
Merge branch 'dev/master' into user/nnz/EddyKineticEnergyAnalysisScript
nikizadehgfdl Jul 22, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 56 additions & 0 deletions tools/analysis/EddyKineticEnergy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There should be an i-average of u^2 from u-points to t-points and a j-average of v^2 from v-pointsto t-points. Use numpy.roll().


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')
7 changes: 5 additions & 2 deletions tools/analysis/m6plot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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).

Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down