File contents

import cdms, cdutil, MA, vcs, cdtime
import string, Numeric, time, MV, sys, os
from regrid import Regridder
#
# From the Climate Research Unit (CRU), University of East Anglia
# available at http://www.cru.uea.ac.uk/cru/data/
# These data represent anomalies from a merged 2-meter surface
# air temperature over land and ssts over the ocean on a
# 5-degree grid.
# Received 03/2004.
#
# Extract ERA40 data over longest matching time period, regrid both
# 2-meter temperature and ssts to a 5 degree grid then apply a land/sea
# mask to separate appropriate data for a merged product.  Add these two
# fields, apply the final "spatial missing mask" and add applicable 
# metadata for ease of use.

# In this example we will be working with the sample data and much
# shorter time spans then you would want to do.
#
# Examples of:
#
# getting slices of data
# getting metadata
# regridding
# simple plotting
# manipulating and applying a land/sea mask
# calculating an annual cycle and departures from it
# the application of a "spatial missing mask"
# creating a new variable
# writing data out to a netcdf file
#
#
file1 = os.path.join(sys.prefix, 'sample_data/hadcrut2_sample.nc')
a = cdms.open(file1)
#
print a.listvariable()

start_time = a.getAxis('time').asComponentTime()[0]
end_time = a.getAxis('time').asComponentTime()[-1]

#
data=a('temanom',time=(start_time, end_time),latitude=(-90,90))
print data.shape
##  (36, 36, 72)

x=vcs.init()
x.setcolormap('default')
x.plot(data)
print "Just plotted CRU data - annual cycle temperature anomaly"
print "Press the Return key to see next plot."
sys.stdin.readline()


# get grid for regridding
grid1=data.getGrid()
# see how it looks like
grid1.info
## <bound method TransientRectGrid.info of Grid has Python id -0x491e3394.
## Gridtype: generic
## Grid shape: (36, 72)
## Order: yx
## >

# get "spatial missing mask"
mask1=data.mask()


# get metadata for final desired data (e.g., latitudes, longitudes, time)
lat=data.getLatitude()
lon=data.getLongitude()
tim=data.getTime()
# close the file
a.close()
#
# Get the ERA40 data for both 2-meter temperature (tas) and ssts (sst) on
# original grid
#
file2 = os.path.join(sys.prefix, 'sample_data/era40_tas_sample.nc')
b = cdms.open(file2)
# get data
print "reading 'tas' data, please wait..."
tas=b('tas')
print tas.shape
# close the file
b.close()

file3 = os.path.join(sys.prefix, 'sample_data/era40_sst_sample.nc')
b = cdms.open(file3)
# get data
print "reading 'sst' data, please wait..."
sst=b('sst')
print sst.shape
# get grid for regridding
grid2=sst.getGrid()
# close the file
b.close()

# setup a regridding function (as: fromgrid, togrid)
#
regridfunc=Regridder(grid2,grid1)
# put in data at the "fromgrid" resolution
sst_new=regridfunc(sst)
tas_new=regridfunc(tas)
# tas_new and sst_new are the data on the 5 degree grid
#

# example plot of original ssts (approx 1.1 degree grid)

print "Just plotted tas original 2m resolution data"
print "Press the Return key to see next plot."
sys.stdin.readline()

# example plot of original ssts (approx 1.1 degree grid)
x.clear()
x.plot(sst)
print "Just plotted sst original 2m resolution data"
print "Press the Return key to see next plot."
sys.stdin.readline()


#
# example plot of regridded ssts (5 degree grid)
y=vcs.init()
y.setcolormap('default')
y.plot(sst_new)
print "Just plotted regridded sst data"
print "Press the Return key to see next plot."
sys.stdin.readline()

#
# extract a land/sea mask and regrid it to our desired 5 degree grid
# (these data are percent land coverage [0-100])
#
file4 = os.path.join(sys.prefix, 'sample_data/geo.1deg.ctl')
c = cdms.open(file4)
fraction=c('sftlf',squeeze=1)
c.close()

# get data grid for regrdding
grid3=fraction.getGrid()
# setup regrid function (as above)
regridfunc=Regridder(grid3,grid1)
# regrid mask values
fraction=regridfunc(fraction)
#
# I selected 50% or more coverage in a box is land and less than or
# equal to 50% coverage is ocean.  All other values in the arrays are zeros
#
land=Numeric.where(Numeric.greater(fraction.filled(),50.),1.,0.)
ocean=Numeric.where(Numeric.less_equal(fraction.filled(),50.),1.,0)

#
masked_sst=Numeric.multiply(sst_new.filled(),ocean)
masked_tas=Numeric.multiply(tas_new.filled(),land)

x.clear()
x.plot(masked_sst)
print "Just plotted masked sst data"
print "Press the Return key to see next plot."
sys.stdin.readline()

y.clear()
y.plot(masked_tas)
print "Just plotted masked tas data"
print "Press the Return key to see next plot."
sys.stdin.readline()

# add land and ocean contributions for the merged product
merged=masked_sst+masked_tas
# add metadata to this numeric array
merged=cdms.createVariable(merged,axes=(tim,lat,lon),typecode='f',id='merged_tas_sst')
merged.id='merged_tas_sst'
merged.set_fill_value(1e20)
cdutil.setTimeBoundsMonthly(merged)

x.clear()
x.plot(merged)
print "Just plotted merged tas and sst masked data"
print "Press the Return key to see next plot."
sys.stdin.readline()


print "writing merged data to output 'era40_merged_tas_sst.nc' NetCDF file"
# write out the total temperature data to a netcdf file
o=cdms.open('era40_merged_tas_sst.nc','w')
o.write(merged)
#
# To match the CRU data we must calculate and subtract an annual cycle
# (from the base period 1991-1993, inclusive)
#
#start_time = cdtime.comptime(1991,1,1)
#end_time   = cdtime.comptime(1993,12,1)
#
# define the annualcycle
#
print "calculating annual cycle anomaly..."
ac=cdutil.ANNUALCYCLE.climatology(merged(time=(start_time, end_time, 'co')))
# use the defined annual cycle and generate anomalies
merged_an=cdutil.ANNUALCYCLE.departures(merged,ref=ac)
#
# add metadata to the new anomaly variable
#
merged_an=cdms.createVariable(merged_an,axes=(tim,lat,lon),typecode='f',id='anomalies_merged_tas_sst')
merged_an.id='anomalies_merged_tas_sst'
#
# Lastly apply the "spatial missing mask" to these data
#
merged_an=MV.masked_where(MV.equal(mask1,1),merged_an)

print "writing merged anomaly data to output 'era40_merged_tas_sst.nc' NetCDF file"
# write out data
o.write(merged_an)
o.close()


x.clear()
x.plot(merged_an)
print ""
print 'Just plotted merged tas+sst AERA40 data annual cycle anomaly'
print "Press the Return key to see next plot."
sys.stdin.readline()


y.clear()
y.plot(data)
print 'Just plotted CRU data annual cycle anomaly for comparison'
print ''
print "---------------"
print "Press the Return key to EXIT."
sys.stdin.readline()

#
# done!
#