- 积分
- 577
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-2-27
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
我的目的是出pdo指数。这是我查看的一篇文献给的的计算方法。
利用全球月平均海表温度资料,对北太平洋 20°~60° N、110° E~ 100° W 范围内求距平并进行 EOF分析,得到的第一时间序列即为 PDO 指数。(不过说的并不详细)
我手里的这份脚本也是对海温进行距平再做eof分析。不过他做的是大西洋nino。脚本如下
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl" ; 如需将台站资料插值至格点资料,则必须加载该ncl文件
load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
neval=3
optEOF=True
optEOF@jopt=0
optETS = False
latS=-30.5
latN=30.5
lonL=-70.5
lonR=16.5
yrStrt = 1979
yrLast = 2014
;读取数据
f = addfile ("/cygdrive/d/ex/data/HadISST_sst.nc", "r")
;选取1979-2014数据
TIME = f->time
YYYY = cd_calendar(TIME,-1)/100
iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
sst = f->sst(iYYYY,{latS:latN},{lonL:lonR})
;printVarSummary(sst)
; sst = lonFlip( sst )
;夏季距平
sst=where(sst.lt.0,sst@_FillValue,sst)
;sst1=dim_rmvmean_n_Wrap(sst,0)
SST1=month_to_season(dim_rmvmean_n_Wrap(sst,0),"JJA")
SST=dtrend_leftdim(SST1,False) ;去除线性增长趋势
;权重
w=sqrt(cos(0.01745329*SST&latitude))
wp=SST*conform(SST,w,1)
copy_VarMeta(SST,wp)
data=wp({latitude|latS:latN},{longitude|lonL:lonR},time|:);调整维度
;EOF分析
ev1=eofunc_Wrap(data,neval,optEOF)
ts1=eofunc_ts_Wrap(data,ev1,optETS)
ev=10*ev1
ts=ts1/10
copy_VarMeta(ev1,ev)
copy_VarMeta(ts1,ts)
;空间数列
wks=gsn_open_wks("eps","Fig1a")
res=True
res@gsnAddCyclic = False
res@cnFillOn=False
res@cnLinesOn=True
res@lbBoxLinesOn = False
res@cnSmoothingOn = True
res@cnSmoothingDistanceF = 0.005
res@gsnSpreadColors = True
res@mpMinLatF=latS
res@mpMaxLatF=latN
res@mpMinLonF = lonL
res@mpMaxLonF = lonR
res@cnLineLabelsOn=True
res@cnLineLabelPlacementMode="Constant"
res@cnLevelSelectionMode="ExplicitLevels"
res@cnLevels=(/-0.2,-0.1,0.0,0.1,0.2,0.3,0.4/)
;res@cnLineLabelDensityF=2
res@cnLineLabelInterval =1
res@cnInfoLabelOn=False
res@gsnContourNegLineDashPattern=16
res@gsnContourPosLineDashPattern=0
res@cnLineDashSegLenF = 0.15
res@cnInfoLabelOn = False ;去除右下角contour from to
res@cnLineLabelFontHeightF = 0.01
res@cnLineLabelFormat ="@*+^sg" ;数值形式 eg:0.1
res@pmTickMarkDisplayMode = "Always"
res@gsnLeftString = "EOF "+(0+1)
res@gsnRightString = sprintf("%5.1f", ev@pcvar(1)) +"%"
symMinMaxPlt(ev, 16, False, res)
plot=gsn_csm_contour_map_ce(wks,smth9_Wrap(ev(1,:,:),0.5,0.25,False),res)
;时间序列
year=ispan(1979,2014,1)
wks_ts=gsn_open_wks("png","Fig1b")
rts = True
rts@gsnScale = True
rts@vpHeightF = 0.40
rts@vpWidthF = 0.85
rts@vpXF = 0.10
rts@vpYF = 0.75
rts@tiYAxisString = "~S~o~N~C"
rts@gsnYRefLine = 0.
rts@gsnXYBarChart = True
rts@gsnAboveYRefLineColor = "grey"
rts@gsnBelowYRefLineColor = "grey"
rts@gsnLeftString = "EOF "+(0+1)
rts@gsnRightString = sprintf("%5.1f", ev@pcvar(1)) +"%"
plot2=gsn_csm_xy(wks_ts,year,ts(1,:),rts)
;;;;;;;;;;;;;;;;;;;;;
;时间序列标准化
ts=dim_standardize_n_Wrap(ts,1,1)
;
end
然后我再讲一个我妄图改程序结果失败的故事。
这份脚本中距平部分做的是季节平均,只选择了一个季节,就相当与是做的年分析。曾经我想改这个程序,直接把他做季节距平的函数部分改成了做距平的函数。结果就由于改动后的距平是逐月的数据,而year是年的。所以
这一句当时就报错了。我的一次瞎改也失败了。
真的急啊!!!希望来个前辈帮帮我!可以有偿的!
微信:lrc756866304 qq756866304
希望有人愿意帮忙,资料我全都有! |
|