Re: MJO clivar

You can try debugging yourself by printing values *prior*
to the statement in question. If the following is in error

print(xBegin)
print(yBegin)
print(pc)

Maybe one of them is _FillValue?

plot@$unique_string("dum")$ = gsn_add_polyline(wks, \
       plot, (/xBegin,pc1(nt)/), (/yBegin,pc2(nt)/), plLine))

There is no "ismissing" in the driver script.

Good luck

On 5/31/12 5:27 PM, juki juki wrote:
> Thanks for nice suggestion Dennis, it works now.
>
> However, allow me to ask another thing in plotting MJOclivar.
>
> Question-1: I try to run mjoclivar_15.ncl, however I got the following
> error:
>
> fatal:The result of the conditional expression yields a missing
> value. NCL can not determine branch, see ismissing function
>
> fatal:Execute: Error occurred at or near line 115 in mjoclivar_15.ncl
> (the following line: plot@$unique_string("dum")$ = gsn_add_polyline(wks,
> plot, (/xBegin,pc1(nt)/), (/yBegin,pc2(nt)/), plLine))
>
> I didn't change much the available code, only the input data and
> start-ending time. The complete code is:
>
> --
> 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/diagnostics_cam.ncl"
>
> begin
> ymdStrt = 20110101 ; start yyyymmdd
> ymdLast = 20111231 ; last
>
> pltDir = "./" ; plot directory
> pltType = "ps"
> pltName = "mjoclivar15" ; yrStrt+"_"+yrLast
>
> pltMovie= False ; animation
> pltTitle= "MJO Phase: 15S-15N: "+ymdStrt+"-"+ymdLast
> ;***********************************************************
> ; Open PC components file created in 'mjo_14.ncl'
> ;***********************************************************
> diri = "e:/RData/MJO_cal/Anom/" ; input directory
> fMJO = "MJO_PC_INDEX.nc" ; created in mjo_14.ncl
> f = addfile (diri+fMJO, "r")
>
> ;***********************************************************
> ; Find the indices corresponding to the start/end times
> ;***********************************************************
> 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)
> delete(YMD )
>
> ;***********************************************************
> ; Read the data for the desired period
> ;***********************************************************
> pc1 = f->PC1(iStrt:iLast)
> pc2 = f->PC2(iStrt:iLast)
> mjo_indx= f->MJO_INDEX(iStrt:iLast)
>
> time = f->time(iStrt:iLast)
> ymdhms = cd_calendar(time, 0)
>
> ntim = dimsizes( time )
> ; print(ymdhms)
>
> ;------------------------------------------------------------
> ; PLOT
> ;------------------------------------------------------------
>
> opt = True ; Background options
> ;opt@gsnMaximize = True ; make large
> opt@tiMainString = pltTitle
>
> plMark = True ; poly marker
> plMark@gsMarkerIndex = 16 ; solid circle
> plMark@gsMarkerSizeF = 0.005
> plMark@gsMarkerThicknessF = 1
>
> plLine = True ; line segments
> plLine@gsLineThicknessF = 1.0 ; 1.0 is default
>
> txres = True
> txres@txFontHeightF = 0.010 ; size of font for printed day
>
> namMonth = (/ "Jan","Feb","Mar","Apr","May","June" \
> ,"July","Aug","Sep","Oct","Nov","Dec" /)
>
> colMonth = (/2,3,5,6,7,8,10,11,12,13,18,20/) ; indices into color table
>
> imon = floattoint( ymdhms(:,1) ) ; convenience
> iday = floattoint( ymdhms(:,2) ) ; sunscripts must be integer
> print(colMonth)
>
> pltPath = pltDir+pltName
> if (.not.pltMovie) then
> wks = gsn_open_wks(pltType, pltPath) ; open workstation
> gsn_define_colormap(wks,"radar_1")
> end if
>
> plot = mjo_phase_background(wks, opt) ; generic phase space background
> xBegin= pc1(0) ; need for start of the line
> yBegin= pc2(0)
> plMark@gsMarkerColor = "black" ; indicate initial point
> plMark@gsMarkerSizeF = 2.5*plMark@gsMarkerSizeF ; make larger
> plot@$unique_string("dum")$ = gsn_add_polymarker(wks, plot, xBegin,
> yBegin, plMark)
> plMark@gsMarkerSizeF = plMark@gsMarkerSizeF/2.5 ; reset
>
> nMon = 0
> monInfo = new ((/12,2/),"string")
> monInfo(nMon,0) = namMonth(imon(0)-1)
> monInfo(nMon,1) = colMonth(imon(0)-1)
> label_opt = 5 ; label every # day 0 = don't label days
> label_color = True ; draw month label with its proper color
> do nt=1,ntim-1
> if (pltMovie) then
> ext = "_"+sprinti("%4.0i", nt )
> wks = gsn_open_wks(pltType,pltName+ext)
> gsn_define_colormap(wks,"radar_1")
> end if
> ; color changes w month
> plLine@gsLineColor = colMonth(imon(nt)-1) ; -1 is cuz NCL is 0-based
> plot@$unique_string("dum")$ = gsn_add_polyline(wks, plot,
> (/xBegin,pc1(nt)/), (/yBegin,pc2(nt)/), plLine)
> if (label_opt.eq.0) then
> plMark@gsMarkerColor = plLine@gsLineColor ; same as line color
> plot@$unique_string("dum")$ = gsn_add_polymarker(wks, plot, xBegin,
> yBegin, plMark)
> else
> if (iday(nt)%label_opt.eq.0) then
> txres@txFontColor = "black"
> plot@$unique_string("dum")$ = gsn_add_text(wks, plot, iday(nt)+"",
> xBegin, yBegin, txres)
> else ; marker only
> plMark@gsMarkerColor = plLine@gsLineColor ; same as line color
> plot@$unique_string("dum")$ = gsn_add_polymarker(wks, plot, xBegin,
> yBegin, plMark)
> end if
> end if
>
> if (iday(nt).eq.1) then
> nMon = nMon+1
> monInfo(nMon,0) = namMonth(imon(nt)-1)
> monInfo(nMon,1) = colMonth(imon(nt)-1)
> end if
>
> if (pltMovie) then
> draw(plot)
> frame(wks)
> end if
>
> xBegin= pc1(nt)
> yBegin= pc2(nt)
> end do
>
> plMark@gsMarkerColor = "black" ; indicate last point
> plMark@gsMarkerSizeF = 2.5*plMark@gsMarkerSizeF ; make larger
> plot@$unique_string("dum")$ = gsn_add_polymarker(wks, plot, xBegin,
> yBegin, plMark)
>
> if (.not.pltMovie) then
> if (label_color) then ; fancy coloring of months
> getvalues plot
> "trXMinF" : xmin
> end getvalues
> xinc = (xmin*-2) / (nMon+1+2.) ; total number months = nMon+1. +2
> ; for spacing on ends
> txres@txJust = "CenterLeft"
> txres@txFontHeightF = 0.013 ; size of font for printed month
> xstart = xmin+(1.25*xinc)
> if (nMon.gt.0) then
> do n=0, nMon
> txres@txFontColor = monInfo(n,1)
> plot@$unique_string("dum")$ = gsn_add_text(wks, plot, monInfo(n,0)
> ,xstart , -3.75, txres)
> xstart = xstart+xinc
> end do
> end if
> end if
> draw(plot)
> frame(wks)
> end if
> end
> ----
>
> Question 2: regarding the mjoclivar_16.ncl. I run the code to plot MJO
> Phase during 2011. it run successfully. However, I think the result is
> little bit abnormal because the negative anomaly of OLR propagates
> westward. It is commonly known that the MJO propagates estward. I send
> the plot example in the email. The code for plotting is:
>
> ;***********************************************************
> ; Generate life cycle composites based upon daily phase space
> ; If the MJO_INDEX is < 1.0 it is not included
> ;***********************************************************
> 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 = -15
> latN = 15
>
> ymdStrt = 20110101 ; start yyyymmdd
> ymdLast = 20111231 ; last
>
> yrStrt = ymdStrt/10000
> yrLast = ymdLast/10000
>
> pltSubTitle = "Anomalous: OLR, U850, V850"
>
> pltDir = "./" ; plot directory
> pltType = "ps"
> pltName = "mjoclivar" ; yrStrt+"_"+yrLast
> diri = "e:/RData/MJO_cal/Anom/" ; input directory
>
> filo = "olr.day.anomalies.2001-2011.nc"
> filu = "uwnd.day.850.anomalies.2001-2011.nc"
> filv = "vwnd.day.850.anomalies.2001-2011.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 indicies (subscripts) corresponding to the start/end times
> ;***********************************************************
>
> f = addfile (diri+filu , "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)
> delete(YMD )
>
> time = f->time(iStrt:iLast) ; days since ...
> u = f->UWND_anom(iStrt:iLast,{latS:latN},:)
> ;***********************************************************
> ; Read anomalies frpm other fields
> ;***********************************************************
> f = addfile (diri+filv , "r")
> v = f->VWND_anom(iStrt:iLast,{latS:latN},:)
>
> f = addfile (diri+filo , "r")
> x = f->OLR_anom(iStrt:iLast,{latS:latN},:)
>
> dimx = dimsizes( x )
> ntim = dimx(0)
> nlat = dimx(1)
> mlon = dimx(2)
> ;************************************************
> ; Apply the band pass filter to the original anomalies
> ;************************************************
> x = wgt_runave_leftdim (x, wgt, 0)
> u = wgt_runave_leftdim (u, wgt, 0)
> v = wgt_runave_leftdim (v, wgt, 0)
>
> ;***********************************************************
> ; Open PC components file created in 'mjo_14.ncl'
> ;***********************************************************
> diri = "e:/RData/MJO_cal/Anom/" ; input directory
> fMJO = "MJO_PC_INDEX.nc" ; created in mjo_14.ncl
> f = addfile (diri+fMJO, "r")
>
> ;***********************************************************
> ; Find the indices corresponding to the start/end times
> ;***********************************************************
> TIME = f->time ; days since ...
> YMD = cd_calendar(TIME, -2) ; entire (time,6)
> printVarSummary(TIME)
>
> iStrt = ind(YMD.eq.ymdStrt) ; index start
> iLast = ind(YMD.eq.ymdLast) ; index last
> delete(TIME)
> delete(YMD )
>
> ;***********************************************************
> ; Read the data for the desired period
> ;***********************************************************
> pc1 = f->PC1(iStrt:iLast)
> pc2 = f->PC2(iStrt:iLast)
> mjo_indx= f->MJO_INDEX(iStrt:iLast)
>
> time = f->time(iStrt:iLast)
> ymdhms = cd_calendar(time, 0)
>
> imon = floattoint( ymdhms(:,1) ) ; convenience
> iday = floattoint( ymdhms(:,2) ) ; subscripts must be integer
>
> ;***********************************************************
> ; Place each array into an appropriate array
> ;***********************************************************
> 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(pc2,pc1)*r2d ; phase space
> nn = ind(ang.lt.0)
> ang(nn) = ang(nn) + 360 ; make 0 to 360
>
> nDays = new (nPhase, "integer")
> pLabel = "P"+ispan(1,nPhase,1)+": "
>
> ;------------------------------------------------------------
> ; PLOTS
> ;------------------------------------------------------------
> pltPath = pltDir+pltName
>
> wks = gsn_open_wks(pltType,pltPath)
> gsn_define_colormap(wks,"ViBlGrWhYeOrRe")
> plot = new(nPhase,graphic) ; create graphic array
>
> 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@mpFillOn = False ; turn off map fill
> res@mpMinLatF = latS ; zoom in on map
> res@mpMaxLatF = latN
> res@mpCenterLonF = 210.
> res@cnFillOn = True ; turn on color fill
> res@cnLinesOn = False ; True is default
> res@cnLineLabelsOn = False ; True is default
> res@lbLabelBarOn = False ; turn off individual lb's
> res@gsnScalarContour = True ; contour 3rd array
> res@gsnMajorLatSpacing = 15
> res@gsnMajorLonSpacing = 60
> res@tmXBLabelFontHeightF = 0.01
> res@tmYLLabelFontHeightF = 0.01
>
> ; common contours
> ;mnmxint = nice_mnmxintvl( min(x) , max(x), 16, False)
> res@cnLevelSelectionMode = "ManualLevels"
> res@cnMinLevelValF = -40 ; -100; mnmxint(0)
> res@cnMaxLevelValF = 40 ; 80; mnmxint(1)
> res@cnLevelSpacingF = 5 ; 20; mnmxint(2)
> ;print(res)
>
> res@vcMinDistanceF = 0.01 ; thin the vector density
> res@vcRefMagnitudeF = 2.0 ; define vector ref mag
> res@vcRefLengthF = 0.025 ; define length of vec ref
> res@vcRefAnnoOrthogonalPosF = -1.0 ; move ref vector
> res@vcRefAnnoArrowLineColor = "black" ; change ref vector color
> res@vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref
>
> ; 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@lbLabelFontHeightF = 0.01
> resP@gsnPanelBottom = 0.05 ; add some space at bottom
> resP@pmLabelBarWidthF = 0.8 ; label bar width
> resP@pmLabelBarHeightF = 0.05
> resP@gsnPanelFigureStringsFontHeightF = 0.0125 ; bit larger than default
> ;resP@pmLabelBarOrthogonalPosF = 0.015 ; move labelbar up a bit
>
> txres = True
> txres@txFontHeightF = 0.01
> txid = gsn_create_text(wks, pltSubTitle, txres)
>
> amres = True
> ;amres@amParallelPosF = 0.575
> amres@amOrthogonalPosF = 0.75
> amres@amJust = "CenterCenter"
> ;amres@amResizeNotify = True
>
> ;*******************************************
> ; Loop over each phase
> ;*******************************************
> res@gsnLeftString = ""
> res@gsnRightString = ""
> do nSeason=1,2
> if (nSeason.eq.1) then
> resP@txString = yrStrt+"-"+yrLast+": Jan to Dec"
> else
> resP@txString = yrStrt+"-"+yrLast+": Jan to Dec"
> end if
> do n=0,nPhase-1
> if (nSeason.eq.1) then
> nt = ind(mjo_indx.gt.1.0 .and. \
> (imon.ge.1 .and. imon.le.12).and. \
> ang.ge.angBnd(0,n) .and. ang.lt.angBnd(1,n))
> else
> nt = ind(mjo_indx.gt.1.0 .and. \
> (imon.ge.1 .or. imon.le. 12).and. \
> ang.ge.angBnd(0,n) .and. ang.lt.angBnd(1,n))
> end if
> if (.not.all(ismissing(nt))) then
> xAvg = dim_avg_Wrap( x(lat|:,lon|:,time|nt) )
> uAvg = dim_avg_Wrap( u(lat|:,lon|:,time|nt) )
> vAvg = dim_avg_Wrap( v(lat|:,lon|:,time|nt) )
> nDays(n) = dimsizes(nt)
>
> res@tmXBLabelsOn = False ; do not draw lon labels
> res@tmXBOn = False ; lon tickmarks
> if (n.eq.(nPhase-1)) then ;
> res@tmXBLabelsOn = True ; draw lon labels
> res@tmXBOn = True ; tickmarks
> end if
>
> plot(n) = gsn_csm_vector_scalar_map_ce(wks,uAvg,vAvg,xAvg,res)
> end if
> delete(nt) ; will change next iteration
> end do
>
> ;;ann3 = gsn_add_annotation(plot(0), txid, amres)
> resP@gsnPanelFigureStrings= pLabel+nDays
> gsn_panel(wks,plot,(/nPhase,1/),resP) ; now draw as one plot
> end do
>
> system("gv "+pltPath+"."+pltType)
> end
>
> --------------------
>
> Best regards,
> JuKY
>
>
>
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu May 31 22:25:11 2012