Python CDO (Climate Data Operations)
In [1]:
from cdo import *
#--Initialize CDO
cdo=Cdo()
#print cdo.operators # query operators
3. If the avove initialization does not work, you need to edit the function “getOperators” in cdo.py¶
- Modify the function “getOperators” as follows.
$vi (YOUR_ANACONDA_DIR)/lib/python2.7/site-packages/cdo.pydef getOperators(self): import os proc = subprocess.Popen([self.CDO,'--operators'],stderr = subprocess.PIPE,stdout = subprocess.PIPE) ret = proc.communicate() ops = ret[1].decode("utf-8")[0:-1].split(os.linesep)[1:-1] opes=[] for text in ret[0].decode("utf-8")[0:-1].split("\n"): s = re.sub("\s+" , " ", text) opes.append(s.split(" ")[0]) return opes
4. Basic syntax¶
cdo.OPERATOR("PARAETERS", input="FILENAME.nc", output="FILENAME.nc")
Parameters and output are not always necessary depending on operators.
In [11]:
#--input file (Monthly HadISST data 1870.01-2016.02)
infile="/work/hnm/dataset/HadISST1_1/netcdf_data_temp.20160412/HadISST1_SST_187001-201602.nc"
#infile="HadISST1_SST_187001-201602.nc"
# sinfo : Display file info
cdo.sinfo(input=infile)
Out[11]:
In [12]:
#--step 1 pick up the data for 1980-2015 from infile
# seldate : give range of dates
cdo.seldate('1980-01-01,2015-12-31', input=infile, output="output1.nc")
# selyear : give range of years (same )
cdo.selyear('1980/2015', input=infile, output="output1.nc")
# showyear : display years in the file
cdo.showyear(input="output1.nc")
Out[12]:
In [13]:
#--step 2: extract domain of interest
# sellonlatbox
cdo.sellonlatbox('100,280,-50,50', input="output1.nc", output="output2.nc")
cdo.griddes(input="output2.nc") # show grid info
Out[13]:
In [14]:
#--step 3 make multiple-year monthly average
# ymonavg input output (output: (12,jmax,imax))
cdo.ymonavg(input="output2.nc", output="cmean.nc")
cdo.info(input="cmean.nc") #show data info
Out[14]:
Piping¶
Before passing the data to the main operator, some operations can be done beforehand.
cdo.OPERATOR3("PARAM3",input="-OPERATOR2,PARAM2 -OPERATOR1,PARAM1 INPUT.nc", output="OUTPUT.nc")
In [15]:
#--steps 1-3 can be combined using "piping".
cdo.ymonavg(input="-sellonlatbox,100,280,-50,50 -seldate,1980-01-01,2015-12-31 %s" % (infile), output="cmean.nc")
cdo.info(input="cmean.nc") #show data info
Out[15]:
Reading output as a netcdf file¶
In [16]:
#--output can be read directory as a netcdf file
import numpy as np
f1 = cdo.ymonavg(input="-sellonlatbox,100,280,-50,50 -seldate,1980-01-01,2015-12-31 %s" % (infile),
options='-f nc', returnCdf=True) # As if f1=Dataset("cmean.nc","r")
clim = f1.variables["temp"][:]
print "np.shape(clim)=",np.shape(clim),"=(12-months, zmax, jmax, imax)"
In [17]:
#--step 4 make anomaly (the raw data is subtracted from the climatological monthly mean)
# ymonsub raw-data mean-data anomaly-data
cdo.ymonsub(input="output2.nc cmean.nc", output="anom.nc")
Out[17]:
In [18]:
#-- Do step1-step4 with one line without temporaly files using "piping".
f2 = cdo.ymonsub(input='-sellonlatbox,100,290,-50,50 -seldate,1980-01-01,2015-12-31 %s \
-ymonavg -sellonlatbox,100,290,-50,50 -seldate,1980-01-01,2015-12-31 %s' % (infile,infile),
options='-f nc', returnCdf=True)
anom = f2.variables["temp"]
print "np.shape(anom)=",np.shape(anom[:])
In [20]:
#-- Draw anomalies for 1997 and 1998
%matplotlib inline
import os
import pandas as pd
from netCDF4 import num2date
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
anom = f2.variables["temp"] #read anomaly from previous resuts
time = f2.variables["time"] #read time
dlons = f2.variables["lon"][:] # read longitudes
dlats = f2.variables["lat"][:] # read latitudes
ndate = num2date(time[:],units=time.units,calendar=time.calendar) #convert time to date
dates = pd.DatetimeIndex(ndate) # convert date to date object in pandas
syear=1997 # start of plot
eyear=1998 # end of plot
fig = plt.figure(figsize=(24,24)) # set figure environemnt
contours=np.arange(-4.0,4.0,0.5) # set contours
cmap=cm.bwr # set colormap
smon={1:"JAN",2:"FEB",3:"MAR",4:"APR",5:"MAY",6:"JUN",7:"JUL",8:"AUG",9:"SEP",10:"OCT",11:"NOV",12:"DEC"} # set dictionary
ii=1
for year in range(syear,eyear+1): #loop for each year
idxs = np.where(dates.year == year)[0]
for im, idx in enumerate(idxs): #loop for each month in a year
ax = plt.subplot(6,4,ii) # set up panel plot
m = Basemap(projection='mill',llcrnrlat=dlats[0], urcrnrlat=dlats[-1], llcrnrlon=dlons[0],urcrnrlon=dlons[-1],
lat_ts=5, resolution='c') # draw basemap
m.drawcoastlines() # draw coast line
flon,flat=np.meshgrid(dlons,dlats) # meshgrid
x,y=m(flon,flat) # lon,lat => x,y
cs = m.contourf(x,y,anom[idx,0,:,:],contours, cmap=cmap, extend="both") # plot contours
ax.text(0.015,0.9, "%4.4i %s" % (year, smon[im+1]), transform=ax.transAxes, fontsize=20, bbox=dict(boxstyle='square',
fc='w', alpha=1.0), zorder=100) # draw label
ii=ii+1
# set up label
cax = fig.add_axes([0.2, -0.050, 0.6, 0.03])
art = plt.colorbar(cs, cax, orientation='horizontal')
art.set_label('SST Anomaly [K]', fontsize=20)
art.ax.tick_params(labelsize=18)
# adjust layout
plt.tight_layout(pad=0.2, w_pad=0.2, h_pad=0.3)
# show
plt.show()
2. EOF using CDO¶
In [22]:
#cdo.eof("3", input="anom.nc", output="eval.nc eof.nc") #error because anom.nc contains missing value
cdo.eof("3", input="-setmisstoc,0 anom.nc", output="eval.nc eof.nc") # compute first 3 modes
cdo.eofcoeff(input="eof.nc anom.nc", output="eof_pc_") # compute principal coefficients
# eof_pc_0000{0,1,2,3}.nc are created
Out[22]:
In [23]:
#--read output data
ff = cdo.seltimestep("1/3", input="eof.nc", options='-f nc', returnCdf=True) # read eigen vectors for the first 3 modes
eof = ff.variables["temp"][:,0,:,:] # read eigen vectors
print "np.shape(eof)=",np.shape(eof),"=(mode,jmax,imax)"
In [24]:
#--plot
fig = plt.figure(figsize=(12,16))
contours=np.arange(-4.0,4.0,0.5) # create contours
cmap=cm.bwr #set color bar
dlons = ff.variables["lon"][:] # get longitudes
dlats = ff.variables["lat"][:] # get latitudes
sign=[1,1,-1] # sign for [1st,2nd,3rd mode]
jj=0
for ii in range(3): # loop for first 3 modes
#----draw eigen vector
jj=jj+1
ax = plt.subplot(5,2,jj)
m = Basemap(projection='mill',llcrnrlat=dlats[0], urcrnrlat=dlats[-1], llcrnrlon=dlons[0],urcrnrlon=dlons[-1],
lat_ts=5, resolution='c')
m.drawcoastlines()
flon,flat=np.meshgrid(dlons,dlats)
x,y=m(flon,flat)
cs = m.contourf(x,y,eof[ii,:,:]*sign[ii],contours, cmap=cmap, extend="both")
ax.text(0.015,0.9, "EOF%i" % (ii+1), transform=ax.transAxes, fontsize=20, bbox=dict(boxstyle='square',
fc='w', alpha=1.0), zorder=100)
#----plot PC
jj=jj+1
ax = plt.subplot(5,2,jj)
f3 = cdo.selname("temp",input="eof_pc_%5.5i.nc" % (ii), options='-f nc', returnCdf=True)
pc = f3.variables["temp"][:,0,0,0]*sign[ii]
pc = (pc - pc.mean())/pc.std() #normalize
xx = np.arange(0,(2015-1980+1)*12,1)
lns3 = ax.plot(xx, pc, "-", ms=10, color="black", linewidth=2)
ax.set_ylabel('PC')
ax.set_ylim((-3.0,3.0))
ax.axhline(y=0,ls='--',linewidth=1, color='k')
xdata=np.arange(0,(2015-1980+1)*12,12*5)
xlabels=np.arange(1980,2016,5)
plt.xticks( xdata, xlabels, rotation=0)
cax = fig.add_axes([0.2, 0.35, 0.6, 0.03])
art = plt.colorbar(cs, cax, orientation='horizontal')
art.set_label('SST Anomaly [K]', fontsize=20)
art.ax.tick_params(labelsize=18)
plt.tight_layout(pad=0.2, w_pad=0.2, h_pad=0.3)
plt.show()
3. Other Userful Operators¶
In [25]:
# Add -270.15 (convert C to K)
cdo.addc("273.15",input="cmean.nc", output="cmean_in_K.nc")
Out[25]:
In [26]:
# Zonal mean
cdo.zonmean(input="cmean.nc", output="cmean_zmean.nc")
Out[26]:
In [27]:
# Bilinear interpolation to 144x73 (2.5deg x 2.5deg) grid
cdo.remapbil("r144x73",input="cmean.nc", output="cmean_144x73.nc")
# To T42 Gaussian Grids
cdo.remapbil("t42grid",input="cmean.nc", output="cmean_T42.nc")
Out[27]:
In [28]:
# Bandpass filter (ex. 2-10 year)
cdo.bandpass("%f,%f" % (12./120.,12./24.), input="-setmisstoc,0 output1.nc", output="out_2-10yr.nc")
Out[28]:
In [29]:
# Lowpass filter (ex. 10 year)
cdo.lowpass("%f" % (12./120.), input="-setmisstoc,0 output1.nc", output="out_low10yr.nc")
Out[29]:
In [30]:
# Detrend
cdo.detrend(input="output1.nc",output="detrend.nc")
Out[30]:
In [32]:
#--plot for Nino-3.4 mean SST
# fldmean : area weighted
raw = cdo.fldmean(input="-sellonlatbox,170,240,-5,5 output1.nc", options='-f nc', returnCdf=True)
low10 = cdo.fldmean(input="-sellonlatbox,170,240,-5,5 out_low10yr.nc", options='-f nc', returnCdf=True)
b2_10 = cdo.fldmean(input="-sellonlatbox,170,240,-5,5 out_2-10yr.nc", options='-f nc', returnCdf=True)
#250,270,0,10
# plot timeseries
fig=plt.figure(figsize=(12.0,6.0), dpi=100, facecolor='w', edgecolor='k')
ax = plt.subplot(1,1,1)
xx = np.arange(0,(2015-1980+1)*12,1)
plt.plot(xx, raw.variables["temp"][:,0,0,0], "-", color="black", lw=3, label="raw data")
plt.plot(xx, low10.variables["temp"][:,0,0,0], "-", color="red", lw=3, label="10-yr lowpass")
plt.plot(xx, b2_10.variables["temp"][:,0,0,0]+raw.variables["temp"][:,0,0,0].mean(axis=0),
"-", color="blue", lw=3, label="2-10-yr bandpass")
xx = np.arange(0,(2015-1980+1)*12,1)
xdata=np.arange(0,(2015-1980+1)*12,12*5)
xlabels=np.arange(1980,2015+1,5)
plt.xticks( xdata, xlabels, rotation=0)
ax.axhline(y=0,ls='--',linewidth=1, color='k')
for x in xdata:
ax.axvline(x=x,ls='--',linewidth=1,color='k')
ax.set_ylim((25,30))
plt.legend(loc=2,ncol=4)
plt.show()
In [33]:
# Compute Nino-3.4 Index
infile="/work/hnm/dataset/HadISST1_1/netcdf_data_temp.20160412/HadISST1_SST_187001-201602.nc"
#infile="HadISST1_SST_187001-201602.nc"
# Nino 3.4 index
f34 = cdo.ymonsub(input='-fldmean -sellonlatbox,170,240,-5,5 -selyear,1980/2015 %s \
-ymonavg -fldmean -sellonlatbox,170,240,-5,5 -selyear,1980/2015 %s' % (infile,infile),
options='-f nc', returnCdf=True)
nino34 = f34.variables["temp"][:,0,0,0]
# plot timeseries
fig=plt.figure(figsize=(12.0,6.0), dpi=100, facecolor='w', edgecolor='k')
ax = plt.subplot(1,1,1)
xx = np.arange(0,(2015-1980+1)*12,1)
xdata=np.arange(0,(2015-1980+1)*12,12*5)
xlabels=np.arange(1980,2015+1,5)
ax.plot(xx, nino34, "-", ms=10, color="black", linewidth=4, label="Nino-3.4")
ax.set_ylabel('Nino-3.4 Index', fontsize=14)
ax.set_ylim((-2.5,2.5))
ax.axhline(y=0,ls='--',linewidth=1, color='k')
for x in xdata:
ax.axvline(x=x,ls='--',linewidth=1,color='k')
xdata=np.arange(0,(2015-1980+1)*12,12*5)
xlabels=np.arange(1980,2016,5)
xticks = plt.xticks(xdata, xlabels, rotation=0)
In [ ]: