热带SST的PC1与SST做回归并检验问题

我是小猪

如图，我现在出的图填色是用student_t对回归检验的P值，不知道为什么会在陆地上有值，与原图差别较大？
其次想知道如何对t检验，负值填色冷色，正值填色暖色，我用的student_t 计算的值都是正值
代码如下(eof_ts的图并没有画出来)
希望得到解答，在这谢谢各位啦！！

1. 
   
2. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   
3. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
   
4. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
   
5. 
   
6. ; ==============================================================
   
7. ; User defined parameters that specify region of globe and
   
8. ; ==============================================================
   
9. begin
   
10.   latS   = -20.
    
11.   latN   =  20.
    
12.   yrStrt = 1979
    
13.   yrLast = 1996
    
14. 
    
15.   var    = "sst"
    
16.   title  = "Regr.on SST's PC1"
    
17. 
    
18.   ymStrt = yrStrt*100 +  1
    
19.   ymLast = yrLast*100 + 12
    
20. 
    
21.   neof   = 1                                  ; Leading EOF only
    
22.   optEOF = True
    
23.   optEOF@jopt = 1      
    
24.   optETS = False
    
25. 
    
26. ; ==============================================================
    
27. ; Open the file: Read only the user specified period and level
    
28. ; ==============================================================
    
29.   f      = addfile ("/disk3/HadISST/HadISST_sst.nc","r")
    
30. 
    
31.   TIME   = f->time
    
32.   YYYY   = cd_calendar(TIME,-1)/100                 ; entire file
    
33.   iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
    
34. 
    
35.   x      = short2flt( f->sst(iYYYY,{latS:latN},:))
    
36.   y      = short2flt(f->sst(iYYYY,:,:))
    
37.   printVarSummary(x)                                ; (time, lat,lon)
    
38. 
    
39. ; ==============================================================
    
40. ; compute climatology and Anomalies
    
41. ; ==============================================================
    
42.   
    
43.   xAnom=month_to_season(x, "DJF")
    
44.   yAnom=month_to_season(y, "DJF")
    
45.   ;xAnom=xAnom1/100.
    
46.   ;copy_VarMeta(xAnom1, xAnom)
    
47.   printVarSummary(xAnom)         
    
48.   printMinMax(xAnom, True)
    
49. 
    
50. ; =================================================================
    
51. ; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
    
52. ; =================================================================
    
53.   rad    = get_r2d("float")
    
54.   clat   = xAnom&latitude           
    
55.   clat   = sqrt( cos(rad*clat) )                 ; gw for gaussian grid
    
56.   printVarSummary(clat)
    
57. 
    
58. ; =================================================================
    
59. ; weight all observations
    
60. ; =================================================================
    
61.   xw     = xAnom*conform(xAnom, clat, 1)
    
62.   copy_VarMeta(xAnom, xw)
    
63.   xw@long_name = "Wgt: "+x@long_name
    
64. 
    
65. ; =================================================================
    
66. ; Reorder (lat,lon,time) the *weighted* input data
    
67. ; Compute EOFs & Standardize time series
    
68. ; =================================================================
    
69.   eof    = eofunc_n_Wrap(xw, neof, optEOF, 0)      
    
70.   eof    = -1*eof                                ; *special* match sign of CPC
    
71.   eof_ts = eofunc_ts_n_Wrap (xw, eof, optETS, 0)
    
72. 
    
73. 
    
74.   eof2    = eofunc_n_Wrap(yAnom, neof, optEOF, 0)      
    
75.   eof2    = -1*eof2
    
76.   printVarSummary( eof )                         ; examine EOF variables
    
77.   printVarSummary( eof_ts )
    
78. 
    
79.   eof_ts = dim_standardize_n( eof_ts, 0, 1)      ; normalize
    
80.   printMinMax(eof_ts, False)
    
81. ; =================================================================
    
82. ; Regress
    
83. ; =================================================================
    
84.   yw2=yAnom(latitude|:,longitude|:,time|:)
    
85.   tval =new((/180,360/),"float")
    
86.   nxy  =new((/180,360/),"integer")
    
87.   lon = lonGlobeF(360, "lon", "lon", "degrees_east")
    
88.   tval!1= "lon"
    
89.   tval&lon = lon
    
90.   lat = latGlobeF(180, "lat", "lat", "degrees_north")
    
91.   tval!0= "lat"
    
92.   tval&lat = lat
    
93. 
    
94.   lon = lonGlobeF(360, "lon", "lon", "degrees_east")
    
95.   nxy!1= "lon"
    
96.   nxy&lon = lon
    
97.   lat = latGlobeF(180, "lat", "lat", "degrees_north")
    
98.   nxy!0= "lat"
    
99.   nxy&lat = lat
    
100.   eof_regres = eof2                               ; create an array w meta data
     
101.   do ne=0,neof-1
     
102.      eof_regres(ne,:,:) = (/ regcoef(eof_ts(ne,:), yw2, tval,nxy) /)/100.
     
103.   end do
     
104.   df=16
     
105.   tval=student_t(tval, df)
     
106.   printVarSummary(tval)
     
107.   printMinMax(tval, False)
     
108.   printMinMax(eof_regres, False)
     
109. ; =================================================================
     
110. ; Extract the YYYYMM from the time coordinate
     
111. ; associated with eof_ts [same as x&time]
     
112. ; =================================================================
     
113. 
     
114.   yyyymm = cd_calendar(eof_ts&time,-1)  
     
115.   yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here
     
116. 
     
117. ;============================================================
     
118. ; PLOTS
     
119. ;============================================================
     
120.   wks = gsn_open_wks("png","/disk1/dzz/study/ncl/TT/ssteof")         ; send graphics to PNG file
     
121.   gsn_define_colormap(wks,"amwg256")
     
122.   plot = new(neof,graphic)                ; create graphic array
     
123.   plot2= new(neof,graphic)                                  ; only needed if paneling
     
124. ; EOF patterns
     
125. 
     
126.   res                      = True         
     
127.   res@gsnDraw              = False        ; don't draw yet
     
128.   res@gsnFrame             = False        ; don't advance frame yet
     
129. 
     
130.   res@mpFillOn             = False        ; turn off map fill
     
131.   res@mpOutlineOn          = True
     
132.   res@mpMaxLatF            = 60
     
133.   res@mpMinLatF            = -60
     
134.   res@mpMaxLonF            = 360
     
135.   res@mpMinLonF            = 0
     
136.   res@mpCenterLonF         = 180
     
137.   res@cnFillOn             = False         ; turn on color fill
     
138.   res@cnLinesOn            = True
     
139.   res@gsnContourNegLineDashPattern = 1
     
140.   res@lbLabelBarOn         = True        ; turn off individual lb's
     
141.   res@cnInfoLabelOn        = False
     
142.   res@lbTitleOn            = False
     
143.   res@cnLevelSelectionMode = "ManualLevels"
     
144.   res@cnLevelSpacingF             = 0.3
     
145.   ;res@cnFillPalette        = "amwg256"  
     
146.   res@cnLineDrawOrder      ="PostDraw"
     
147.   res@cnRasterSmoothingOn        = True                                  ; set symmetric plot min/max
     
148.   ;symMinMaxPlt(eof_regres, 0.3, False, res)       ; contributed.ncl
     
149. ; res@cnLevelSpacingF      = 0.3          ; *special* match CPC
     
150.   res2                       = True
     
151.   res2@gsnDraw               = False
     
152.   res2@gsnFrame              = False
     
153.   res2@cnFillOn              = True        ; turn on color fill
     
154.   res2@cnLinesOn            = False       ; True is default
     
155.   res2@cnLineLabelsOn       = False        ; True is default
     
156.   res2@lbLabelBarOn         = False        ; turn off individual lb's
     
157.   res2@cnLevelSelectionMode = "ExplicitLevels"
     
158.   res2@cnLevels             = (/0.05,0.01/)
     
159.   res2@cnFillColors         =  (/225,177,0/)
     
160.   res2@cnInfoLabelOn        = False
     
161.   res2@lbTitleOn            = False
     
162.   res2@cnSmoothingOn        = True
     
163. ; panel plot only resources
     
164. 
     
165.   resP                     = True         ; modify the panel plot
     
166.   resP@gsnMaximize         = True         ; large format
     
167.   resP@gsnPanelLabelBar    = True         ; add common colorbar
     
168.   ;resP@gsnPaperOrientation = "portrait"   ; force portrait
     
169. 
     
170.   resP@txString            = title
     
171. 
     
172. ;*******************************************
     
173. ; first plot
     
174. ;*******************************************
     
175.   do n=0,neof-1
     
176.      res@gsnLeftString  = "(c)1979-1996"
     
177.      res@gsnRightString = ""
     
178.      plot(n)=gsn_csm_contour_map(wks,eof_regres(n,:,:),res)
     
179.      plot2(n)=gsn_csm_contour    (wks,tval,res2)  
     
180.      overlay(plot, plot2)
     
181.      draw(plot)
     
182.      frame(wks)
     
183.   end do
     
184.   gsn_panel(wks,plot,(/neof,1/),resP)     ; now draw as one plot
     
185. end

复制代码

我是小猪

解决啦，主要是tval的经纬信息（lat,lon）与eof_regres的经纬信息（latitude，longitude）不一致，我统一改成（latitude，longitude）就可以了，但是不知道NCL为什么会这么识别

talkd

填色图错位了180度吧

我是小猪

另：EOF是（20S，20N）SST对其做EOF1，回归是全球SST与以上EOF1求得的标准化时间序列做的相关（再次感谢所有解答的人）

rpf

cos函数必须对应的是弧度吧，你用get_r2d函数将弧度值转化成度之后，再应用cos函数是错的吧？

我是小猪

talkd 发表于 2017-4-4 09:13
填色图错位了180度吧

哦 对对 这个倒是没注意，谢谢提醒！！

我是小猪

rpf 发表于 2017-4-4 11:16
cos函数必须对应的是弧度吧，你用get_r2d函数将弧度值转化成度之后，再应用cos函数是错的吧？

这个函数我是在官网计算EOF上看到的，不过我觉得 直接用eofunc函数来计算就可以，这个自带权重计算这几行不懂是为什么

zihan6688

学习了，谢谢大家讨论

deemo7

亲 修改的部分可以一下吗

deemo7

亲 修改的部分可以贴一下吗