Re: Panel plot with loop

Liz,

The error is indicating something is wrong with the "nfp" index of this line:

WrfPrdiff(nfp,:,:)=xdiff ; assign xdiff to WrfPr variable

nfp is an index that goes from 0 to nfilp-1. What is "nfilp-1" equal to?

The leftmost dimension of wrfPrdiff is 6, so my guess is that the do loop is reaching "nfp=6" which will give you an error, since the leftmost dimension can only have indexes from 0 to 5.

--Mary

On Feb 11, 2014, at 11:17 AM, burakows@ucar.edu 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"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
>
> ;======================================================================
> ; The main code
> ;======================================================================
>
> begin
>
> ;---Specify interpolation to be used
> method = "bilinear"
>
> wgtFileDir = "/glade/p/work/burakows/plot/NCL/ESMF_regridding/"
> wgtFileName = "PRISM_4kmM2_to_WRF."+method+"_wgt.nc"
> wgtFilePath = wgtFileDir+wgtFileName
> print("wgtFilePath="+wgtFilePath)
>
> ;---PRISM files [input]
> varp = "tmax" ; PRISM name
> dirp = "/glade/p/work/burakows/prism/data/"
> filp = systemfunc("cd "+dirp+" ; ls PRISM_"+varp+"*bil.nc")
> nfilp = dimsizes(filp)
> print(filp(0:nfilp-1)) ; 1st 5 file
>
> ;---WRF directory
> varw = "T2" ; WRF name
> dirw =
> "/glade/p/work/burakows/plot/albedo_CLMV351current_200811-200904/"
>
> varm = "T2max" ; WRF monthly name
> prtFlg = True
>
> ;---Create new ncl variable to hold panel plot contents
> plot_diff = new(6,graphic) ; from p.80 of NCL manual, right
> before looping (panel1e.ncl)
> WrfPrdiff = new((/6,192,132/),float,1e20) ; create new variable
> WrfPr(time,y,x)
>
> print(dimsizes(WrfPrdiff))
>
> ;----------------------------------------------------------------------
> ; Loop through six PRISM and WRF monthly files to calculate difference
> ;----------------------------------------------------------------------
>
> do nfp=0,nfilp-1
> yyyymm = toint( str_get_field(filp(nfp),5,"_.") ) ; parse
> PRISM name for time
>
> yyyy = yyyymm/100 ; current
> year as integer
> mm = yyyymm-(yyyy*100) ; current month
>
> fp = addfile(dirp+filp(nfp), "r")
> XP = fp->z ; original
> PRISM variable (all 'z')
> xp = ESMF_regrid_with_weights(XP,wgtFilePath,False); create new
> variable on WRF grid
>
> printVarSummary(XP)
> printVarSummary(xp)
>
> ; matching WRF file
> filw = "avg_"+varm+"_"+yyyy+"-"+sprinti("%0.2i",mm)+".nc"
>
> if (isfilepresent(dirw+filw)) then
> fw = addfile(dirw+filw, "r")
> xw = fw->$varm$ ; WRF
> monthly variable
> ;xw@lat2d = fw->XLAT(0,:,:)
> ;xw@lon2d = fw->XLONG(0,:,:)
>
> printVarSummary(xw)
> ;----------------------------------------------------------------------
> ; Subtract PRISM from WRF
> ;----------------------------------------------------------------------
>
> fwo = addfile(dirw+"wrfout_d03_2008-11-01_00:00:00.nc","r"); Open any
> wrfout daily file to obtain
> ; lat and lon
> xwo = fwo->T2 ; read T2 and preserve metadata
>
> xdiff = xw(0,:,:)-xp ; Avg monthly WRF minus Prism for month i
> print(dimsizes(xdiff))
> xdiff@long_name = "difference: "+yyyymm ; Assign long name
> to difference
> xdiff@units = "Degrees C" ; Assign units to difference
> printVarSummary(xdiff)
>
>
> xdiff@lat2d = fwo->XLAT(0,:,:) ; Assign lat to xdiff
> xdiff@lon2d = fwo->XLONG(0,:,:) ; Assign lon to xdiff
>
> ;print("xdiff: min="+min(xdiff(nfp,:,:))+" max="+max(xdiff(nfp,:,:)))
>
> minlat = min((/min(xdiff@lat2d),min(fp&y)/))
> maxlat = max((/max(xdiff@lat2d),max(fp&y)/))
>
> minlon = min((/min(xdiff@lon2d),min(fp&x)/))
> maxlon = min((/max(xdiff@lon2d),max(fp&x)/))
>
> bndadd = 0.10 ; spacing around the plot edges
>
> WrfPrdiff(nfp,:,:)=xdiff ; assign xdiff to WrfPr variable
>
> plot_diff(nfp) = gsn_csm_contour_map(wks,WrfPrdiff(nfp),res)
> ; create plot of month nfp in plot_diff
>
> if (prtFlg) then ; first time only
> printVarSummary(XP)
> printVarSummary(xp)
> printVarSummary(xw)
> prtFlg = False
> end if ; nfp
> end if ; isfilepresent
>
> end do
>
> printVarSummary(WrfPrdiff)
> ;----------------------------------------------------------------------
> ; Plotting options section
> ;----------------------------------------------------------------------
>
> pltType = "ps" ; plot type
> pltDir = "./" ; plot directory
> pltName = "Diff_WRF_PRISM_2009" ; plot name (ps file)
> pltPath = pltDir+pltName ; plot path
>
> wks = gsn_open_wks(pltType,pltPath) ; create
> workstation for ps file
>
> gsn_define_colormap(wks,"temp_diff_18lev") ; define color table
>
> res = True
>
> res@gsnDraw = False
> res@gsnFrame = False
>
> res@cnFillOn = True ; color plot desired
> res@cnLinesOn = False ; turn off contour lines
> res@cnLineLabelsOn = False ; turn off contour labels
> res@cnFillMode = "RasterFill" ; turn raster on
>
> res@lbLabelBarOn = False ; Will turn on in panel
> later
> res@lbOrientation = "Horizontal" ; Horizontal label bar
>
> res@mpFillOn = False
> res@mpOutlineOn = True
> res@mpOutlineBoundarySets = "AllBoundaries"
> res@mpProjection = "CylindricalEquidistant"
>
> res@mpLimitMode = "LatLon" ; required
> res@mpMinLatF = minlat-bndadd
> res@mpMaxLatF = maxlat+bndadd
> res@mpMinLonF = minlon-bndadd
> res@mpMaxLonF = maxlon+bndadd
> res@mpCenterLonF = (minlon + maxlon)*0.5
> res@mpCenterLatF = (minlat + maxlat)*0.5
>
> res@gsnMajorLatSpacing = 2
> ;res@gsnMajorLonSpacing = 2
> res@gsnMinorLonSpacing = 2
>
> res@cnLinesOn = False
> res@cnLineLabelsOn = False
> res@cnFillOn = True
> res@cnFillMode = "RasterFill"
>
> res@gsnAddCyclic = False ; regional grid (changes central meridian)
>
> ;----------------------------------------------------------------------
> ; Panel Plotting options section
> ;----------------------------------------------------------------------
>
> pres = True
> pres@gsnMaximize = True
> pres@gsnPanelLabelBar = True
> pres@gsnCenterString = "WRF minus PRISM" + yyyymm
> pres@txString = "Monthly averaged WRF T2 minus PRISM"
>
> gsn_panel(wks,plot_diff,(/2,3/),pres)
> end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Feb 12 09:56:59 2014