- 积分
- 56
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2020-7-25
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2021-4-7 16:57:37
|
显示全部楼层
成功运行的程序如下:
begin
f = addfile("HadISST_sst.nc","r")
timeARR = cd_calendar(f->time,-5)
yrARR = timeARR(:,0)
timeIND = ind(yrARR.ge.1949 .and. yrARR.le.2019);IND948:1799
yrst = yrARR(timeIND);1949-2019,每个storm的年份
year = ispan(1949,2019,1)
sst0 = f->sst(timeIND,:,:) ; [time|852]x[latitude|180]x[longitude|360]
sstd = where(sst0.eq.-1000,sst0@_FillValue,sst0) ;min=-1.8 max=35.2085
copy_VarCoords(sst0,sstd)
m=dimsizes(sst0&latitude)
n=dimsizes(sst0&longitude)
sstdd = new((/71,m,n/),float)
do i=0,70 ;将一年的12个月做平均
term = 1949+i
tInd = ind(yrst.eq.term)
sstdd(i,:,:) = dim_avg_n_Wrap(sstd(tInd,:,:),0)
delete(tInd)
end do
sstAnom = dtrend_msg_n(year,sstdd,True,False,0)
;(包含sstdd第0维坐标的一维数组,去趋势的变量,是否去除均值,是否传回y的斜率截距att,去趋势的维)
copy_VarCoords(sstdd,sstAnom)
sstAnom!0 = "time" ;rmvmean函数需要变量有维名
;=======EOF分析===========
optEOF = True
optEOF@jopt = 0 ;1表示用相关矩阵计算EOFs,0表示用协方差矩阵计算(默认)
; 得到 EOF 模态
eof = eofunc_n_Wrap(sstAnom,10,optEOF,0) ;计算前10个模态;第四个参数dim,表示除dim维,其他维均参与EOF
; 投影得到主分量
pc = eofunc_ts_n_Wrap(sstAnom,eof,0,0) ;时间序列,第2,3个变量为optEOF,dim
; 提取 EOF1 模态
eofPlot = eof(0,:,:) ; 经过标准化
eofPlot = (/eofPlot*sqrt(eof@eval(0))/) ;去标准化,数值大小与原data一样,对应的时间系数除以相应特征值的开方
;标准后数值非原data,只能分析分布型,因此乘上相应特征值的开方去标准化
; 提取 PC1 的时间序列
pcPlot = dim_standardize(pc(0,:),0) ;[852]无缺省
;=======画图==================
wks = gsn_open_wks("png","eof")
res = True
res@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
res@gsnMaximize = False
;res@gsnAddCyclic = True ; 预设
res@vpWidthF = 0.7
res@vpHeightF = 0.6
res@mpCenterLonF = -150. ; 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@mpFillOn = True
res@mpLandFillColor = "grey"
res@mpOutlineOn = True
res@cnFillOn = True ; turn on color fill
res@cnLinesOn = False ; True is default
res@cnLineLabelsOn = False
res@cnFillPalette = "BlueWhiteOrangeRed"
res@cnInfoLabelOn = False
res@lbBoxEndCapStyle = "TriangleBothEnds"
;res@pmLabelBarOrthogonalPosF = 0.15
res@lbLabelBarOn = True
res@cnInfoLabelOn = False
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnFillColors = (/2,29,49,81,108,161,180,209,245/) ;颜色比数值多一个
res@cnLevels = (/-0.6,-0.5,-0.3,-0.1,0,0.1,0.2,0.3/)
res@tiMainString = "EOF1, " + sprintf("%.2f",eof@pcvar(0)) + "%"
res@tiMainFontHeightF = 0.026
res@lbLabelFontHeightF = 0.018
;=======pcPlot===============
resPC = True
resPC@gsnDraw = False
resPC@gsnFrame = False
resPC@gsnMaximize = False
resPC@vpWidthF = 0.7
resPC@vpHeightF = 0.3
;resPC@vpXF = .12
resPC@gsnCenterString = "PC1"
resPC@trXMinF = 1949
resPC@trXMaxF = 2019
resPC@trYMinF = -3.0
resPC@trYMaxF = 3.0
resPC@xyMarkLineMode = "Lines"
resPC@gsnXYBarChart = True
resPC@gsnYRefLine = 0.
resPC@gsnAboveYRefLineBarColors = "red"
resPC@gsnBelowYRefLineBarColors = "blue"
resPC@gsnCenterStringFontHeightF = 0.028
; panel plot only resources
resP = True
resP@gsnMaximize = True
;resP@gsnPanelXWhiteSpacePercent = 2.5
resP@gsnPanelYWhiteSpacePercent = 0.2
plot = new(2,graphic)
plot(0) = gsn_csm_contour_map_ce(wks,eofPlot,res)
plot(1) = gsn_csm_xy(wks,year,pcPlot,resPC)
gsn_panel(wks,plot,(/2,1/),resP)
end |
|