【NCL】ENSO冷暖事件海温异常合成与t检验_ncl 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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244

合成效果如上图所示，结果尚不完美，有待改进，欢迎指正