【NCL】ENSO冷暖事件海温异常合成与t检验

本代码参考了施宁老师《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)

合成效果如上图所示，结果尚不完美，有待改进(～￣▽￣)～