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"
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 -rotate '-90<' "+case+"_ivt.png"
cmd = "convert "+case+"_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
Received on Thu Sep 19 15:00:16 2013
This archive was generated by hypermail 2.1.8 : Tue Oct 01 2013 - 14:41:43 MDT