[ncl-talk] Time vs. height (Pressure) plot not overlaying wind barbs
Scott Capps
scott at allvertum.com
Tue Dec 23 16:11:55 MST 2014
Happy holidays,
I am creating a time vs. height (in pressure coords) plot and am having
a problem with adding wind barbs. The wind barbs are not being plotted
while all other variables are and I am not getting any error messages.
Attached is my ncl script (time_height_plot.ncl) and associated script
(time_conversion_module.ncl). This script reads in post-processed WRF
output. I am using NCL version 6.2.1 running Ubuntu 64-bit OS with gcc
version 4.8.2. I can provide some sample upon request.
Thank you,
Scott
-------------- next part --------------
;***********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCL_FXN_PATH/time_conversion_module.ncl"
; USAGE:
; ncl 'stnlat=38.5556' 'stnlon=-121.4689' 'fil_inv="/home/scapps/projects/pg_e/wrf/wdt_domain/invariant_d01.nc"' 'dir_in="/data2/scapps/wrf_pge_nam_12z/"' 'fl_search_str="wrfout_d01_2014-12-13*.nc"' 'dir_out="./"' 'fl_out="timehgt_sacramento_plot"' 'fl_typ="ps"' time_height_plot.ncl
;***********************************************
begin
;***********************************************
; search for existing model output at this time
str_tmp = "ls "+dir_in+fl_search_str
all_files = systemfunc (str_tmp)
fall = addfiles (all_files, "r")
ListSetType (fall, "cat") ; concatenate or "merge" (default)
;
; INVARIANT FILE
file1 = addfile(fil_inv,"r")
lat2d = file1->XLAT(0,:,:)
lon2d = file1->XLONG(0,:,:)
lat2d at units= "degrees_north"
lon2d at units= "degrees_east"
;
; Find the ij location for the point if interest
llres = True
llres at ReturnInt = True ; Return integer values
locij = wrf_user_ll_to_ij(file1, stnlon, stnlat, llres)
locij = locij - 1 ; array pointers in NCL space
locX = locij(0)
locY = locij(1)
;
; Handle Times
offset = 0
tzone_offset_get(offset)
times = fall[:]->Times
tmp1 = chartostring(times(:,0:3))
tmp2 = chartostring(times(:,5:6))
tmp3 = chartostring(times(:,8:9))
tmp4 = chartostring(times(:,11:12))
;
utc_yyyy = stringtointeger(tmp1)
utc_mm = stringtointeger(tmp2)
utc_dd = stringtointeger(tmp3)
utc_hh = stringtointeger(tmp4)
;
delete(tmp1)
delete(tmp2)
delete(tmp3)
delete(tmp4)
ntimes = dimsizes(utc_yyyy)
lst_yy = new((/ntimes/),integer)
lst_mm = new((/ntimes/),integer)
lst_dd = new((/ntimes/),integer)
lst_hh = new((/ntimes/),integer)
utc_to_local_time(offset,utc_yyyy,utc_mm,utc_dd,utc_hh,lst_yy,lst_mm,lst_dd,lst_hh)
delete(utc_mm)
delete(utc_dd)
delete(utc_yyyy)
delete(utc_hh)
delete(times)
;
minxval = new(ntimes,integer)
lbl_count = 0
do i=0,ntimes-1
if (lst_hh(i).eq.0.or.lst_hh(i).eq.12.or.lst_hh(i).eq.18.or.lst_hh(i).eq.6) then
lbl_count = lbl_count+1
end if
end do
majxval = new(lbl_count,integer)
majxlbl = new(lbl_count,string)
idx = 0
do i=0,ntimes-1
minxval(i) = i
if (lst_hh(i).eq.0.or.lst_hh(i).eq.12.or.lst_hh(i).eq.18.or.lst_hh(i).eq.6) then
majxval(idx) = i
if (lst_hh(i).eq.0.or.lst_hh(i).eq.12) then
majxlbl(idx) = sprinti("%0.2i",lst_dd(i))+"/"+sprinti("%0.2i",lst_hh(i))
else
majxlbl(idx) = ""
end if
idx = idx+1
end if
end do
;
wks = gsn_open_wks(fl_typ,dir_out+fl_out)
colors = (/ (/255,255,255/), (/ 0, 0, 0/), (/255,255,255/), \
(/255,255,255/), (/ 80,255, 80/), (/ 50,200, 50/), (/ 20,150, 20/) /) / 255.0
gsn_define_colormap(wks,colors) ; create colormap out of above colors
;
;-----------------------------------------------------------------------
geohgt = fall[:]->GEOHGT ; [Time | 1] x [bottom_top_stag | 51] x [south_north | 320] x [west_east | 256]
dimW = dimsizes(geohgt(0,:,0,0))
tptk = fall[:]->T ; [K] Temperature
prsfull= fall[:]->P ; [Pa] Full model pressure
qwv = fall[:]->QVAPOR ; [kg/kg]
geohgt_unst = 0.5*(geohgt(:,:dimW-2,:,:)+geohgt(:,1:dimW-1,:,:))
;
wcomp = fall[:]->W
tmpomg = wrf_omega(qwv,tptk,wcomp,prsfull) ; [Pa/s]
delete(wcomp)
tmpomg = 10.*tmpomg ; [ub/s]
copy_VarAtts(tptk,tmpomg)
omega = tmpomg(bottom_top|0:30,Time|:,south_north|locY,west_east|locX)
delete(tmpomg)
; 4D RH
rh = wrf_rh(qwv,prsfull,tptk)
rhrh = rh(bottom_top|0:30,Time|:,south_north|locY,west_east|locX)
delete(rh)
;
; THETA-E (Equivalent Potential Temperature)
thetae = wrf_eth(qwv,tptk,prsfull)
thetae_out = thetae(bottom_top|0:30,Time|:,south_north|locY,west_east|locX)
delete(thetae)
; Get Pressure values for y-axis
prs_lvl = prsfull(bottom_top|0:30,Time|:,south_north|locY,west_east|locX)*0.01 ; Convert to MB
; Create 2D variable for plotting
vfarr = new((/31,ntimes/),integer)
do i=0,ntimes-1
do i2=0,30
vfarr(i2,i) = i
end do
end do
;
delete(geohgt)
delete(prsfull)
delete(qwv)
delete(geohgt_unst)
;
; 4D Winds - Earth-relative and unstaggered
umet = fall[:]->U ;m/s
vmet = fall[:]->V ;m/s
delete(fall)
;
ugrid = umet(bottom_top|0:30,Time|:,south_north|locY,west_east|locX)
delete(umet)
ugrid = ugrid*1.943846
vgrid = vmet(bottom_top|0:30,Time|:,south_north|locY,west_east|locX)
delete(vmet)
vgrid = vgrid*1.943846
delete(vgrid at coordinates)
delete(ugrid at coordinates)
delete(vgrid at stagger)
delete(ugrid at stagger)
;-----------------------------------------------------------------------
res2D = True
res2D at gsnMaximize = True
res2D at gsnDraw = False
res2D at gsnFrame = False
res2D at vpWidthF = 0.70
res2D at vpHeightF = 0.40
res2D at tiXAxisString = ""
res2D at tiXAxisFontHeightF = 0.016
res2D at trGridType = "TriangularMesh"
;
res2D at tmXBMode = "Explicit"
res2D at tmXBValues = majxval ;taus
res2D at tmXBLabels = majxlbl ;times
res2D at tmXBLabelJust = "CenterCenter"
res2D at tmXBLabelFontHeightF = .010
res2D at trYReverse = True ; Reverse the Y values.
res2D at tmXMajorGrid = True
res2D at tmXMajorGridThicknessF = 0.8
res2D at tmXMajorGridLineDashPattern = 1
res2D at tmXBMinorValues = minxval
;
;res2D at tiYAxisString = "Pressure (mb)"
res2D at tmYMajorGrid = True
res2D at tmYMajorGridThicknessF = 0.8
res2D at tmYMajorGridLineDashPattern = 1
res2D at tmYLMode = "Explicit"
res2D at tmYLValues = (/1000.,925.,850.,700.,500.,400.,300./)
res2D at tmYLLabels = (/1000.,925.,850.,700.,500.,400.,300./)
;
; Omega
omg_res = res2D
omg_res at sfXArray = vfarr
omg_res at sfYArray = prs_lvl
omg_res at cnLineLabelsOn = True
omg_res at cnLineColor = "Red" ; color of contour lines
omg_res at cnLineLabelFontColor = "Red"
omg_res at cnInfoLabelOn = False
omg_res at cnLineLabelBackgroundColor = -1
omg_res at cnLineLabelFontHeightF = .010
omg_res at cnLevelSelectionMode = "ManualLevels"
omg_res at cnMinLevelValF = -8.
omg_res at cnMaxLevelValF = 0.
omg_res at cnLevelSpacingF = 2.
;
tt_res = res2D
tt_res at sfXArray = vfarr
tt_res at sfYArray = prs_lvl
tt_res at cnLineLabelsOn = True ; no contour labels
tt_res at cnLineThicknessF = 3.0 ; line thickness
tt_res at cnLineColor = "Black" ; color of contour lines
tt_res at cnLineLabelFontColor = "Black"
tt_res at cnInfoLabelOn = False
tt_res at cnLineLabelBackgroundColor = -1
tt_res at cnLineLabelFontHeightF = .010
tt_res at cnLevelSelectionMode = "ManualLevels"
tt_res at cnMinLevelValF = 200.
tt_res at cnMaxLevelValF = 400.
tt_res at cnLevelSpacingF = 3.
rh_res = res2D
rh_res at sfXArray = vfarr
rh_res at sfYArray = prs_lvl
rh_res at gsnSpreadColors = True ; use full range of colors
rh_res at gsnSpreadColorEnd = -2
rh_res at cnLineLabelsOn = True;False
rh_res at cnLineLabelBackgroundColor = -1
rh_res at cnFillOn = True ; turns on color fill
rh_res at cnInfoLabelOn = False
rh_res at cnLevelSelectionMode = "ManualLevels"
rh_res at cnMinLevelValF = 50.
rh_res at cnMaxLevelValF = 100.
rh_res at cnLevelSpacingF = 10
rh_res at cnLineLabelFontHeightF = .010
rh_res at cnLineColor = "DarkGreen" ; color of contour lines
rh_res at cnLineLabelFontColor = "DarkGreen"
;
uv_res = res2D
uv_res at vfXArray = vfarr
uv_res at vfYArray = prs_lvl
uv_res at vcRefMagnitudeF = 10.
uv_res at vcRefAnnoOn = False ; turns off the ref vector
uv_res at vcRefLengthF = 0.04 ; set length of ref vector
uv_res at vcGlyphStyle = "WindBarb" ; turn on wind barbs
uv_res at vcMinDistanceF = 0.002
uv_res at vcMapDirection = False
uv_res at vcMonoWindBarbColor = True
;-----------------------------------------------------------------------
;-----------------------------------------------------------------------
rhfill = gsn_contour(wks,rhrh,rh_res)
ttfill = gsn_contour(wks,thetae_out,tt_res)
overlay(rhfill,ttfill)
omgfill = gsn_contour(wks,omega,omg_res)
overlay(rhfill,omgfill)
windlayer = gsn_vector(wks,ugrid,vgrid,uv_res)
overlay(rhfill,windlayer)
draw(rhfill)
frame(wks)
;-----------------------------------------------------------------------
end
-------------- next part --------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.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"
procedure tzone_offset_get(offset)
begin
tmp = systemfunc("date '+%z' | awk '{print substr($0,0,3)}'")
if (strlen(tmp).lt.3) then
tmp = systemfunc("date '+%z' | awk '{print substr($0,0,4)}'")
end if
offset = stringtoint(tmp)
end
procedure utc_to_local_time(local_time_offset,utc_yyyy,utc_mm,utc_dd,utc_hh,yyyy,mm,dd,hh)
begin
; local_time_offset is passed in as: -07 or -08, etc..
; from linux date command: date +%Z
; timezone_offset=`date +%z`
; timezone=`date +%Z`
; offset=${timezone_offset:0:3}
; ncl 'offset="'${offset}'"' 'yyyy=2012' 'mm=1' 'dd=01' 'hh=0' test.ncl
; OR, from another NCL script:
; offset = systemfunc("date '+%z' | awk '{print substr($0,0,3)}'")
; ncl 'offset=offset' 'yyyy=2012' 'mm=1' 'dd=01' 'hh=0' test.ncl
; *****************************************************************
ltime_offset = local_time_offset
;
nidx = dimsizes(utc_yyyy)
yyyy = new((/nidx/),integer)
mm = new((/nidx/),integer)
dd = new((/nidx/),integer)
hh = new((/nidx/),integer)
;
; Check to see if hours from 0 to 23
if (max(utc_hh).gt.23)then
print("ERROR: Hours must range from 0 to 23, "+utc_hh(idx)+"...")
exit
end if
do idx=0,nidx-1
if (utc_hh(idx).ge.0.and.utc_hh(idx).le.abs(ltime_offset)-1)then
; Decrement Day and, maybe MM, YYYY
; Determine LST hour
hh(idx) = utc_hh(idx) + (24+ltime_offset)
if (utc_dd(idx).eq.1) then
; Decrement Month
if (utc_mm(idx).ge.2) then
mm(idx) = utc_mm(idx)-1
dd(idx) = days_in_month(utc_yyyy(idx),mm(idx))
yyyy(idx) = utc_yyyy(idx)
else
mm(idx) = 12
dd(idx) = days_in_month(utc_yyyy(idx),mm(idx))
yyyy(idx) = utc_yyyy(idx)-1
end if
else
dd(idx) = utc_dd(idx)-1
mm(idx) = utc_mm(idx)
yyyy(idx) = utc_yyyy(idx)
end if
else
hh(idx) = utc_hh(idx)+ltime_offset
dd(idx) = utc_dd(idx)
mm(idx) = utc_mm(idx)
yyyy(idx) = utc_yyyy(idx)
end if
end do
end
procedure date_hr_add(hr_add_in,utc_yyyy,utc_mm,utc_dd,utc_hh,yyyy,mm,dd,hh)
begin
yyyy = 0
mm = 0
dd = 0
hh = 0
; Hours added must be between 0 and 23
hr_add = stringtoint(hr_add_in)
;
; Check to see if hours from 0 to 23
if (utc_hh.gt.23)then
print("ERROR: Hours must range from 0 to 23...")
exit
end if
;
if (hr_add.gt.23)then
print("ERROR: Hours added must range from 0 to 23...")
exit
end if
dd_in_mm = days_in_month(utc_yyyy,utc_mm)
if (utc_hh.lt.(24-hr_add))then
hh = utc_hh+hr_add
dd = utc_dd
mm = utc_mm
yyyy = utc_yyyy
else
hh = hr_add-(24-utc_hh)
if (utc_dd.eq.dd_in_mm) then
dd = 1
if (utc_mm.eq.12) then
mm = 1
yyyy = utc_yyyy+1
else
mm = utc_mm+1
yyyy = utc_yyyy
end if
else
dd = utc_dd+1
mm = utc_mm
yyyy = utc_yyyy
end if
end if
end
procedure date_day_add(utc_yyyy,utc_mm,utc_dd,utc_hh,yyyy,mm,dd,hh)
begin
; Add one day
yyyy = 0
mm = 0
dd = 0
hh = 0
; Check to see if hours from 0 to 23
if (utc_hh.gt.23)then
print("ERROR: Hours must range from 0 to 23...")
exit
end if
dd_in_mm = days_in_month(utc_yyyy,utc_mm)
if (utc_dd.eq.dd_in_mm)then
hh = utc_hh
dd = 1
if (utc_mm.eq.12) then
mm = 1
yyyy = utc_yyyy+1
else
mm = utc_mm+1
yyyy = utc_yyyy
end if
else
hh = utc_hh
dd = utc_dd+1
mm = utc_mm
yyyy = utc_yyyy
end if
end
More information about the ncl-talk
mailing list