EOF 问题

qq469015280

对热带印度洋的SSTA距平场做EOF
程序：参考这里的第1个例子http://met.sysu.edu.cn/GloCli/Team/ncl-mirror/Applications/eof.shtml

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
  latS   =  -20.
  latN   =  20.
  lonL   = 40.
  lonR   = 110.
  yrStrt = 1982
  yrLast = 2015
  neof   =3       ; number of EOFs
  optEOF = True      
  optEOF@jopt = 0   ; This is the default; most commonly used; no need to specify.

  optETS = False
; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
  
     
      f2 =addfile ("ssta.nc", "r")               ;sst距平数据          ; variable overview
sst = f2->sst  
   printVarSummary(sst)  
slp    =  dim_avg_n_Wrap(sst(mon|:,time|:,lat|:,lon|:)  ,0)
print(slp)
   
  printVarSummary(slp)                              ; variable overview
  SLP=slp(time|:,lat|:,lon|:)
  printVarSummary(SLP)   
time=ispan(1,34,1)
time@units="years since 1982-01-01 00:00:0.0"
  SLP&time=time
  nyrs   = dimsizes(SLP&time)
  printVarSummary(SLP)
printVarSummary(nyrs)
print(nyrs)
; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
; =================================================================

   rad    = 4.*atan(1.)/180.
   
  clat   =lat  
  
  clat!0="lat"
    clat&lat=lat
  clat   = sqrt( cos(rad*clat) )                 ; gw for gaussian grid
; =================================================================
; weight all observations
; =================================================================
  wSLP   = SLP                                   ; copy meta data
  wSLP   = SLP*conform(SLP, clat, 1)
  wSLP@long_name = "Wgt: "+wSLP@long_name
; =================================================================
; Reorder (lat,lon,time) the *weighted* input data
; Access the area of interest via coordinate subscripting
; =================================================================
  xw     = wSLP({lat|latS:latN},{lon|lonL:lonR},time|:)
    printVarSummary( xw )
  x      = wSLP(time|:,{lat|latS:latN},{lon|lonL:lonR})
  eof      = eofunc_Wrap(xw, neof, optEOF)      
  eof_ts   = eofunc_ts_Wrap (xw, eof, optETS)
  printVarSummary( eof )                         ; examine EOF variables
  printVarSummary( eof_ts )
; =================================================================
; Normalize time series: Sum spatial weights over the area of used
; =================================================================
  dimxw  = dimsizes( xw )
  mln    = dimxw(1)
  sumWgt = mln*sum( clat({lat|latS:latN}) )
  eof_ts = eof_ts/sumWgt
; =================================================================
; Extract the YYYYMM from the time coordinate
; associated with eof_ts [same as SLP&time]
; =================================================================
  yyyymm = cd_calendar(eof_ts&time,-2)/100
;============================================================
; PLOTS

;============================================================
  wks = gsn_open_wks("ps","eof")
  plot = new(neof,graphic)                ; create graphic array
                                          ; only needed if paneling
; EOF patterns
  res                      = True         
  res@gsnDraw              = False        ; don't draw yet
  res@gsnFrame             = False        ; don't advance frame yet
  res@gsnAddCyclic         = False        ; plotted dataa are not cyclic

  res@mpFillOn             = True        ; turn off map fill
      res@mpMinLatF              = -20
  res@mpMaxLatF              = 20
        res@mpMinLonF              =40
  res@mpMaxLonF              =110
   res@mpCenterLonF           = 180     ; This is necessary to get the
   
  res@gsnAddCyclic =False
  res@cnFillOn             = True         ; turn on color fill
  res@cnLinesOn            = False        ; True is default
;res@cnLineLabelsOn       = False        ; True is default
  res@cnFillPalette        = "BlWhRe"     ; set color map
  res@lbLabelBarOn         = False  
   res@cnLineLabelsOn      = False      ; turn off individual lb's
   res@cnInfoLabelOn=False
                                          ; set symmetric plot min/max
  symMinMaxPlt(eof, 16, False, res)       ; contributed.ncl
; panel plot only resources
  resP                     = True         ; modify the panel plot
  resP@gsnMaximize         = True         ; large format
  resP@gsnPanelLabelBar    = True         ; add common colorbar
  resP@lbLabelAutoStride   = True         ; auto stride on labels
  yStrt                    = 1982
  yLast                    = 2015
  resP@txString            = "SST: "+": "+yStrt+"-"+yLast
;*******************************************
; first plot
;*******************************************
  do n=0,neof-1
     res@gsnLeftString  = "EOF "+(n+1)
     res@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
     plot(n)=gsn_csm_contour_map_ce(wks,eof(n,:,:),res)
  end do
  gsn_panel(wks,plot,(/neof,1/),resP)     ; now draw as one plot
;*******************************************
; second plot
;*******************************************
; EOF time series  [bar form]
  rts           = True
  rts@gsnDraw   = False       ; don't draw yet
  rts@gsnFrame  = False       ; don't advance frame yet
  rts@gsnScale  = True        ; force text scaling               
; these four rtsources allow the user to stretch the plot size, and
; decide exactly where on the page to draw it.
  rts@vpHeightF = 0.40        ; Changes the aspect ratio
  rts@vpWidthF  = 0.85
  rts@vpXF      = 0.10        ; change start locations
  rts@vpYF      = 0.75        ; the plot

  rts@tiYAxisString = ""                    ; y-axis label      
  rts@gsnYRefLine           = 0.              ; reference line   
  rts@gsnXYBarChart         = True            ; create bar chart
  rts@gsnAboveYRefLineColor = "red"           ; above ref line fill red
  rts@gsnBelowYRefLineColor = "blue"          ; below ref line fill blue
; panel plot only resources
  rtsP                      = True            ; modify the panel plot
  rtsP@gsnMaximize          = True            ; large format
  rtsP@txString             = "SST: "+": "+yStrt+"-"+yLast
year=ispan(1,34,1)
  print(eof_ts)
; create individual plots
  do n=0,neof-1
     rts@gsnLeftString  = "EOF "+(n+1)
     rts@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
     plot(n) = gsn_csm_xy (wks,year,eof_ts(n,:),rts)
  end do
  gsn_panel(wks,plot,(/neof,1/),rtsP)     ; now draw as one plot
end


使用ssta数据：

结果：数值量值太小，请教为什么出现这种状况？ （与书中数值量级不一样）

songwei

用特征向量乘以时间系数得到的量级和参考文献中才是一样的，其实这里你算的是标准化以后进行了EOF，可以尝试将EOF得到的时间系数标准化以后回归到原始场，这样量级就和参考文献中一样了。

songwei

我自己以前写的 你参考一下

qq469015280

songwei 发表于 2017-12-11 09:12
用特征向量乘以时间系数得到的量级和参考文献中才是一样的，其实这里你算的是标准化以后进行了EOF，可以尝 ...

可是参考文献中 也是对sst异常做的EOF的结果，并不是回归后的结果。

songwei

qq469015280 发表于 2017-12-11 14:55
可是参考文献中 也是对sst异常做的EOF的结果，并不是回归后的结果。

得到的时间系数标准化后再对原始场做回归得到的结果和用距平做EOF得到的空间模态是一致的。
如果你不信可以自己试试喽