Blank Image with wrf_map_overlays() ???

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