- 积分
- 15746
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2015-3-8
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
;min_lat = -10.0
; max_lat =40.0
;min_lon = 90.0
; max_lon = 160.0
latS = -20
latN = 50
lonL =60
lonR =180
neof = 3 ; number of EOFs
optEOF = True
optEOF@jopt = 0 ; This is the default; most commonly used; no need to specify.
;;optEOF@jopt = 1 ; **only** if the correlation EOF is desired
optETS = False
diri = "H:/twentyc/"
fili1 = "uchunelnino.grd"
fili2 = "vchunelnino.grd"
f1=fbindirread(diri+fili1,0,(/29,91,180/),"float")
f2=fbindirread(diri+fili2,0,(/29,91,180/),"float")
nlat = 91 ; YDEF
mlon = 180 ; XDEF
year = 1878 ; TDEF
ntim = 29 ; time steps
nmos = 1
; not required
time = new (ntim, float ,"No_FillValue") ; generate a "time" variable
;date = new (ntim, integer,"No_FillValue") ; generate a "date" variable
time = ispan(0,ntim-1,1)*1. + 1878
;n = -1
;do nmo=1,nmos
; YRM= year*10000 + nmo*100
;ndm= days_in_month(year, nmo)
;do ndy=1,ndm
;n = n+1
;time(n) = n
;date(n) = YRM + ndy ; YYYYMMDD
;end do
;end do
time!0 = "time"
time@long_name = "time"
time@units = "yr" ; "yyyy.fraction_of_year"
time&time = time
;date!0 = "time"
;date@long_name = "date"
; date@units = "dy"
; date&time = time
; generate lat/lon
lon = ispan(0,mlon-1,1)*2. + 0.
lon!0 = "lon"
lon@long_name = "longitude"
lon@units = "degrees_east"
lat = ispan(0,nlat-1,1)*2. -90.
lat!0 = "lat"
lat@long_name = "latitude"
lat@units = "degrees_north"
;time=ispan(1878,1906,1)
; lat=fspan(-90,90,91)
; lon=fspan(0,360,180)
;; time@units="years since 1878-03-01 00:00:0.0"
; lat@units="degrees_north"
; lon@units="degrees_east"
u=f1(:,:,:)
v=f2(:,:,:)
u!0="time"
u!1="lat"
u!2="lon"
u&time=time
u&lat=lat
u&lon=lon
v!0="time"
v!1="lat"
v!2="lon"
v&time=time
v&lat=lat
v&lon=lon
printVarSummary(u)
printVarSummary(v)
;YYYY = cd_calendar(time,-1)
s1 = uv2vr_cfd (u(:,:,:),v(:,:,:),lat,lon, 3)
s1=s1*1000000
copy_VarMeta(u,s1)
printVarSummary(s1)
s=dtrend_leftdim(s1,False)
copy_VarMeta(s1,s)
s = lonFlip(s)
printVarSummary(s) ; note the longitude coord
rad = 4.*atan(1.)/180.
clat =lat
clat = sqrt( cos(rad*clat) ) ; gw for gaussian grid
; =================================================================
; weight all observations
; =================================================================
wlh = s ; copy meta data
wlh = s*conform(s, clat, 1)
; wlh@long_name = "Wgt: "+wlh@long_name
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
x = wlh({lat|latS:latN},{lon|lonL:lonR},time|:)
eof = eofunc_Wrap(x, neof, optEOF)
eof_ts = eofunc_ts_Wrap (x, eof, optETS)
printVarSummary( eof ) ; examine EOF variables
printVarSummary( eof_ts )
; =================================================================
; Normalize time series: Sum spatial weights over the area of used
; =================================================================
dimx = dimsizes( x )
mln = dimx(1)
sumWgt = mln*sum( lat(0:36) )
eof_ts = eof_ts/sumWgt
; =================================================================
; Extract the YYYYMM from the time coordinate
; associated with eof_ts [same as SLP&time]
; =================================================================
yyyy = cd_calendar(eof_ts&time,-2)/100
;yyyymm = cd_calendar(time,-2)/100
;;yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here
;============================================================
; PLOTS
;============================================================
wks = gsn_open_wks("ps","eofvor29")
gsn_define_colormap(wks,"BlWhRe") ; choose colormap
plot = new(neof,graphic) ; create graphic array
; only needed if paneling
; EOF patterns
res = True
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
;---This resource not needed in V6.1.0
res@gsnSpreadColors = True ; spread out color table
res@gsnAddCyclic = False ; plotted dataa are not cyclic
res@mpFillOn = False ; turn off map fill
res@mpMinLatF = latS ; zoom in on map
res@mpMaxLatF = latN
res@mpMinLonF = lonL
res@mpMaxLonF = lonR
res@cnFillOn = True ; turn on color fill
res@cnLinesOn = False ; True is default
;res@cnLineLabelsOn = False ; True is default
res@lbLabelBarOn = False ; turn off individual lb's
; set symmetric plot min/max
symMinMaxPlt(eof, 16, False, res) ; contributed.ncl
; panel plot only resources
resP = True ; modify the panel plot
resP@gsnMaximize = True ; large format
resP@gsnPanelLabelBar = True ; add common colorbar
resP@lbLabelAutoStride = True ; auto stride on labels
yStrt = yyyy(0)
yLast = yyyy(ntim-1)
; yStrt = 1878
; yLast = 1906
;resP@txString = "lh: "+season+": "+yStrt+"-"+yLast
;*******************************************
; first plot
;*******************************************
do n=0,neof-1
res@gsnLeftString = "EOF "+(n+1)
res@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
plot(n)=gsn_csm_contour_map_ce(wks,eof(n,:,:),res)
end do
gsn_panel(wks,plot,(/neof,1/),resP) ; now draw as one plot
;*******************************************
; second plot
;*******************************************
; EOF time series [bar form]
rts = True
rts@gsnDraw = False ; don't draw yet
rts@gsnFrame = False ; don't advance frame yet
rts@gsnScale = True ; force text scaling
; these four rtsources allow the user to stretch the plot size, and
; decide exactly where on the page to draw it.
rts@vpHeightF = 0.40 ; Changes the aspect ratio
rts@vpWidthF = 0.85
rts@vpXF = 0.10 ; change start locations
rts@vpYF = 0.75 ; the plot
rts@tiYAxisString = "" ; y-axis label
rts@gsnYRefLine = 0. ; reference line
rts@gsnXYBarChart = True ; create bar chart
rts@gsnAboveYRefLineColor = "red" ; above ref line fill red
rts@gsnBelowYRefLineColor = "blue" ; below ref line fill blue
; panel plot only resources
rtsP = True ; modify the panel plot
rtsP@gsnMaximize = True ; large format
;rtsP@txString = "lh: "+season+": "+yStrt+"-"+yLast
;year1 = yyyymm/100
; create individual plots
do n=0,neof-1
rts@gsnLeftString = "EOF "+(n+1)
rts@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
plot(n) = gsn_csm_xy (wks,year,eof_ts(n,:),rts)
end do
gsn_panel(wks,plot,(/neof,1/),rtsP) ; now draw as one plot
end
空间场没什么大问题，但是时间场出错：
我知道是eof_ts这块出错了，但是因为是grd的我也不知道应该怎么改，就自己试了好久还是出错，希望做过这方面的大神指点一下吧
|
|
窥视卡
雷达卡
发表于 2016-3-6 12:06:31
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
{:eb303:}{:eb303:}
发表于 2016-3-6 16:47:17
{:eb511:}