QBO位相划分作图问题求助

hanxiang0721

我最近刚开始用ncl绘图，按照官网上的MJO划分位相的例子，自己写了一个脚本，但是画出来的图不理想，想问一下有没有人知道哪里错了。
下面是脚本
;************************************************************************
;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
        neof=2
        latS=-20
        latN=30
        lonL=60
        lonR=180
        ymdStrt=19790101
        ymdLast=20131231
        yrStrt=ymdStrt/10000
        yrLast=ymdLast/10000
;*************绘图路径****************************************************
        pltDir="./"
        pltType="png"
        pltName="MJO_Phase"
;*************资料路径*****************************************************
        diro="/home/hanxiang/dachuang/"
        filo="olr.day.mean.nc"
;*************设置lanczos滤波器********************************************
        ihp=2       ;带通滤波
        nWgt=201
        sigma=1.0
        fca=1./20       ;滤波上限
        fcb=1./10       ;滤波下限
        wgt=filwgts_lanczos(nWgt,ihp,fca,fcb,sigma)
;*************读入资料******************************************************
        f=addfile(diro+filo,"r")
        TIME=f->time
        YMD=cd_calendar(TIME,0)
        iyear=floattoint(YMD(:,0))
        imon=floattoint(YMD(:,1))
        nt=ind(iyear.ge.1979.and.iyear.le.2013.and.imon.ge.5.and.imon.le.10)  ;挑选1979年至2013年5~10月的资料
        OLR=f->olr(nt,{latS:latN},{lonL:lonR})
        time=f->time(nt)
        dimx=dimsizes(OLR)
        ntim=dimx(0)
        nlat=dimx(1)
        nlon=dimx(2)
;************滤波**************************************************************
        olr_filter_1=wgt_runave_Wrap(OLR(lat|:,lon|:,time|:),wgt,0)
        olr_filter=dim_rmvmean_Wrap(olr_filter_1)
;************EOF分解**********************************************************
        eof=eofunc_Wrap(olr_filter,neof,False)
        eof_ts=eofunc_ts_Wrap(olr_filter,eof,False)      ;时间序列
        pc1=eof_ts(0,:)
        pc2=eof_ts(1,:)
;        index=pc1^2+pc2^2
;*************位相划分************************************************************
        nPhase=8
        angBnd=new((/2,nPhase/),"float")
        angBnd(0,:)=fspan(0,315,nPhase)
        angBnd(1,:)=fspan(45,360,nPhase)
        r2d=180./(4.*atan(1.0))
        ang=atan2(pc1,pc2)*r2d
        nn=ind(ang.lt.0)
        ang(nn)=ang(nn)+360
        nDays=new(nPhase,"integer")
        pLabel="P"+ispan(1,nPhase,1)+": "
;***************绘图要素设置********************************************************
        pltPath=pltDir+pltName
        wks=gsn_open_wks(pltType,pltPath)
        gsn_define_colormap(wks,"ViBlGrWhYeOrRe")
        plot=new(nPhase,graphic)
        res=True
        res@gsnDraw=False
        res@gsnFrame=False
        res@gsnSpreadColors=True
        res@gsnAddCyclic=False
        ;res@gsnSpreadColors= True  
        res@mpFillOn=False
;**************设置绘图经纬度范围******************************************************
        res@mpMinLatF=latS
        res@mpMaxLatF=latN
        res@mpMinLonF=lonL
        res@mpMaxLonF=lonR
        res@mpCenterLonF=120
;******************************************************************************************
        res@cnFillOn=True
        res@cnLinesOn=False
        res@cnLineLabelsOn=False
        res@lbLabelBarOn=False
        symMinMaxPlt(eof,32,False,res)
        resP=True
        resP@gsnMaximize=True
        resP@gsnPanelLabelBar=True
        resP@lbLabelAutoStride=True
        resP@gsnMaximize= True  
;        resP@lbLabelFontHeightF=0.0125
        res@gsnLeftString=""
        res@gsnRightString=""
        resP@txString=yrStrt+"-"+yrLast+": May to Oct"
        do n=0,nPhase-1
                na=n+4
                if(na.gt.7) then
                        na=na-8
                end if
          nt1= ind(ang.ge.angBnd(0,n) .and. ang.lt.angBnd(1,n))       
          if (.not.all(ismissing(nt1))) then
         xAvg = dim_avg_Wrap(olr_filter(lat|:,lon|:,time|nt1) )
         nDays(na) = dimsizes(nt1)
           res@tmXBLabelsOn         = False   
         res@tmXBOn               = False  
             res@tmXBLabelsOn     = True     
             res@tmXBOn           = True   
         end if
          plot(na) = gsn_csm_contour_map_ce(wks,xAvg,res)
         delete(nt1)               
      resP@gsnPanelFigureStrings= pLabel+nDays
  gsn_panel(wks,plot,(/nPhase,2/),resP)     
end do
end

绘出的图是这个样子的

lllllys

请问楼主关于QBO位相划分作图的问题解决了嘛想请教一下

wsx

lllllys 发表于 2021-1-10 23:00
请问楼主关于QBO位相划分作图的问题解决了嘛想请教一下

同问，问题解决了吗？