#--step 1 pick up the data for 1980-2015 from infile# seldate : give range of datescdo.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 filecdo.showyear(input="output1.nc")

#--output can be read directory as a netcdf fileimportnumpyasnpf1=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)"

#-- 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[:])

np.shape(anom)= (432, 1, 100, 190)

In [20]:

#-- Draw anomalies for 1997 and 1998%matplotlib inline
importosimportpandasaspdfromnetCDF4importnum2dateimportnumpyasnpimportmatplotlib.pyplotaspltimportmatplotlib.cmascmfrommpl_toolkits.basemapimportBasemapanom=f2.variables["temp"]#read anomaly from previous resutstime=f2.variables["time"]#read timedlons=f2.variables["lon"][:]# read longitudesdlats=f2.variables["lat"][:]# read latitudesndate=num2date(time[:],units=time.units,calendar=time.calendar)#convert time to datedates=pd.DatetimeIndex(ndate)# convert date to date object in pandassyear=1997# start of ploteyear=1998# end of plotfig=plt.figure(figsize=(24,24))# set figure environemntcontours=np.arange(-4.0,4.0,0.5)# set contourscmap=cm.bwr# set colormapsmon={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 dictionaryii=1foryearinrange(syear,eyear+1):#loop for each yearidxs=np.where(dates.year==year)[0]forim,idxinenumerate(idxs):#loop for each month in a yearax=plt.subplot(6,4,ii)# set up panel plotm=Basemap(projection='mill',llcrnrlat=dlats[0],urcrnrlat=dlats[-1],llcrnrlon=dlons[0],urcrnrlon=dlons[-1],lat_ts=5,resolution='c')# draw basemapm.drawcoastlines()# draw coast lineflon,flat=np.meshgrid(dlons,dlats)# meshgridx,y=m(flon,flat)# lon,lat => x,ycs=m.contourf(x,y,anom[idx,0,:,:],contours,cmap=cmap,extend="both")# plot contoursax.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 labelii=ii+1# set up labelcax=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 layoutplt.tight_layout(pad=0.2,w_pad=0.2,h_pad=0.3)# showplt.show()