爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 2012|回复: 1

[求助] NCL制作mveof的图代码报错,求帮忙

[复制链接]

新浪微博达人勋

发表于 2021-11-10 16:50:26 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
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   =  -40.
  latN   =   40.
  lonL   =  0.
  lonR   =  180.
  neof=2
  ymdStrt = 19810101                         ; start yyyymmdd
  ymdLast = 20101231                         ; last  
  yrStrt  = ymdStrt/10000
  yrLast  = ymdLast/10000

  pltDir  = "./"                             ; plot directory
  pltType = "ps"
  pltName = "monsoon_mveof"   

  diri    = "./"                             ; input directory

  filolr  = "olr.day.anomalies.1981-2010.nc"
  ;filu850 = "uwnd.day.850.anomalies.1981-2010.nc"
  ;filv850 = "vwnd.day.850.anomalies.1981-2010.nc"

;************************************************
; create BandPass Filter
;************************************************
; ihp      = 2                             ; bpf=>band pass filter
;nWgt     = 201
;sigma    = 1.0                           ; Lanczos sigma
;fca      = 1./100.
;fcb      = 1./20.
;wgt      = filwgts_lanczos (nWgt, ihp, fca, fcb, sigma )

;***********************************************************
; Find the indices corresponding to the start/end times
;***********************************************************
   f       = addfile (diri+filolr , "r")                        
   TIME    = f->time                          ; days since ...

   YMD     = cd_calendar(TIME, -2)            ; entire (time,6)

   iStrt   = ind(YMD.eq.ymdStrt)              ; index start
   iLast   = ind(YMD.eq.ymdLast)              ; index last
   delete([/ TIME, YMD /])
;***********************************************************
; Read anomalies
;***********************************************************
   ;f =addfile("olr.day.anomalies.1981-2010.nc","r")
   fu=addfile("uwnd.day.850.anomalies.1981-2010.nc","r")
   fv=addfile("vwnd.day.850.anomalies.1981-2010.nc","r")
   TIME    = f->time
   olr=f->olr({iStrt:iLast},{latS:latN},{lonL:lonR})
   olr=dim_standardize_n_Wrap(olr, 0, 0)

   delete(TIME)

   TIME    = fu->time
   uwnd=short2flt(fu->uwnd({iStrt:iLast},{latS:latN},{lonL:LonR}))
   u=dim_standardize_n_Wrap(uwnd, 0, 0)

   delete(TIME)

   TIME    = fv->time
   vwnd=short2flt(fv->vwnd({iStrt:iLast},{latS:latN},{lonL:LonR}))
   v=dim_standardize_n_Wrap(vwnd, 0, 0)

   delete(TIME)

;*******************Write into one dataset*****************************
;               Combine the normalized data into one variable                                                       ;
;**********************************************************************
  dimw    = dimsizes(olr)
  mtim    = dimw(0)
  mlat    = dimw(1)
  mlon    = dimw(2)
  cdata   = new ( (/3*mlat,3*mlon,mtim/), typeof(olr), getFillValue(olr))
  do m1=0,mlat-1
    do m2=0,mlon-1
        cdata(m1,       m2,       :)=(/olr(:,m1,m2)/)
        copy_VarMeta(olr,cdata(m1,       m2,       :) )
        cdata(m1+mlat,  m2+mlon,  :)=(/u     (:,m1,m2)/)
        copy_VarMeta(olr, cdata(m1+mlat,  m2+mlon,  :))
        cdata(m1+2*mlat,m2+2*mlon,:)=(/v     (:,m1,m2)/)
        copy_VarMeta(olr, cdata(m1+2*mlat,m2+2*mlon,:))
    end do
  end do
;***************************************t*****************************
;               Compute **combined** EOF; Sign of EOF is arbitrary                                                      ;
;**********************************************************************
eof_cdata    = eofunc(cdata   , neof, False)
eof_ts_cdata = eofunc_ts(cdata,eof_cdata,False)
printVarSummary(eof_cdata)
printVarSummary(eof_ts_cdata)
nvar=3                                 ;olr,u,v
ceof = new( (/nvar,neof,mlat,mlon/), typeof(cdata), getFillValue(cdata))
do n=0,neof-1
     ceof(0,n,:,:) = eof_cdata(n,0:mlat-1,0:mlon-1)           ; olr
     ceof(1,n,:,:) = eof_cdata(n,mlat:mlat*2-1,mlon:2*mlon-1) ; u
     ceof(2,n,:,:) = eof_cdata(n,2*mlat:,2*mlon:)             ; v
  end do
ceof!0   = "var"
ceof!1   = "eof"
ceof!2   = "lat"
ceof!3   = "lon"   
ceof&lat = olr&lat
ceof&lon = olr&lon
printVarSummary(ceof)
printMinMax(ceof(0,1,:,:),False)

;***************************************t*****************************
;               plot now                                              ;
;**********************************************************************
wks=gsn_open_wks("ps","MVEOF")
plot = new(neof,graphic)
plot2= new(neof,graphic)
gsn_define_colormap(wks,"MPL_RdYlGn")
  res                      = True         
  res@gsnDraw              = False        ; don't draw yet
  res@gsnFrame             = False        ; don't advance frame yet
; res@gsnPolar             = "SH"
  res@gsnAddCyclic         = False
  res@gsnMaximize          = True
  res@mpFillOn             = True        ; turn off map fill
  res@mpOutlineOn          = False
  res@mpMinLatF            = -40
  res@mpMaxLatF            = 40
  res@mpMaxLonF            = 180
  res@mpMinLonF            = 0

  res@cnFillOn             = True         ; turn on color fill
  ;res@cnFillPalette        = "BlueWhiteOrangeRed"  
  res@cnLinesOn            = False        ; True is default
  res@cnLineLabelsOn       = False        ; True is default
  res@lbLabelBarOn         = False        ; turn off individual lb's
  ;res@cnLevelSelectionMode = "ManualLevels"
  ;res@cnMaxLevelValF       = -0.02
  ;res@cnMinLevelValF       = 0.02
  ;res@cnLevelSpacingF      = 0.005
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/-0.02,-0.015,-0.01,-0.005,0,0.005,0.01,0.015,0.02/)
  res@cnFillColors         = (/5,25,38,56,0,0,80,94,108,121/)
; panel plot only resources

  vcres                    = True
  vcres@gsnDraw            = False
  vcres@gsnFrame           = False
  vcres@vcRefLengthF       =0.045
  vcres@vcMinDistanceF     =0.017
  ;vcres@vcGlyphStyle     = "CurlyVector"

  resP                     = True         ; modify the panel plot
  resP@gsnMaximize         = True         ; large format
  resP@gsnPanelLabelBar    = True         ; add common colorbar
  resP@gsnPaperOrientation = "portrait"   ; force portrait

  resP@txString            = ""

;*******************************************
; first plot
;*******************************************
  do n=0,neof-1
     res@tiMainString   = "MVEOF "+(n+1)
     res@gsnLeftString  = "olr and 850hPa winds"
     res@gsnRightString = ""
     plot(n) =gsn_csm_contour_map(wks,ceof(0,n,:,:),res)
     plot2(n)=gsn_csm_vector     (wks,ceof(1,n,:,:), ceof(2,n,:,:), vcres)
     overlay(plot(n),plot2(n))
     draw(plot(n))
     frame(wks)
  end do
  gsn_panel(wks,plot,(/neof,1/),resP)     ; now draw as one plot
end





下载.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2023-1-31 20:48:22 | 显示全部楼层
求问你的报错修改完了吗,能不能请教几个问题
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表