ncl 画全球slp的eof（官网学习的例子）

liyuan3970

官网上给的eof分析的例子我试着画一下，发现空间模态和时间序列是一致的，但是方差贡献不一样不知道是为什么？经纬度和时间截取的也是一致的，

liyuan3970

begin
f1=addfile("HadISST_sst.nc","r")

sst=f1->sst(:,:,:)
sst@missing_value=-1e+30
y=month_to_annual(sst, 1)
b=y(88:140,:,:)
copy_VarCoords(sst(0,:,:),b(0,:,:))
b!0="time"
time=ispan(1958,2010,1)
b&time=time
b&time@long_name="year"
b&time@units="year"
b@missing_value=-1e+30
;x=b
x=runave_n_Wrap(b, 3, 0, 0)
copy_VarCoords(b,x)
printVarSummary(x)
optEOF = False
optETS = False      
neval  =2
eof= eofunc_Wrap(x({latitude|:},{longitude|:},{time|:}),neval,optEOF)
eof_ts=eofunc_ts_Wrap (x({latitude|:},{longitude|:},{time|:}),eof, optETS)
eof_ts=eof_ts/10000
eof=eof*10000
;f1->ts=eof_ts
;print(max(eof(0,:,:)))
;print(max(eof(1,:,:)))
;print(min(eof(0,:,:)))
;print(min(eof(1,:,:)))
wks = gsn_open_wks("png","1950_2014hwf_std_eof(1)")
res2=True
res2@gsnDraw = False
res2@gsnFrame = False
res2@gsnAddCyclic        = False
res2@mpCenterLonF         = 100  
res2@mpLimitMode       = "LatLon"
res2@mpMinLatF         = -40         
res2@mpMaxLatF         = 60
res2@mpMinLonF         = -180
res2@mpMaxLonF         = 180
res2@cnFillOn       = True
;res2@cnFillPattern ="Explicit"
res2@cnLevelSelectionMode="ManualLevels"
res2@cnMinLevelValF=-0.6
res2@cnMaxLevelValF=0.6
res2@cnLevelSpacingF=0.1
;res2@cnFillColors=(/129,145,161,177,193,200,215,225/)
res2@mpDataBaseVersion="Ncarg4_1"       ;ÖDμè·Ö±æÂê
res2@mpDataSetName="Earth..4"           ;μú4°æμØí¼￡¬óDÖD1ú±ß½çêy¾Y
res2@mpOutlineOn            = True
res2@mpOutlineSpecifiers=(/"China:states","Taiwan"/)   ;ÖD1ú±ß½çóDÎêìa￡¬è±2ØÄÏ¡￠ì¨íå
res2@mpOutlineBoundarySets ="NoBoundaries"

res2@cnSmoothingOn      =True
res2@cnRasterSmoothingOn= True
;gsn_define_colormap(wks,"leng")
;res2@tiMainString    = "hwfeof(0)1950-2014dtrend"
;res2@tmXBMode = "Explicit"
;res2@tmXBValues = (/"70","80","90","100","110","120","130","140"/)
;res2@tmXBLabels  = (/"70E","80E", "90E", "100E","110E","120E","130E","140E"/)

res2@cnMissingValFillColor="white"


rts2           = True
rts2@gsnDraw = False
rts2@gsnFrame = False
rts2@vpHeightF = 0.40        ; Changes the aspect ratio
rts2@vpWidthF  = 0.80
rts2@vpXF      = 0.80        ; change start locations
rts2@vpYF      = 0.75        ; the plot
rts2@tiYAxisString = ""                  
rts2@tiXAxisString = "year"
rts2P                      = True            ; modify the panel plot
rts2P@gsnMaximize          = True            ; large format
;rts2@trYMaxF=2
;rts2@trYMinF=-1.2
  
rts2@gsnYRefLine           = 0.              ; reference line   
rts2@gsnXYBarChart         = True            ; create bar chart
rts2@gsnAboveYRefLineColor = "red"           ; above ref line fill red
rts2@gsnBelowYRefLineColor = "blue"          ; below ref line fill blue

rts2@tmXBMode = "Explicit"
rts2@tmXBValues=(/"1950","1960","1970","1980","1990","2000","2010"/)
  rts2@tmXBLabels=(/"1950","1960","1970","1980","1990","2000","2010"/)
;rts2@tiMainString    = "AP_Ts_02:1950-2014"
rts2@gsnLeftString  = "EOF "+(0+1)
;rts2@gsnRightString = sprintf("%5.1f", eof@pcvar(1)) +"%"


plots = new(2,graphic)

  res2@tiMainString = "1950_2014hwf-std-eof1"
  res2@gsnLeftString  = "EOF1 "
  ;cmap2 = read_colormap_file("14")
;res2@cnFillColors = cmap2
  ;gsn_define_colormap(wks,"14")
res2@gsnRightString = sprintf("%5.1f", eof@pcvar(0)) +"%"
  plots(0) = gsn_csm_contour_map(wks,eof(0,:,:),res2)



plots(1)= gsn_csm_xy (wks,time,eof_ts(0,:),rts2)   






  pres = True
  ;pres@tiMainString = "1950_2014hwf_eof"
  pres@gsnMaximize = True
  pres@gsnPanelDebug=True
  pres@gsnPanelRowSpec = True
  pres@gsnPanelLabelBar = False
  pres@lbBoxLinesOn = False
  gsn_panel(wks,plots,(/1,2/),pres)
end
我自己的

liyuan3970

; ==============================================================
; eof_1.ncl
;
; Concepts illustrated:
;   - Calculating EOFs
;   - Drawing a time series plot
;   - Using coordinate subscripting to read a specified geographical region
;   - Rearranging longitude data to span -180 to 180
;   - Calculating symmetric contour intervals
;   - Drawing filled bars above and below a given reference line
;   - Drawing subtitles at the top of a plot
;   - Reordering an array
; ==============================================================
; NCL V6.4.0 has new functions eofunc_n_Wrap and
; eofunc_ts_n_Wrap that allow you to calculate the EOFs without
; first having to first reorder the data. See eof_1_640.ncl.
; ==============================================================
; Calculate EOFs of the Sea Level Pressure over the North Atlantic.
; ==============================================================
; The slp.mon.mean file can be downloaded from:
; http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.surface.html
; ==============================================================
; These files are loaded by default in NCL V6.2.0 and newer
; 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
; ==============================================================
; User defined parameters that specify region of globe and
; ==============================================================
  latS   =  25.
  latN   =  80.
  lonL   = -70.
  lonR   =  40.

  yrStrt = 1979
  yrLast = 2003

  season = "DJF"    ; choose Dec-Jan-Feb seasonal mean

  neof   = 3        ; number of EOFs
  optEOF = True      
  optEOF@jopt = 0   ; This is the default; most commonly used; no need to specify.
;;optEOF@jopt = 1   ; **only** if the correlation EOF is desired

  optETS = False

; ==============================================================
; Open the file: Read only the user specified period
; ==============================================================
  f      = addfile ("slp.mon.mean.nc", "r")

  TIME   = f->time
  YYYY   = cd_calendar(TIME,-1)/100                 ; entire file
  iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)

  slp    = f->slp(iYYYY,:,:)
  printVarSummary(slp)                              ; variable overview

; ==============================================================
; dataset longitudes span 0=>357.5
; Because EOFs of the North Atlantic Oscillation are desired
; use the "lonFlip" (contributed.ncl) to reorder
; longitudes to span -180 to 177.5: facilitate coordinate subscripting
; ==============================================================
  slp    = lonFlip( slp )
  printVarSummary(slp)                              ; note the longitude coord

; ==============================================================
; compute desired global seasonal mean: month_to_season (contributed.ncl)
; ==============================================================
  SLP    = month_to_season (slp, season)
  nyrs   = dimsizes(SLP&time)
  printVarSummary(SLP)

; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
; =================================================================
  rad    = 4.*atan(1.)/180.
  clat   = f->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|:)
  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  
;;yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here

;============================================================
; PLOTS
;============================================================
  wks = gsn_open_wks("png","eof")         ; send graphics to PNG file
  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             = False        ; turn off map fill
  res@mpMinLatF            = latS         ; zoom in on map
  res@mpMaxLatF            = latN
  res@mpMinLonF            = lonL
  res@mpMaxLonF            = lonR

  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        ; turn off individual lb's

                                          ; 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

  yStrt                    = yyyymm(0)/100
  yLast                    = yyyymm(nyrs-1)/100
  resP@txString            = "SLP: "+season+": "+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 = "Pa"                    ; 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             = "SLP: "+season+": "+yStrt+"-"+yLast

  year = yyyymm/100

; 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

liyuan3970

上面是官网的脚本

liyuan3970

; =================================================================
; create weights:  sqrt(cos(lat))   [or sqrt(gw) ]
; =================================================================
  rad    = 4.*atan(1.)/180.
  clat   = f->lat           
  clat   = sqrt( cos(rad*clat) )                 ; gw for gaussian grid
官网上这个地方我没有看懂，查了atan是弧度单位也不知道为什么要这么计算

liyuan3970

。。。。。。自己画的eof贴错了，sorry，但是基本上是这样写的就变了下数据而已

明德格物

学习学习了

柚子皮

gw for gaussian grid，似乎是求了个纬向权重。如果整个程序只改了数据的话，应该是不会有问题的，可能是不同资料的出入，整体的模态一致就应该没问题了。

liyuan3970

柚子皮 发表于 2017-7-6 20:15
gw for gaussian grid，似乎是求了个纬向权重。如果整个程序只改了数据的话，应该是不会有问题的，可能是不 ...

我也觉得只要模态一致就不会有问题了，谢谢您，我再看看这个纬向权重是什么意思

东东的小画家

感谢分享，这个好棒找了好久