HADISST 海温数据 EOF分解 北太平洋，奇葩图

guizi

最近两天在用  ncl  对北太平洋 PDO 进行 EOF 分解 ，
几乎相同的程序，对 ERSST资料和 HADISST资料进行 EOF 分解，结果差别很大，

我的处理 方法是，先去个趋势，然后做个距平，最后进行EOF分解
下面是 ERSST 的 资料结果，正常


除了 数据读取的头脚本不一样，其他计算和画图脚本完全一样，如下是 HADISST的 结果
数据来自 http://www.metoffice.gov.uk/hadobs/hadisst/data/download.html




HADISST的 时间系数也比较怪，都 10**30 次方了？？？？
求解释：
如下是 ncl脚本比较多 我就单独占个 沙发吧

Aires

HadISST的数据我没用过，你给问下你师兄了。我刚用ERSST得到的PDO模态是这个样子的。

guizi

;************************************************
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl"
;************************************************
begin
;*************  set up  some parameter        *********   
  lonL   = 100
  lonR   = 280

  latB   = -10
  latT   = 60

  yrStrt = 1950
  yrLast = 2013
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  neof = 8
;************  open files and  read the data ***********  
  f     = addfile("../HadISST_sst.nc","r")
  time   = f->time                   ; days since 0000-09-01 00:00:00
  YYYY   = cd_calendar(time,-1)/100                 ; entire file
  iYYYY  = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
;  sst  = f->sst(iYYYY,{latT:latB},:)
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  g_sst  = f->sst(iYYYY,:,:)
  g1_sst = lonPivot(g_sst,0.5)
  sst    = g1_sst(:,{latB:latT},{lonL:lonR})
;*************  deal with the data     *****

  sst    = where(abs(sst).gt.500,sst@_FillValue,sst)

  printVarSummary(sst)
  printMinMax(sst,True)
  
  f_sst   = sst

  printVarSummary(f_sst)
  printMinMax(f_sst,True)

  ntime   = dimsizes(sst&time)
  ny      = dimsizes(sst&latitude)
  nx      = dimsizes(sst&longitude)
;******************************************************


  copy_VarCoords(sst,f_sst)
;  printVarSummary(f_sst)

;  y_sst = dtrend_msg_n(f_sst&time,f_sst,True,False,0)
;  copy_VarCoords(sst,y_sst)

  y_sst = dtrend_leftdim(f_sst,False)
;  printVarSummary(y_sst)
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   temp   = y_sst
   nt     = ntime/12
   print(nt)
   month  = 12
   temp_mon  = new((/month,nt,ny,nx/),float)
   temp_d    = new((/ntime,ny,nx/),float)

   do imon=0,month-1
      temp_mon(imon,:,:,:) =  dim_rmvmean_n_Wrap(temp(imon:ntime-1:12,:,:),0)
   end do

   do it=0,nt-1
   do imon=0,month-1
      temp_d(it*12+imon,:,:)  =  temp_mon(imon,it,:,:)
   end do
   end do
      copy_VarCoords(temp,temp_d)


;~~~~~~~~~~~~~~~~~~~~~~~~~~~
  x    = temp_d
  printVarSummary(x)

;*******************************************
;  EOF
;*******************************************
  X = x(latitude|:,longitude|:,time|:)              ; Space x Time

  optEof = True
  eof    = eofunc_Wrap( X, neof, optEof)
  eof_ts = eofunc_ts_Wrap( X, eof, False)

  printVarSummary(eof)
  printVarSummary(eof_ts)

  yyyymm = cd_calendar(eof_ts&time,-2)/100
  year   = yyyymm/100


;*******************************************
;  plots
;*******************************************
  wks = gsn_open_wks("pdf","pdo_eof")
  gsn_define_colormap(wks,"posneg_1")     ; choose colormap
;
; This will not be necessary in V6.1.0 and later. Named colors can
; be used without having to first add them to the color map.
;
  dum = NhlNewColor(wks,0.7,0.7,0.7)      ; add gray to color table

  plot = new(neof,graphic)                ; create graphic array
                                          ; only needed if paneling
  res                      = True
  res@gsnDraw              = False        ; don't draw yet
  res@gsnFrame             = False        ; don't advance frame yet
  res@gsnSpreadColors      = True         ; spread out color table
  res@gsnSpreadColorEnd    = -2           ; don't use added gray
  res@gsnAddCyclic         = False        ; data not cyclic


  res@mpCenterLonF         = 190.         ; defailt is 0 [GM]
  res@mpMinLatF            = min(x&latitude)
  res@mpMaxLatF            = max(x&latitude)
  res@mpMinLonF            = min(x&longitude)
  res@mpMaxLonF            = max(x&longitude)
  res@mpFillDrawOrder      = "PostDraw"

  res@cnFillOn             = True         ; turn on color fill
  res@cnLinesOn            = True         ; True is default
  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
  resP@lbLabelAutoStride   = True         ; auto stride on labels


  resP@txString            = "Kaplan_sst (Pierce)"
  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(0:3),(/2,2/),resP)     ; draw all 'neof' as one plot

;*******************************************
; time series (principal component) plot
;*******************************************
  eof_ts@long_name = "Amplitude"

  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 resources 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@gsnYRefLine           = 0.              ; reference line
  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

;  year = cd_calendar(eof_ts&time,1)


  resP@txString            = "KAplan_sst: Pierce"
  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(0:3),(/2,2/),rtsP)        ; draw all 'neof' as one plot

end


                                       

Aires

检查你HadISST原始SST数据是否是正确的，另外，PDO是对20N以北做EOF分解，你的区域是不是大了？你用ERSST分解的EOF1也不太像是PDO啊，PC1看不出来是decadal oscillation

guizi

Aires 发表于 2014-3-10 09:44

                               
登录/注册后可看大图

检查你HadISST原始SST数据是否是正确的，另外，PDO是对20N以北做EOF分解，你的区域是不是大了？你用ERSST分 ...

                               
登录/注册后可看大图

http://jisao.washington.edu/pdo/
上面是  google 的 pdo 差不多呀

guizi

guizi 发表于 2014-3-10 10:01

                               
登录/注册后可看大图

http://jisao.washington.edu/pdo/
上面是  google 的 pdo 差不多呀

google 的

如下是 脚本读取的 ERRSST和 HADISST的 相对应 的 海温图，也就是说数据读取的结果，是正确的......??????

   hadisst    ERSST


hadisst ERSST

Aires

guizi 发表于 2014-3-10 10:01

                               
登录/注册后可看大图

http://jisao.washington.edu/pdo/
上面是  google 的 pdo 差不多呀

但是你的PC1也不像是年代际振荡啊

guizi

Aires 发表于 2014-3-10 10:15

                               
登录/注册后可看大图

但是你的PC1也不像是年代际振荡啊

画图的时候 平滑一下，是不是就可以了，时间系数不是很像，模态很像， 可是 hadisst 很奇葩呀

Aires

guizi 发表于 2014-3-10 10:51

                               
登录/注册后可看大图

画图的时候 平滑一下，是不是就可以了，时间系数不是很像，模态很像， 可是 hadisst 很奇葩呀

你再检查一下你程序，或者用20度以北去做做看，我一会给你看下我用ERSST做的结果。

guizi

Aires 发表于 2014-3-10 11:08

                               
登录/注册后可看大图

你再检查一下你程序，或者用20度以北去做做看，我一会给你看下我用ERSST做的结果。

ERSST是 正常的 ，hadisst 有点不正常呢，我刚刚问了师兄，那个 ERSST的结果，第一模态叫 类enso模态，第二模态叫 北太平洋模态，第三模态才是 pdo，因为我去的范围比较大，

hadisst的 资料比较 难搞，结果莫名其妙，

我看hadisst的 数据，里面是不是有海冰数据呢，我打印最大最小值的时候，有-1000 的值，但是缺测是 1×10**30 方，然后我用了 sst    = where(abs(sst).gt.500,sst@_FillValue,sst)  但是结果图，没变

guizi

Aires 发表于 2014-3-10 11:08

                               
登录/注册后可看大图

你再检查一下你程序，或者用20度以北去做做看，我一会给你看下我用ERSST做的结果。

http://www.metoffice.gov.uk/hadobs/hadisst/data/download.html

nc 格式的那个，，210m的