- 积分
- 20209
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-4-24
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
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”我就不再标上了,因为这个较为简单,官网一查即可。)
|
-
原图:FIG. 2. Spatial patterns of the (a) first and (b) second MV-EOF modes of the East Asian summer ...
-
自己做的图:mveof1
-
MVEOF2
|