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()