MV-EOF的作图程序分享（NCL）+两个小问题

我是小猪 · 发表于 2017-4-7 00:59:06

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
; ==============================================================
begin
latS = 0.
latN = 50.
lonL = 100.
lonR = 140.
neof=2
yrStrt = 1979
yrLast = 2006
ymStrt = yrStrt*100 + 1
ymLast = yrLast*100 + 12
f =addfile("/disk3/gpcp/precip.mon.mean.nc","r")
fu=addfile("/disk3/ncep2/monthly/pressure/uwnd.mon.mean.nc","r")
fv=addfile("/disk3/ncep2/monthly/pressure/vwnd.mon.mean.nc","r")
TIME = f->time
YYYY = cd_calendar(TIME,-1)/100
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
pre=f->precip(iYYYY,{latS:latN},{lonL:lonR})
pre=dim_standardize_n_Wrap(pre, 0, 0)
precip=month_to_season(pre, "JJA")
delete(TIME)
delete(YYYY)
delete(iYYYY)
TIME = fu->time
YYYY = cd_calendar(TIME,-1)/100
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
uwnd=short2flt(fu->uwnd(iYYYY,{850},{latS:latN},{lonL:lonR}))
uwnd=dim_standardize_n_Wrap(uwnd, 0, 0)
u =month_to_season(uwnd, "JJA")
delete(TIME)
delete(YYYY)
delete(iYYYY)
TIME = fv->time
YYYY = cd_calendar(TIME,-1)/100
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
vwnd=short2flt(fv->vwnd(iYYYY,{850},{latS:latN},{lonL:lonR}))
vwnd=dim_standardize_n_Wrap(vwnd, 0, 0)
v =month_to_season(vwnd,"JJA")
printVarSummary(v)
printMinMax(v,False)
;*******************Write into one dataset*****************************
; Combine the normalized data into one variable ;
;**********************************************************************
dimw = dimsizes(precip)
mtim = dimw(0)
mlat = dimw(1)
mlon = dimw(2)
cdata = new ( (/3*mlat,3*mlon,mtim/), typeof(v), getFillValue(v))
do m1=0,mlat-1
do m2=0,mlon-1
cdata(m1, m2, :)=(/precip(:,m1,m2)/)
copy_VarMeta(precip,cdata(m1, m2, :) )
cdata(m1+mlat, m2+mlon, :)=(/u (:,m1,m2)/)
copy_VarMeta(precip, cdata(m1+mlat, m2+mlon, :))
cdata(m1+2*mlat,m2+2*mlon,:)=(/v (:,m1,m2)/)
copy_VarMeta(precip, cdata(m1+2*mlat,m2+2*mlon,:))
end do
end do
;***************************************t*****************************
; Compute **combined** EOF; Sign of EOF is arbitrary ;
;**********************************************************************
eof_cdata = eofunc(cdata , neof, False)
eof_ts_cdata = eofunc_ts(cdata,eof_cdata,False)
printVarSummary(eof_cdata)
printVarSummary(eof_ts_cdata)
nvar=3 ;precip,u,v
ceof = new( (/nvar,neof,mlat,mlon/), typeof(cdata), getFillValue(cdata))
do n=0,neof-1
ceof(0,n,:,:) = eof_cdata(n,0:mlat-1,0:mlon-1) ; precip
ceof(1,n,:,:) = eof_cdata(n,mlat:mlat*2-1,mlon:2*mlon-1) ; u
ceof(2,n,:,:) = eof_cdata(n,2*mlat:,2*mlon:) ; v
end do
ceof!0 = "var"
ceof!1 = "eof"
ceof!2 = "lat"
ceof!3 = "lon"
ceof&lat = precip&lat
ceof&lon = precip&lon
printVarSummary(ceof)
printMinMax(ceof(0,1,:,:),False)
;***************************************t*****************************
; plot now ;
;**********************************************************************
wks=gsn_open_wks("pdf","/disk1/dzz/study/ncl/TT/MVEOF")
plot = new(neof,graphic)
plot2= new(neof,graphic)
gsn_define_colormap(wks,"MPL_RdYlGn")
res = True
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
; res@gsnPolar = "SH"
res@gsnAddCyclic = False
res@gsnMaximize = True
res@mpFillOn = True ; turn off map fill
res@mpOutlineOn = False
res@mpMinLatF = 0
res@mpMaxLatF = 50
res@mpMaxLonF = 140
res@mpMinLonF = 100
res@cnFillOn = True ; turn on color fill
;res@cnFillPalette = "BlueWhiteOrangeRed"
res@cnLinesOn = False ; True is default
res@cnLineLabelsOn = False ; True is default
res@lbLabelBarOn = False ; turn off individual lb's
;res@cnLevelSelectionMode = "ManualLevels"
;res@cnMaxLevelValF = -0.02
;res@cnMinLevelValF = 0.02
;res@cnLevelSpacingF = 0.005
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = (/-0.02,-0.015,-0.01,-0.005,0,0.005,0.01,0.015,0.02/)
res@cnFillColors = (/5,25,38,56,0,0,80,94,108,121/)
; panel plot only resources
vcres = True
vcres@gsnDraw = False
vcres@gsnFrame = False
vcres@vcRefLengthF =0.045
vcres@vcMinDistanceF =0.017
;vcres@vcGlyphStyle = "CurlyVector"
resP = True ; modify the panel plot
resP@gsnMaximize = True ; large format
resP@gsnPanelLabelBar = True ; add common colorbar
resP@gsnPaperOrientation = "portrait" ; force portrait
resP@txString = ""
;*******************************************
; first plot
;*******************************************
do n=0,neof-1
res@tiMainString = "MVEOF "+(n+1)
res@gsnLeftString = "precipitation and 850hPa winds"
res@gsnRightString = ""
plot(n) =gsn_csm_contour_map(wks,ceof(0,n,:,:),res)
plot2(n)=gsn_csm_vector (wks,ceof(1,n,:,:), ceof(2,n,:,:), vcres)
overlay(plot(n),plot2(n))
draw(plot(n))
frame(wks)
end do
gsn_panel(wks,plot,(/neof,1/),resP) ; now draw as one plot
end

复制代码

代码如上，亲测可用，重复的图为B.Wang于2008年的一篇文章“How to Measure the Strength of the East Asian Summer Monsoon”，附上MV-EOF1-2和自己用上述程序出的MV-EOF1-2(关于原图和出图略有所差别是色标定义的问题（误差可忽略）以及没有平滑的缘故，但是整体上可表现出原图的模态)。
  接下来，就是我在画图中遇到的两个问题：
  问题1：一开始我将三个变量标准化然后写入一个矩阵的时候，是按全球来写的，接着对该矩阵进行EOF变换，这样出来的图和原图差别较大，之后我用东亚夏季风区域（0-50N，100-140E），也就是原图的边界，然后再写入矩阵作EOF的时候就是和原图一样的了，所以我的问题是：NCL中的EOF函数是否和数据的经纬度信息有关？
  问题2：大家应该能看出来，我的底图和填色图之间有空隙，也就是填色图不能完全覆盖底图，我觉得问题出在数据这快，我的precip降水数据是（time,lat(72),lon(144)）,然后uwnd和vwndde 数据是（time,lev,lat(73),lon(144)），也就是说降水的纬度比风场少一个格点，然后我在写入矩阵的时候，该新的矩阵维数信息定义为降水维数，所以出现这种填色图和底图有空隙，但是我要是将矩阵维数信息定义为风场维数的话，系统又会提示我：cdata(m1,    m2,    :)=(/precip(:,m1,m2)/)  这一行维数信息不对，我不知道该怎么解决，希望能有高人指点下，谢谢大家。
（风场的气旋和反气旋的“A”，“C”我就不再标上了，因为这个较为简单，官网一查即可。）