Y-Axis Label Problem On WRF Cross Section

From: Mike Pirhalla <glacierbayncl_at_nyahnyahspammersnyahnyah>
Date: Mon Dec 02 2013 - 23:48:14 MST

Hi NCL Community,

I have a working cross section of PM data, wind, and temperature. However,
I'd like to cut off the graph at 3km in height and below. When I do this by
putting in a value at opts@trYMaxF, the Yaxis labels get messed up, and I
think it's inaccurately labeling the y-axis because the terrain height and
other contours don't match as the once did. I'd also like the make the
tickmarks more prevalent on the y-axis. Attached is my code and the plots,
one with a cut off yaxis at opts@trYMaxF=24 (? Was experimenting with
values) and other with the entire WRF column. Any help is appreciated,
thanks!

 

Michael Pirhalla

Graduate Research Assistant

University of Alaska Fairbanks

 

 

;--------------------------------------------------------------------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.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/shea_util.ncl"

load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

;load "color_map.ncl"

;---------------------------------------------------------------------------

 

begin

;What type of plot do we want?

type = "x11"

;type = "ps"

;type = "pdf"

 

wks = gsn_open_wks(type,"WRF_PM_Cross_Section")

gsn_define_colormap(wks,"BlAqGrYeOrReVi200")

 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; READ WRF OUTPUT FILE AND INCLUDE ALL CASES

 

; all_files = systemfunc ("ls
/center/d/molders/WITH/wrfout_d01_200*")+".nc"

; print(all_files)

; a = addfiles(all_files,"r")

; print(a)

; ListSetType (a, "cat") ; Combine all files

; idim=dimsizes(all_files)

 

a =
addfile("/center/d/molders/WITH/wrfout_d01_2008-07-19_00:00:00"+".nc","r")

 

;---------------------------------------------------------------------------

  FirstTime = True

  times = wrf_user_getvar(a,"times",-1) ; get times in the file

  ntimes = dimsizes(times) ; number of times in the file

  print("ntimes"+ntimes)

;---------------------------------------------------------------------------
-

; Set some basic resources

  res = True

  res@MainTitle = "REAL-TIME WRF" ; main title of plot

  res@Footer = True ; footer information

 

  pltres = True

;---------------------------------------------------------------------------

--
; HOURLY LOOP OF PARTICULATE MATTER
 
do it = 0,ntimes-1  ; loop over all hours
;do it=11,11 ; choose specific hour
 
  print("Working on time: " + times(it) )
  res@TimeLabel = times(it)   ; Set Valid time to use on plots
 
; First get the variables we will need        
    tc  = wrf_user_getvar(a,"tc",it)     ; T in C
    z   = wrf_user_getvar(a,"z",it)     ; grid point height
    hgt = wrf_user_getvar(a,"HGT",it)
 
   if(FirstTime) then       ; get height for labels
                zmin = 0. 
                zmax = max(z) / 1000  ; make to km
                nz = floattoint(zmax/2 + 1)
                FirstTime = False
   end if
;z_values=fspan(zmin,zmax,100)
;print(z_values)
;exit
 
; Get PM10 data from WRF/Chem
    pm = a->PM10(it,:,:,:)
    pm@description = "PM~B~10~N~ Concentration"
    pm@units = "~F33~m~F21~g/m~S~3~N~"
 
; Get the wind data
uu=wrf_user_getvar(a,"ua",it) ; u component of wind on mass points
vv=wrf_user_getvar(a,"va",it) ; v component of wind on mass points
 
;uu = a->U10(it,:,:) ; u wind component at 10m
;vv = a->V10(it,:,:) ; v wind component at 10m
 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;maxlat=59.1125
;maxlon=-137.228
;minlat=58.3728
;minlon=-135.739
;lats =(/maxlat,minlat/) 
;lons = (/maxlon,minlon/)
;    loc = wrf_user_latlon_to_ij(a,lats,lons)
;           ys = loc(0,0)                        ; Starting y (j) location
;                  ye = loc(1,0)                        ; Ending y (j)
location
;                  xs = loc(0,1)                         ; Starting x (i)
location
;                  xe = loc(1,1)                        ; Ending x (i)
location
 
do ip = 3,3    ; we are doing 3 plots, specifying different start and end
points
  opts =  True           ; setting start and end times/points
  plane = new(4,float)
 
  if(ip .eq. 1) then
  plane = (/51,42, 52,62/) ; start x;y & end x;y point (X-sec top left
diagonal)
;The values should be ordered as (xstart,xend,ystart,yend). ; plane =
(/xs,ys,  xe,ye/) 
 end if
        if(ip .eq. 2) then
         plane = (/48,49, 50,66/) ; start x;y & end x;y point (X-sec top
right diagonal)
       end if
               if(ip .eq. 3) then
                                plane = (/43,55,   57,57/)
       end if
 
; Interpolate WRF Model data vertically 
  angle = 0. 
 
; Establish starting and end points in domain
;  plane = new(4,float)
;  plane = (/  xs,ys,  xe,ye  /) ; start x;y & end x;y point
 pm_plane = wrf_user_intrp3d(pm,z,"v",plane,angle,opts)   ; "v" vertical
cross section
tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts)
 u_plane = wrf_user_intrp3d(uu,z,"v",plane,angle,opts)
v_plane = wrf_user_intrp3d(vv,z,"v",plane,angle,opts)
 
u_plane=u_plane*1.94386 ; convert to kts
v_plane=v_plane*1.94386 
u_plane@units = "kts"
v_plane@units = "kts"
 
  xcoord = ispan(0,dimsizes(pm_plane(0,:)), 1) ;takes the size of the plane
(in this case 12 grid points)
  dim = dimsizes(pm_plane)  ; find the data span for use in labels
zspan = dim(0)
 
; Plotting Options for PM
  opts_xy = res
  opts_xy@tiYAxisString          = "Altitude (km)"
  opts_xy@tiXAxisString          = "Cross Section Distance (km)"
  opts_xy@tiXAxisFontHeightF      = 0.015
  opts_xy@tiYAxisOn              = True
  opts_xy@tmXBMode                = "Explicit" ;turn manual X axis labels on
  opts_xy@tmXBValues             =  xcoord     ;labels the start to end grid
points as stated above
  opts_xy@tmXBLabels            =  xcoord*7    ;take values and mult by 7
since grid cells are 7km width
  opts_xy@tmXBLabelFontHeightF = 0.01
  opts_xy@tmXTOn            = False  ;turns off top plot tickmarks
  opts_xy@tmYROn            = False  ;turns off right tickmarks
  opts_xy@gsnMaximize          = True         ;Plot maximized in workstation
 
; yxis options for Height/PM
  opts_xy@tmYLMode = "Explicit"
  opts_xy@tmYLValues = fspan(0,zspan,nz)
  opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz))
;  opts_xy@trYMaxF = 24  ; not sure how this works?
  ; Specifies a minimum Y coordinate value that defines the lower bound of
the Y-Axis coordinate system.
  
; PM10 Plotting Commands
  opts_xy@cnFillOn         = True  ;Fill in the PM Contours
  opts_xy@cnFillMode       = "AreaFill" ; Smooth them
  opts_xy@cnLinesOn          = False    ; turn off PM contour lines
  opts_xy@cnLineLabelsOn     = False    ; turn off PM contour labels
 
;options for size of plot dimensions
  opts_xy@vpXF = 0.105   ; left edge of the view objects bounding box
  opts_xy@vpWidthF=0.75  ; width of the viewing box
  opts_xy@vpYF = 0.84355 ; top edge location
  opts_xy@vpHeightF = 0.51 ; height of the viewing box
 
;--Temperature lotting Options
  opts_tc = res
  opts_tc@cnInfoLabelOn = False   ; turn off contour label
  opts_tc@ContourParameters  = (/ 2. /) ; label temp every 5 degrees
opts_tc@cnLineLabelInterval = 1 ; every contour labeled
opts_tc@cnLineLabelDensityF = .5  ; only one label per contour line
; opts_tc@gsnMaximize = True
;  opts_tc@tiXAxisOn = False
 
;options for size of plot dimensions NECESSARY
  opts_tc@vpXF = 0.105
  opts_tc@vpWidthF=.75
  opts_tc@vpYF = 0.84355
  opts_tc@vpHeightF = 0.51
 
;Same Y Plotting options to restrict data to certain level
  opts_tc@tmYLMode = "Explicit"
  opts_tc@tmYLValues = fspan(0,zspan,nz)
  opts_tc@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz))
;  opts_tc@trYMaxF = 24  ; cut off plot at this level?
 
 
 
;---Wind Plotting Options
opts_w=res
opts_w@FieldTitle = "Wind Speed"        ; Overwrite Field Title
opts_w@vcMinDistanceF        = 0.05      ; thin out windbarbs
opts_w@NumVectors =   20          ; Wind barb density
opts_w@gsnMaximize = True
 
;options for size of plot dimensions NECESSARY
  opts_w@vpXF = 0.105
  opts_w@vpWidthF=.75
  opts_w@vpYF = 0.84355
  opts_w@vpHeightF = 0.51
 
;Same Y Plotting options to restrict data to certain level
  opts_w@tmYLMode = "Explicit"
  opts_w@tmYLValues = fspan(0,zspan,35)
  opts_w@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,35))
; opts_w@trYMaxF = 24  ; cut off plot at this level?
 
; Contour information for PM, Temp, Wind
contour_pm = wrf_contour(a,wks,pm_plane,opts_xy)
contour_tc = wrf_contour(a,wks,tc_plane,opts_tc) 
 vector=wrf_vector(a,wks,u_plane,v_plane,opts_w)
 
; RESOURCES FOR FINAL OVERLAY PLOT
pltres@gsnMaximize = True
pltres@vpXF = 0.105
pltres@vpWidthF=0.75
pltres@vpYF = 0.84355
pltres@vpHeightF = 0.51
 
; Plot the graph
plot = wrf_overlays(a,wks,(/contour_pm,contour_tc,vector/),pltres)
 
; delete option fields to avoid carry over
delete(opts_tc)
delete(opts_xy)
delete(tc_plane)
delete(pm_plane)
 
end do  ; make the next cross section
;-----------------------------------------
end do ; end of time loop
;----------------------------------------
end


_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

image005.jpg image006.jpg
Received on Mon Dec 2 23:48:38 2013

This archive was generated by hypermail 2.1.8 : Wed Dec 04 2013 - 20:42:38 MST