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
----
thank you
JuKY
________________________________
From: Dennis Shea <shea@ucar.edu>
To: juki juki <juky_emc2@yahoo.com>
Cc: "ncl-talk@ucar.edu" <ncl-talk@ucar.edu>
Sent: Thursday, May 31, 2012 1:06 PM
Subject: Re: MJO clivar
Whenever you have an error message like the one you had,
please use "printVarSummary" on the variable(s) involved.
Look at the dimension names and, especially, the sizes.
fatal:Dimension sizes on right hand side of assignment do not match
dimension sizes of left hand side
========================================================
Change
lat = f[:]->lat
lon = f[:]->lon
to
lat = f[0]->lat
lon = f[0]->lon
On 5/30/12 7:31 PM, juki juki wrote:
> 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 specifications
> ;-------------------------------------------------------------
> NC = True ; create netCDF?
> PLOT = False ; sample plots?
>
> ymdStrt = 20010101 ; start yyyymmdd
> ymdLast = 20111231 ; last
> yrStrt = ymdStrt/10000
> yrLast = ymdLast/10000
> lev = 850 ; desired level
> nhar = 4 ; number of fourier comp
>
> var = "uwnd" ; name of file
> vName = "UWND" ; name for plots
> fili = var+".day.mean.nc" ; input file
>
> ;========================
> ; get list of all files and open as "one big file"
> ;========================
> all_files = systemfunc ("ls e:/Rdata/MJO_cal/Data/uwnd*.nc")
> f = addfiles (all_files, "r") ; note the "s" of addfile
> print(all_files)
> ;========================
> ; choose how files are combined and read in variable across files
> ;========================
> ListSetType (f, "cat") ; concatenate or "merge" (default)
>
> if (NC) then
> diro= "e:/RData/MJO_cal/Anom/" ; output dir
> filo= var+".day."+lev+".anomalies."+yrStrt+"-"+yrLast+".nc" ;
> output file
> end if
>
> if (PLOT) then
> wksType = "ps"
> wksName = "mjoclivar" ; "mjo."+yrStrt+"_"+yrLast
> end if
>
> ;***********************************************************
> ; Read user specified time and create required yyyyddd
> ;***********************************************************
> ; f = addfile (diri+fili , "r")
>
> time = f[:]->time ; all times on file
> ymd = cd_calendar(time, -2) ; yyyymmdd
> iStrt = ind(ymd.eq.ymdStrt) ; index start
> iLast = ind(ymd.eq.ymdLast) ; index last
> delete(time)
> delete(ymd)
>
> ;***********************************************************
> ; Read user specified time and create required yyyyddd
> ;***********************************************************
> time = f[:]->time(iStrt:iLast) ; time:units = "hours
> since"
> TIME = cd_calendar(time, 0) ; type float
> year = floattointeger( TIME(:,0) )
> month = floattointeger( TIME(:,1) )
> day = floattointeger( TIME(:,2) )
> ddd = day_of_year(year, month, day)
> yyyyddd = year*1000 + ddd ; needed for input
>
> ;***********************************************************
> ; Read data: short2flt
> ;***********************************************************
> x = short2flt(f[:]->$var$(iStrt:iLast,{lev},:,:)) ;
> convert to float
> printVarSummary(x)
>
> ;***********************************************************
> ; Compute daily climatology: raw and then 'smoothed'
> ;***********************************************************
> xClmDay = clmDayTLL(x, yyyyddd) ; daily climatology at each grid
> point
> printVarSummary(xClmDay)
>
> ;***********************************************************
> ; Compute smoothed daily climatology using 'nhar' harmonics
> ;***********************************************************
> xClmDay_sm = smthClmDayTLL(xClmDay, nhar)
> printVarSummary(xClmDay_sm)
>
> ;***********************************************************
> ; Compute daily anomalies using raw and smoothed climatologies
> ;***********************************************************
> xAnom = calcDayAnomTLL (x, yyyyddd, xClmDay)
> printVarSummary(xAnom)
> printMinMax(xAnom, True)
>
> xAnom_sm = calcDayAnomTLL (x, yyyyddd, xClmDay_sm)
> xAnom_sm@long_name = "Anomalies from Smooth Daily Climatology"
> printVarSummary(xAnom_sm)
> printMinMax(xAnom_sm, True)
>
> delete( x ) ; no longer needed
>
>
> ;***********************************************************
> ; Create netCDF: convenience use 'simple' method
> ;***********************************************************
>
> dimx = dimsizes(xAnom)
> ntim = dimx(0)
> nlat = dimx(1)
> mlon = dimx(2)
>
> if (NC) then
>
> lat = f[:]->lat
> lon = f[:]->lon
>
> system("/bin/rm -f "+diro+filo) ; rm any pre-exist file,
> if any
> fnc = addfile (diro+filo, "c")
>
> filAtt = 0
> filAtt@title = vName+": Daily Anomalies:
> "+yrStrt+"-"+yrLast
> filAtt@source_file = fili
> filAtt@creation_date = systemfunc("date")
> fileattdef( fnc, filAtt ) ; copy file attributes
>
> setfileoption(fnc,"DefineMode",True)
> varNC = vName+"_anom"
> varNC_sm = vName+"_anom_sm"
>
> dimNames = (/"time", "lat", "lon"/)
> dimSizes = (/ -1 , nlat, mlon/)
> dimUnlim = (/ True , False, False/)
> filedimdef(fnc,dimNames,dimSizes,dimUnlim)
>
> filevardef(fnc, "time" ,typeof(time),getvardims(time))
> filevardef(fnc, "lat" ,typeof(lat) ,getvardims(lat))
> filevardef(fnc, "lon" ,typeof(lon) ,getvardims(lon))
> filevardef(fnc, varNC ,typeof(xAnom) ,getvardims(xAnom))
> filevardef(fnc, varNC_sm,typeof(xAnom) ,getvardims(xAnom))
>
> filevarattdef(fnc,"time" ,time) ; copy time
> attributes
> filevarattdef(fnc,"lat" ,lat) ; copy lat
> attributes
> filevarattdef(fnc,"lon" ,lon) ; copy lon
> attributes
> filevarattdef(fnc,varNC ,xAnom)
> filevarattdef(fnc,varNC_sm,xAnom_sm)
>
> fnc->time = (/time/)
> fnc->lat = (/lat/)
> fnc->lon = (/lon/)
> fnc->$varNC$ = (/xAnom /)
> fnc->$varNC_sm$ = (/xAnom_sm/)
> end if
>
> ;************************************************
> ; plotting parameters
> ;************************************************
> if (PLOT) then
> LAT = (/ 60, 45, 5, -5, -45, 60 /)
> LON = (/270, 30, 90, 90, 180, 0 /)
> nPts = dimsizes( LAT )
>
> plot = new ( nPts, graphic)
> data = new ( (/2,366/), typeof(xClmDay), getFillValue(xClmDay))
>
> wks = gsn_open_wks (wksType,wksName)
>
> res = True ; plot mods
> desired
> res@gsnDraw = False
> res@gsnFrame = False
> res@trXMinF = 1
> res@trXMaxF = 366
> ;res@tiMainString = ""
> res@xyLineThicknesses = (/1.0, 2.0/) ; make 2nd
> lines thicker
> res@xyLineColors = (/"blue","red"/) ; change line
> color
> res@xyMonoDashPattern = True ; all solid
>
> do np=0,nPts-1
> data(0,:) = xClmDay(:,{LAT(np)},{LON(np)})
> data(1,:) = xClmDay_sm(:,{LAT(np)},{LON(np)})
> res@gsnCenterString = "lat="+LAT(np)+" lon="+LON(np)
> plot(np) = gsn_csm_y (wks,data,res) ; create plot
> end do
> resP = True ; modify the
> panel plot
> resP@txString = vName+": Raw/Smooth Daily
> Climatology: "+yrStrt+"-"+yrLast
> resP@gsnMaximize = True
> resP@gsnPaperOrientation = "portrait"
> gsn_panel(wks,plot,(/(nPts/2),2/),resP) ; now
> draw as one plot
>
> ;==========
> ; Plot anomalies for an arbitrarily selected near equatorial
> location
> ; Time: Oct 1, 1996 to April 1,1997 [arbitrary selection]
> ;==========
> LATX = 0
> LONX = 90
>
> yyyymmdd = cd_calendar(time, -2)
> ;;yrfrac = yyyymmdd_to_yyyyfrac (yyyymmdd, 0)
> ;;delete(yrfrac@long_name)
>
> xAnom@long_name = "Anomalies from Raw" ; short labels for plot
> xAnom_sm@long_name = "Anomalies from Smooth"
>
> ntBegin = ind(yyyymmdd.eq.19961001)
> ntEnd = ind(yyyymmdd.eq.19970401)
> monthLabels = (/1,4,7,10/)
> monNam = (/"Jan","Feb","Mar","Apr","May","Jun" \
> ,"Jul","Aug","Sep","Oct","Nov","Dec" /)
> xVal = new(ntim, typeof(xAnom&time) , "No_FillValue") ;
> bigger than
> xLab = new(ntim, "string", "No_FillValue") ; needed
> xValm = new(ntim, typeof(xAnom&time) , "No_FillValue") ;
> bigger than
>
> ntm = -1
> cr = inttochar(10) ; carriage
> return
> do nt=ntBegin,ntEnd
> if (day(nt).eq.1) then
> ntm = ntm + 1
> xVal(ntm) = xAnom&time(nt)
> xLab(ntm) = monNam(month(nt)-1)
> if (month(nt).eq.1) then
> xLab(ntm) = xLab(ntm) + cr +sprinti("%0.4i", year(nt))
> end if
> end if
> end do
>
> rxy = True
> rxy@gsnDraw = False
> rxy@gsnFrame = False
> rxy@gsnYRefLine = 0.0 ; create a reference line
> rxy@gsnAboveYRefLineColor = "red" ; above ref line fill red
> rxy@gsnBelowYRefLineColor = "blue" ; below ref line fill blue
> rxy@xyLineThicknessF = 2.0
> rxy@vpHeightF = 0.4 ; resize
> rxy@vpWidthF = 0.8
> rxy@tmXBMode = "Explicit"
> rxy@tmXBValues = xVal(0:ntm)
> rxy@tmXBLabels = xLab(0:ntm)
>
> plot(0) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \
> ,xAnom(ntBegin:ntEnd,{0},{90}),rxy)
> plot(1) = gsn_csm_xy (wks,time(ntBegin:ntEnd) \
> ,xAnom_sm(ntBegin:ntEnd,{0},{90}),rxy)
> resP@txString = vName+": Anomalies: (0,90E)"
> gsn_panel(wks,plot(0:1),(/2,1/),resP)
>
> end if
> end
>
>
>
> ------------------------------------------------------------------------
> *From:* Dennis Shea <shea@ucar.edu>
> *To:* juki juki <juky_emc2@yahoo.com>
> *Cc:* "ncl-talk@ucar.edu" <ncl-talk@ucar.edu>
> *Sent:* Thursday, May 17, 2012 10:16 AM
> *Subject:* Re: [ncl-talk] MJO clivar
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu May 31 17:30:32 2012
This archive was generated by hypermail 2.1.8 : Wed Jun 06 2012 - 15:17:44 MDT