Re: Blank Image with wrf_map_overlays() ???

From: Adam Phillips <asphilli_at_nyahnyahspammersnyahnyah>
Date: Fri Sep 20 2013 - 16:01:05 MDT

Hi A.J.,
Typically we recommend that all questions with the WRF functions be sent
to wrfhelp@ucar.edu. But before you write to them it might be best to
separate and plot out each of your overlaid plots to see which ones work
and which ones don't. Hopefully doing that will provide clues as to
where the problem is. You might also want to use printMinMax, print, and
printVarSummary on the arrays going into the plotting functions to see
if anything looks amiss.

If the above does not help, I'd recommend you contact wrfhelp@ucar.edu..
Adam

On 09/19/2013 03:00 PM, A.J. Eiserloh wrote:
> Hi all,
>
> I am trying to do a simple wrf_map_overlays, but I keep getting a
> blank image and I do not understand why. The field that I am trying to
> contour and use for vectors I get from a seperate .nc file. Also I
> should mention that I am not using winds for the vectors. They are the
> x- and y- components of my vector value (IVT : Integrated vapor
> transport). All the data from the .nc files is there because I checked
> it with ncdump. Also, I get no errors out of the output from NCL.
>
> here is my script below:
>
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
>
> begin
>
> case ="3dvar_t1"
> if (case.eq."run3.1" .or. case.eq."n3") then
> cmd ="ls -1 /data2/eiserloh/data/AR/" + case + "/wrfout_d02_2012*"
> ;frost
> else
> cmd ="ls -1 /data2/eiserloh/data/" + case + "/wrfout_d02_2012*"
> ;thunder
> end if
>
> files = systemfunc(cmd) + ".nc"
>
> a = addfiles(files,"r")
> time = wrf_user_list_times(a) ; get times in the files
>
> lat2d=wrf_user_getvar(a[0],"lat",0)
> lon2d=wrf_user_getvar(a[0],"lon",0)
>
> ivt_file=case+"_ivt_6hrly.nc <http://ivt_6hrly.nc>"
> b = addfile(ivt_file,"r")
>
> times = b->times
> IVT = b->IVT
> IVT_x = b->IVT_x
> IVT_y = b->IVT_y
>
> ; IVT@lat2d = lat2d
> ; IVT@lon2d = lon2d
> ; IVT_x@lat2d = lat2d
> ; IVT_x@lon2d = lon2d
> ; IVT_y@lat2d = lat2d
> ; IVT_y@lon2d = lon2d
>
> ntimes = dimsizes(times)
>
> ; We generate plots, but what kind do we prefer?
> type = "ps"
> ; type = "png"
> type@wkOrientation = "landscape"
>
> wks = gsn_open_wks(type,case+"_ivt")
> gsn_define_colormap(wks,"BlAqGrYeOrRe")
>
> ; Set some basic resources
> res = True
> res@gsnDraw = False
> res@gsnFrame = False
>
> res@MainTitle = " "
> res@Footer = False
> MTOPosF = -0.05
>
> pltres = True
> mpres = True
> colr = "black"
> mpres@gsnMaximize = True
> mpres@mpGeophysicalLineColor = colr
> mpres@mpNationalLineColor = colr
> mpres@mpUSStateLineColor = colr
> mpres@mpGridLineColor = colr
> mpres@mpLimbLineColor = colr
> mpres@mpPerimLineColor = colr
> mpres@mpOutlineBoundarySets = "GeophysicalandUSStates"
> mpres@mpGeophysicalLineThicknessF = 1.5
> mpres@mpUSStateLineThicknessF = 1.5
> mpres@mpGridLineThicknessF = 1.0
> mpres@mpLimbLineThicknessF = 1.0
> mpres@mpNationalLineThicknessF = 1.5
> mpres@mpUSStateLineThicknessF = 1.5
> mpres@mpGridLineDashPattern = 1.5
> mpres@mpDataBaseVersion = "MediumRes"
> mpres@mpDataSetName = "Earth..4"
> pltres@PanelPlot = True
> pltres@LatLonOverlay = True
>
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> ; What times and how many time steps are in the data set?
>
> hour_int=6
>
> ; do it = 0,ntimes-1 ; TIME LOOP
> do it =0,1
>
> f=it*hour_int
> print("Working on time: " + time(f) )
>
> res@TimeLabel = times(it) ; Set Valid time to use on plots
>
>
> ; u = wrf_user_getvar(a[f],"ua",0) ; u averaged to mass
> points
> ; v = wrf_user_getvar(a[f],"va",0) ; v averaged to mass
> points
> ; ws = ((u^2)+(v^2))^(0.5) ; Wind speed
> ter = wrf_user_getvar(a[f],"ter",0)
> ; p = wrf_user_getvar(a[f],"pressure",0)
>
> ;mandatory pressure levels
> ; mplevs=(/1000.,925.,850.,700.,500.,300./)
>
> ;get u at mandatory plevs for vectors for plot
> ; u_mpls=wrf_user_intrp3d(u,p,"h",mand_pres_levs,0.,False)
> ; v_mpls=wrf_user_intrp3d(v,p,"h",mand_pres_levs,0.,False)
>
> ; u_bar=dim_sum_n(u_mpls,0)
> ; v_bar=dim_sum_n(v_mpls,0)
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ; Plotting Time
> slp = wrf_user_getvar(a[f],"slp",0) ; slp
> res@InitTime = False
> res@ValidTime = False
> res@NoHeaderFooter = True
> tt=slp
> tt@description = " Valid " + times(it)
> tt@units = "UTC"
> opts=res
> opts@ContourParameters = (/1., 2., 1./)
> opts@cnFillColors = -1
> opts@lbTitleOn = False
> opts@cnInfoLabelOn = False
> opts@sfXArray = lon2d
> opts@sfYArray = lat2d
> title = wrf_contour(a[f],wks,tt,opts)
> delete(opts)
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> ; Plotting options IVT contour fill
> opts = res
> opts@cnFillOn = True
> opts@ContourParameters = (/ 0., 800., 50. /)
> opts@lbOrientation = "Vertical"
> opts@pmLabelBarSide = "Right"
> opts@pmLabelBarOrthogonalPosF= -0.04
> opts@lbTitleOn = False
> opts@FieldTitle = IVT@description
> opts@sfXArray = lon2d
> opts@sfYArray = lat2d
> contour = wrf_contour(a[f],wks,IVT(it,:,:),opts)
> delete(opts)
>
> ; Plotting Opts for Terrain Contours
> opts=res
> opts@cnFillOn=False
> opts@ContourParameters=(/0.,5000.,250./)
> opts@FieldTitle = "Elevation (m)"
> opts@cnLineLabelInterval = 0
> opts@sfXArray = lon2d
> opts@sfYArray = lat2d
> terrain = wrf_contour(a[f],wks,ter,opts)
> delete(opts)
>
> ; Overplot Vectors
> opts = res
> ; opts@vcGlyphStyle = "LineArrow"
> opts@FieldTitle = "IVT Vectors" ; overwrite Field Title
> opts@NumVectors = 25 ; wind barb density
> opts@vfXArray = lon2d
> opts@vfYArray = lat2d
> vector = wrf_vector(a[f],wks,IVT_x(it,:,:),IVT_y(it,:,:),opts)
> delete(opts)
>
>
> ; MAKE PLOTS
> plot =
> wrf_map_overlays(a[f],wks,(/contour,vector,terrain,title/),pltres,mpres)
> draw(plot)
> frame(wks)
>
> end do ; END OF TIME LOOP
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ; cmd = "convert "+case+"_ivt.ps <http://ivt.ps> -rotate '-90<'
> "+case+"_ivt.png"
> cmd = "convert "+case+"_ivt.ps <http://ivt.ps> "+case+"_ivt.png"
> system(cmd)
>
> end
>
> Thanks,
>
> --
> Arthur J. Eiserloh, Jr.
> San Jose State University
> Graduate Student
> Dept. of Meteorology and Climate Science
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
______________________________________________________________
Adam Phillips                                asphilli@ucar.edu
NCAR/Climate and Global Dynamics Division       (303) 497-1726
P.O. Box 3000				
Boulder, CO 80307-3000    http://www.cgd.ucar.edu/cas/asphilli

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Sep 20 16:01:17 2013

This archive was generated by hypermail 2.1.8 : Tue Oct 01 2013 - 14:41:43 MDT