合成分析法，700百帕风场差异，t检验，以及过检验的地方打点

恦之兰泽

前面已经得到突变点为2003年，并且过了99%的显著水平。现在进行合成分析，可是我研究的地区却没有过检验，求教一下程序是否有问题。

1. #1979-2020年5-10月的月平均、700百帕水平风场资料
   
2. import numpy as np
   
3. import Ngl,Nio
   
4. 
   
5. f=Nio.open_file('/mnt/e/xinan/uv700_1979-2020_5-10.nc','r')
   
6. #print(f)
   
7. 
   
8. lon  = f.variables['longitude']            # (361,)
   
9. lat  = f.variables['latitude']             # (161,)
   
10. 
    
11. u    = f.variables['u']                    # (252,161,361)  
    
12. v    = f.variables['v']                    # (252,161,361)
    
13. u    = u*u.scale_factor+u.add_offset
    
14. v    = v*v.scale_factor+v.add_offset
    
15. 
    
16. #突变年份为2003年，分成1979-2003和2004-2020两段做合成分析
    
17. #(2003-1979+1)*6=150 & (2020-2004+1)*6=102
    
18. ave_u1  = np.mean(u[0:150][:][:],0)
    
19. ave_u2  = np.mean(u[150:252][:][:],0)
    
20. ave_v1  = np.mean(v[0:150][:][:],0)
    
21. ave_v2  = np.mean(v[150:252][:][:],0)
    
22. dif_u   = ave_u2-ave_u1
    
23. dif_v   = ave_v2-ave_v1
    
24. 
    
25. # t 检验
    
26. aveX = np.mean(v[0:150][:][:],0)
    
27. aveY = np.mean(v[150:252][:][:],0)
    
28. varX = np.var(v[0:150][:][:],0)
    
29. varY = np.var(v[150:252][:][:],0)
    
30. n1   = 25
    
31. n2   = 17
    
32. df   = n1+n2-2
    
33. d    = (n1-1)*varX+(n2-1)*varY
    
34. t    = (aveX-aveY)/(np.sqrt(d/df)*np.sqrt(1/n1+1/n2))
    
35. #print(t)_(161,361)
    
36. print(np.min(t),np.max(t))
    
37. 
    
38. #得到通过t检验的格点下标
    
39. t1=[]
    
40. t2=[]
    
41. for i in range(0,161):
    
42.     for j in range(0,361):
    
43.         if np.abs(t[i][j])>=2.021:       #自由度=40，α=0.05时，t=2.021
    
44.             t1.append(i)
    
45.             t2.append(j)
    
46. 
    
47. print(len(t1),len(t2))
    
48. 
    
49. #得到通过检验的格点经纬度
    
50. lat1=[]
    
51. lon1=[]
    
52. for i in range(0,len(t1)):
    
53.     lat2=lat[t1[i]]
    
54.     lon2=lon[t2[i]]
    
55.     lat1.append(lat2)
    
56.     lon1.append(lon2)
    
57. 
    
58. #print(len(lon1))
    
59. 
    
60. wks   = Ngl.open_wks('x11','/mnt/f/uv_700(ttest)')
    
61. 
    
62. vcres = Ngl.Resources()
    
63. mpres = Ngl.Resources()
    
64. 
    
65. mpres.nglDraw   =False
    
66. mpres.nglFrame  =False
    
67. vcres.nglDraw   =False
    
68. vcres.nglFrame  =False
    
69. 
    
70. vcres.vfXArray  =lon[::6]
    
71. vcres.vfYArray  =lat[::6]
    
72. 
    
73. mpres.mpLimitMode                 = "LatLon"       # Change the area of the
    
74. mpres.mpMinLatF                   = 10             # map viewed.
    
75. mpres.mpMaxLatF                   = 40
    
76. mpres.mpMinLonF                   = 85
    
77. mpres.mpMaxLonF                   = 120
    
78. 
    
79. #需要的省界
    
80. mpres.mpDataSetName               = "Earth..4"
    
81. mpres.mpDataBaseVersion           = 'MediumRes'
    
82. mpres.mpMaskAreaSpecifiers        = ["Yunnan","Sichuan","Guizhou"]
    
83. mpres.mpOutlineSpecifiers         = ["Yunnan","Sichuan","Guizhou"]
    
84. mpres.mpProvincialLineThicknessF  = 2
    
85. mpres.mpNationalLineThicknessF    = 2
    
86. 
    
87. mpres.mpOutlineBoundarySets       = "Geophysical"
    
88. mpres.mpGeophysicalLineThicknessF = 2.
    
89. 
    
90. #箭头设置
    
91. vcres.vcMinFracLengthF            = 0.35         # Increase length of
    
92. vcres.vcMinMagnitudeF             = 0.001         # vectors.
    
93. vcres.vcRefLengthF                = 0.045
    
94. vcres.vcRefMagnitudeF             = 2.0
    
95. vcres.vcGlyphStyle                = 'CurlyVector'
    
96. 
    
97. #题目
    
98. mpres.tiMainString                = "uv_700(ttest)"
    
99. 
    
100. map_plot            = Ngl.map(wks,mpres)
     
101. vector_plot         = Ngl.vector(wks,dif_u[::6,::6],dif_v[::6,::6],vcres)
     
102. Ngl.overlay(map_plot,vector_plot)
     
103. Ngl.draw(map_plot)
     
104. 
     
105. #标记通过检验的位置
     
106. mres                = Ngl.Resources()
     
107. mres.gsMarkerIndex  = 16
     
108. mres.gsMarkerSizeF  = 0.002
     
109. mres.gsMarkerColor  = 'black'
     
110. 
     
111. Ngl.polymarker(wks,map_plot,lon1,lat1,mres)
     
112. 
     
113. Ngl.frame(wks)

复制代码

Sharon_Liu

风场是矢量，包括u、v，这里只对v做t检验了，是不是不太合适呢？

灰色节能君

你似乎用的是两总体的t检验，这是检验特殊年和一般年均值差异的。在气象数据统计分析方法(2015)黄嘉佑的附录K9中给出了专门用于检验相邻时段突变的检验统计量

恦之兰泽

灰色节能君 发表于 2021-4-18 10:36
你似乎用的是两总体的t检验，这是检验特殊年和一般年均值差异的。在气象数据统计分析方法(2015)黄嘉佑的附 ...

这个检验方法应该是没问题的，和我一组的一个同学也是根据突变年份做合成分析，就是用的这种t检验，效果很好

SunJiaming

感谢楼主大大分享！

807727227

感谢，正好在找合成分析显著性检验的代码呢

文盲易

{:eb502:}{:eb502:}{:eb502:}

892282469

感谢楼主分享

puck66

后来呢？处理成啥样了，分享一下喽

qxtlyf

学习了，谢谢了