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
