本代码参考了施宁老师《NCL数据处理与绘图实习教程》一书,按实际情况有所增删。主要效果为合成ENSO年厄尔尼诺和拉尼娜海温异常,并进行t检验,将置信度为0.05的区域打点。(由于没有专门的NCL高亮代码块,特别注意分号后面全部为注释部分)。
另外,代码书写过程中断断续续,导致每个步骤独立性很强。重复进行的步骤可以缩减,例如在绘图部分可以直接res2=True/res2=res1,而后更改图名即可。错处欢迎指正~~
year=ispan(1940,2018,1)
it_s=194012;;海温数据
it_e=201811
;;;;;;;;;;;;;;;;;;;;;;;;;;;read data;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;HadIsst——1870-2019
f_sst=addfile("/data/home//Data/SST/sst.mon.mean.nc","r")
time =f_sst->time;;读取日期数据
YYYYMM=cd_calendar(time,-1);;转换为公历日期
rec_s=ind(it_s.eq.YYYYMM);;开始日期的记录号
rec_e=ind(it_e.eq.YYYYMM);;终止日期的记录号
sst=f_sst->sst
sst_34=sst(rec_s:rec_e,:,:);;截取开始——终止的时段资料
sst_DJF=month_to_season(sst_34, "JFM");求JFM的季节平均
; printVarSummary(sst_34)
copy_VarMeta(sst_34(0,:,:),sst_DJF(0,:,:));copy_VarMeta(var_from, var_to)将前者的属性给后者
sst_DJF!0="year";;;第一维的变量名字是year
sst_DJF&year=year;;
;printVarSummary(sst_DJF)
sst_ano=dim_rmvmean_n_Wrap(sst_DJF,0)
copy_VarCoords(sst_DJF, sst_ano)
ensoi=wgt_areaave(sst_DJF(:,{-5:5},{190:240}),1.0,1.0,0)
;;0表示仅用非缺省数值计算
;;标准化函数
ensoi=dim_standardize(ensoi, 1)
; ;;1表示标准化时除以[N],1表示除以[N-1]
;printVarSummary(ensoi)
;;;;;;;;;;;;;;;;;;;;;;composite;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
irec_positive=ind(ensoi.gt.0.8)
irec_negative=ind(ensoi.lt.-0.8)
;irec_negative=ind(ensoi.gt.-0.8);;
nnumb=dimsizes(irec_positive)
nnuma=dimsizes(irec_negative)
nino_comp=dim_avg_n_Wrap(sst_ano(irec_positive,:,:),0)
nina_comp=dim_avg_n_Wrap(sst_ano(irec_negative,:,:),0)
copy_VarCoords(sst_ano(0,:,:), nino_comp)
copy_VarCoords(sst_ano(0,:,:), nina_comp)
;printVarSummary(sst_comp)
;;;;;;;;;;;;;;;;;;;(t-test);;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;sst
nino_std=dim_variance_n_Wrap(sst_ano(irec_positive,:,:), 0)
nina_std=dim_variance_n_Wrap(sst_ano(irec_negative,:,:), 0)
nino_std=sqrt(nino_std/nnumb);;标准差
nino_std=where(nino_std.eq.0,nino_std@_FillValue,nino_std);
nina_std=sqrt(nina_std/nnuma);;标准差
nina_std=where(nina_std.eq.0,nina_std@_FillValue,nina_std)
t_nino=nino_comp/nino_std;;;
confi_nino=nino_comp
confi_nino=student_t(t_nino, nnumb-1)
t_nina=nina_comp/nina_std
confi_nina=nina_comp
confi_nina=student_t(t_nina, nnuma-1)
;print(sst_comp)
copy_VarCoords(nino_comp, confi_nino)
copy_VarCoords(nina_comp, confi_nina)
;printVarSummary(confi_sst)
;;;第一种方法分EP和CP
enso3i=wgt_areaave(sst_DJF(:,{-5:5},{210:270}),1.0,1.0,0)
enso3i_s=dim_standardize(enso3i, 1)
enso4i=wgt_areaave(sst_DJF(:,{-5:5},{160:210}),1.0,1.0,0)
enso4i_s=dim_standardize(enso4i, 1)
irec_EP=ind((enso3i_s.gt.1).and.(enso3i.gt.enso4i))
irec_CP=ind((enso4i_s.gt.1).and.(enso4i.gt.enso3i))
nnume=dimsizes(irec_EP)
nnumc=dimsizes(irec_CP)
EP_comp=dim_avg_n_Wrap(sst_ano(irec_EP,:,:),0)
CP_comp=dim_avg_n_Wrap(sst_ano(irec_CP,:,:),0)
copy_VarCoords(sst_ano(0,:,:), EP_comp)
copy_VarCoords(sst_ano(0,:,:), CP_comp)
EP_std=dim_variance_n_Wrap(sst_ano(irec_EP,:,:), 0)
CP_std=dim_variance_n_Wrap(sst_ano(irec_CP,:,:), 0)
EP_std=sqrt(EP_std/nnume);;标准差
EP_std=where(EP_std.eq.0,EP_std@_FillValue,EP_std);;;
CP_std=sqrt(CP_std/nnumc);;标准差
CP_std=where(CP_std.eq.0,CP_std@_FillValue,CP_std)
t_EP=EP_comp/EP_std;;;
confi_EP=EP_comp
confi_EP=student_t(t_EP, nnume-1)
t_CP=CP_comp/CP_std
confi_CP=CP_comp
confi_CP=student_t(t_CP, nnumc-1)
copy_VarCoords(EP_comp, confi_EP)
copy_VarCoords(CP_comp, confi_CP)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
wks=gsn_open_wks("PDF", "/data/home/yuancx/wangjy/nino_compo_test/try6")
gsn_define_colormap(wks, "MPL_coolwarm")
plot=new(4,"graphic")
base=new(4,"graphic")
res =True
res@gsnAddCyclic =True
res@gsnDraw =False
res@gsnFrame =False
res@gsnLeftString =""
res@gsnRightString =""
res@mpCenterLonF =180;中心经度为180°
res@mpGeophysicalLineThicknessF=0.2;地图边界
res@mpFillOn =False;地图不填色
res@mpGridAndLimbOn =True;绘制经纬度线
res@mpGridLatSpacingF =15;纬度线的间隔
res@mpGridLonSpacingF =15;经度线的间隔
res@mpGridLineDashPattern =2
res@mpGridLineThicknessF =0.05;经纬线的粗细
res@cnFillOn =True;填充图
res@lbLabelBarOn =True;色标打开
res@lbOrientation ="Horizontal"
res@cnLinesOn =False;不要等值线怪丑的
res@cnLineLabelsOn =False;等值线上不要标签
res@cnLevelSelectionMode ="ExplicitLevels"
res@cnLevels =(/-1.4,-1.2,-1.0,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1.0,1.2,1.4/)
res@cnLineThicknessF =2
res@cnInfoLabelOn =False
res@gsnCenterString ="El-nino"
res@gsnCenterStringFontHeightF =0.02
plot(0)=gsn_csm_contour_map(wks, nino_comp, res)
res3 =True
res3@gsnAddCyclic =True
res3@gsnDraw =False
res3@gsnFrame =False
res3@gsnLeftString =""
res3@gsnRightString =""
res3@mpCenterLonF =180;中心经度为180°,全球图要从20度经度开始,后面注意!
res3@mpGeophysicalLineThicknessF=0.2;地图边界
res3@mpFillOn =False;地图不填色
res3@mpGridAndLimbOn =True;绘制经纬度线
res3@mpGridLatSpacingF =15;纬度线的间隔
res3@mpGridLonSpacingF =15;经度线的间隔
res3@mpGridLineDashPattern =2
res3@mpGridLineThicknessF =0.05
res3@cnFillOn =True;填充图
res3@lbLabelBarOn =True;色标打开
res3@lbOrientation ="Horizontal"
res3@cnLinesOn =False;不要等值线怪丑的
res3@cnLineLabelsOn =False;等值线上不要标签
res3@cnLevelSelectionMode ="ExplicitLevels"
res3@cnLevels =(/-1.4,-1.2,-1.0,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1.0,1.2,1.4/)
res3@cnLineThicknessF =2
res3@cnInfoLabelOn =False
res3@gsnCenterString ="La-nina"
res3@gsnCenterStringFontColor =0.02
plot(1)=gsn_csm_contour_map(wks, nina_comp, res3)
res1 =True
res1@gsnAddCyclic =True
res1@gsnDraw =False
res1@gsnFrame =False
res1@gsnLeftString =""
res1@gsnRightString =""
res1@mpCenterLonF =180;中心经度为180°
res1@mpGeophysicalLineThicknessF=0.2;地图边界
res1@mpFillOn =False;地图不填色
res1@mpGridAndLimbOn =True;绘制经纬度线
res1@mpGridLatSpacingF =15;纬度线的间隔
res1@mpGridLonSpacingF =15;经度线的间隔
res1@mpGridLineDashPattern =2
res1@mpGridLineThicknessF =0.05;经纬线的粗细
res1@cnFillOn =True;填充图
res1@lbLabelBarOn =True;色标打开
res1@lbOrientation ="Horizontal"
res1@cnLinesOn =False;不要等值线怪丑的
res1@cnLineLabelsOn =False;等值线上不要标签、
res1@cnLevelSelectionMode ="ExplicitLevels"
res1@cnLevels =(/-1.4,-1.2,-1.0,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1.0,1.2,1.4/)
res1@cnLineThicknessF =2
res1@cnInfoLabelOn =False
res1@gsnCenterString ="EP(nino3>1std&nino3>nino4)"
res1@gsnCenterStringFontHeightF =0.02
plot(2)=gsn_csm_contour_map(wks, EP_comp, res1)
res2 =True
res2@gsnAddCyclic =True
res2@gsnDraw =False
res2@gsnFrame =False
res2@gsnLeftString =""
res2@gsnRightString =""
res2@mpCenterLonF =180;中心经度为180°
res2@mpGeophysicalLineThicknessF=0.2;地图边界
res2@mpFillOn =False;地图不填色
res2@mpGridAndLimbOn =True;绘制经纬度线
res2@mpGridLatSpacingF =15;纬度线的间隔
res2@mpGridLonSpacingF =15;经度线的间隔
res2@mpGridLineDashPattern =2
res2@mpGridLineThicknessF =0.05;经纬线的粗细
res2@cnFillOn =True;填充图
res2@lbLabelBarOn =True;色标打开
res2@lbOrientation ="Horizontal"
res2@cnLinesOn =False;不要等值线怪丑的
res2@cnLineLabelsOn =False;等值线上不要标签、
res2@cnLevelSelectionMode ="ExplicitLevels"
res2@cnLevels =(/-1.4,-1.2,-1.0,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1.0,1.2,1.4/)
res2@cnLineThicknessF =2
res2@cnInfoLabelOn =False
res2@gsnCenterString ="CP(nino4>1std&nino4>nino3)"
res2@gsnCenterStringFontHeightF =0.02
plot(3)=gsn_csm_contour_map(wks, CP_comp, res2)
;;;打点部分施工中——请出示健康码;;;;;;;
opt =True
opt@gsnShadeFillType ="pattern"
opt@gsnShadeLow =17;打点
opt@gsnAddCyclic =True
opt@cnFillDotSizeF =0.001;改变点大小的默认值
sres =True
sres@gsnDraw =False
sres@gsnFrame =False
sres@cnLinesOn =False
sres@gsnLeftString =""
sres@gsnRightString =""
sres@cnLevelSelectionMode="ExplicitLevels"
sres@cnLevels =(/0.05,0.01/)
sres@cnFillPalette ="GMT_gray"
sres@cnFillColors =(/5,7,-1/)
sres@cnLineLabelsOn =False
sres@cnInfoLabelOn =False
sres@lbLabelBarOn =False
base(0)=gsn_csm_contour(wks,confi_nino, sres)
base(0)=gsn_contour_shade(base(0),0.05,1,opt)
overlay(plot(0),base(0))
base(1)=gsn_csm_contour(wks,confi_nina, sres)
base(1)=gsn_contour_shade(base(1),0.05,1,opt)
overlay(plot(1),base(1))
base(2)=gsn_csm_contour(wks,confi_EP, sres)
base(2)=gsn_contour_shade(base(2),0.05,1,opt)
overlay(plot(2),base(2))
base(3)=gsn_csm_contour(wks,confi_CP, sres)
base(3)=gsn_contour_shade(base(3),0.05,1,opt)
overlay(plot(3),base(3))
resP =True
resP@txstring ="ENSO composite"
resP@txFontHeightF =0.03
resP@gsnPanelFigureStrings =(/"(a)","(b)","(c)","(d)","(e)","(f)","(g)","(h)"/)
resP@gsnPanelFigureStringsFontHeightF=0.008
resP@amJust ="TopLeft"
resP@gsnPanelRowSpec =True
gsn_panel(wks,plot,(/2,2,2,2/),resP)
合成效果如上图所示,结果尚不完美,有待改进,欢迎指正