- 积分
- 1132
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2012-11-2
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 guizi 于 2014-3-10 15:59 编辑
最近两天在用 ncl 对 HADISST 做EOF 分解
遇到个 神器错误，并且用ncl 老遇到神器错误，
这个 不解就是 加拿大北部，巴芬岛 那块的数据有问题，把范围选小点 就可以，
下面的脚本，发过一次，ncl脚本，稍作修改，就贡献在ncl 板块了
脚本如下
;************************************************
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl"
;************************************************
begin
;************* set up some parameter *********
lonL = 100
lonR = 240
latB = -10
latT = 40
yrStrt = 1950
yrLast = 2013
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
neof = 8
;************ open files and read the data ***********
f = addfile("../HadISST_sst.nc","r")
time = f->time ; days since 0000-09-01 00:00:00
YYYY = cd_calendar(time,-1)/100 ; entire file
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
; sst = f->sst(iYYYY,{latT:latB},:)
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
g_sst = f->sst(iYYYY,:,:)
g1_sst = lonPivot(g_sst,0.5)
sst = g1_sst(:,{latB:latT},{lonL:lonR})
;************* deal with the data *****
sst = where(abs(sst).gt.500,sst@_FillValue,sst)
printVarSummary(sst)
printMinMax(sst,True)
f_sst = sst
printVarSummary(f_sst)
printMinMax(f_sst,True)
ntime = dimsizes(sst&time)
ny = dimsizes(sst&latitude)
nx = dimsizes(sst&longitude)
;******************************************************
copy_VarCoords(sst,f_sst)
; printVarSummary(f_sst)
; y_sst = dtrend_msg_n(f_sst&time,f_sst,True,False,0)
; copy_VarCoords(sst,y_sst)
y_sst = dtrend_leftdim(f_sst,False)
; printVarSummary(y_sst)
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
temp = y_sst
nt = ntime/12
print(nt)
month = 12
temp_mon = new((/month,nt,ny,nx/),float)
temp_d = new((/ntime,ny,nx/),float)
do imon=0,month-1
temp_mon(imon,:,:,:) = dim_rmvmean_n_Wrap(temp(imon:ntime-1:12,:,:),0)
end do
do it=0,nt-1
do imon=0,month-1
temp_d(it*12+imon,:,:) = temp_mon(imon,it,:,:)
end do
end do
copy_VarCoords(temp,temp_d)
;~~~~~~~~~~~~~~~~~~~~~~~~~~~
x = temp_d
printVarSummary(x)
;*******************************************
; EOF
;*******************************************
X = x(latitude|:,longitude|:,time|:) ; Space x Time
optEof = True
eof = eofunc_Wrap( X, neof, optEof)
eof_ts = eofunc_ts_Wrap( X, eof, False)
printVarSummary(eof)
printVarSummary(eof_ts)
yyyymm = cd_calendar(eof_ts&time,-2)/100
year = yyyymm/100
;*******************************************
; plots
;*******************************************
wks = gsn_open_wks("pdf","pdo_eof")
gsn_define_colormap(wks,"posneg_1") ; choose colormap
;
; This will not be necessary in V6.1.0 and later. Named colors can
; be used without having to first add them to the color map.
;
dum = NhlNewColor(wks,0.7,0.7,0.7) ; add gray to color table
plot = new(neof,graphic) ; create graphic array
; only needed if paneling
res = True
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
res@gsnSpreadColors = True ; spread out color table
res@gsnSpreadColorEnd = -2 ; don't use added gray
res@gsnAddCyclic = False ; data not cyclic
res@mpCenterLonF = 190. ; defailt is 0 [GM]
res@mpMinLatF = min(x&latitude)
res@mpMaxLatF = max(x&latitude)
res@mpMinLonF = min(x&longitude)
res@mpMaxLonF = max(x&longitude)
res@mpFillDrawOrder = "PostDraw"
res@cnFillOn = True ; turn on color fill
res@cnLinesOn = True ; 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
resP@txString = "Kaplan_sst (Pierce)"
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(0:3),(/2,2/),resP) ; draw all 'neof' as one plot
;*******************************************
; time series (principal component) plot
;*******************************************
eof_ts@long_name = "Amplitude"
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 resources 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@gsnYRefLine = 0. ; reference line
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
; year = cd_calendar(eof_ts&time,1)
resP@txString = "KAplan_sst: Pierce"
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(0:3),(/2,2/),rtsP) ; draw all 'neof' as one plot
end
自己还有 上面的结果还是不好看，主要在 时间系数的 坐标问题
图正常，
要是 这么设
，图不正常，坐标正常
谁还有个好的法子撒？？？
？？？？？
|
-
|
窥视卡
雷达卡
发表于 2014-3-10 15:59:38
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
发表于 2014-3-10 16:40:35
楼主
谢谢 搞定了，