[ncl-talk] two color bars, different color tables
Matt Flournoy
mflournoy37 at gmail.com
Tue Jun 20 11:34:31 MDT 2017
Hi all,
I'm trying to plot some WRF output with surface reflectivity shading
overlaid with parcel trajectories. I'd like the reflectivity to be shaded
in a grayscale, and the parcel trajectories to be colored according to
their height. I currently have a color bar for the shaded reflectivity, but
am having trouble getting a second color bar to show the parcel heights. I
think my issue might come from my use of a different color table for the
trajectories than for the reflectivity, which is what I defined my
workstation colormap with (MPL_Greys).
I've attached the script I'm currently using and the resulting plot. I'm
trying to get the height color bar (using BlAqGrYeOrReVi200) on the right
side of the plot. When I run the script as is, I get this message in
between my "error1" and "error2" print lines:
fatal:CvtStringGenArrayToColorIndexGenArray: Unable to convert string
"BlAqGrYeOrReVi200" to requested type
warning:Error retrieving resource lbFillColors from args - Ignoring Arg
If you remove the semicolon before the "frame(xwks2)" line, the plot will
show the parcel height color bar (without the labels working...), but not
the rest of the plot.
Any help would be greatly appreciated!
Thank you,
Matt
--
Matthew Flournoy
M.S. Meteorology Candidate, University of Oklahoma
B.S. Meteorology, Penn State University[image: Inline image 1]
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170620/5c88fa83/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: traj_refl_back.png
Type: image/png
Size: 51707 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170620/5c88fa83/attachment.png
-------------- 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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;*******************************************************************************************************
; Script to read in WRF trajectory output and plot trajectories overlaid with lowest model level
; reflectivity.
;*******************************************************************************************************
begin
print("reading in trajectory text files...")
num_trajs = 2 ; <-- set number of trajectories to read in here
num_tsteps = 121 ; <-- enter number of timesteps here
;----------------------------------------------------------------------;
; Allocate trajectory variable arrays here.
;----------------------------------------------------------------------;
traj_time = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_lat = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_lon = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_x = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_y = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_z = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_u = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_v = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_w = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_pres = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
traj_refl = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
; traj_thp = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
; traj_qv = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
; traj_qc = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
; traj_qr = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
; traj_qi = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
; traj_qs = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
; traj_qg = new( (/num_trajs,num_tsteps/) , "float" , -999.0)
start = 0
count_index = 0
flow = "trajs_backwards"
do m = start , start + num_trajs - 1
index = (m) + 1
if (index .lt. 10) then
traj_file = "/home/matt.flournoy/wrf_trajectories/" + flow + "/trajwrf00" + index + ".out"
else if (index .ge. 10 .and. index .lt. 100) then
traj_file = "/home/matt.flournoy/wrf_trajectories/" + flow + "/trajwrf0" + index + ".out"
else if (index .ge. 100) then
traj_file = "/home/matt.flournoy/wrf_trajectories/" + flow + "/trajwrf" + index + ".out"
end if
end if
end if
;----------------------------------------------------------------------;
; Read in all lines in the file as strings, then read them into
; variable arrays by column.
;----------------------------------------------------------------------;
lines = asciiread(traj_file,-1,"string")
header = lines(0)
data = lines(1:)
delim = " "
traj_time(count_index,:) = tofloat(str_get_field(data,1,delim))
traj_lat(count_index,:) = tofloat(str_get_field(data,2,delim))
traj_lon(count_index,:) = tofloat(str_get_field(data,3,delim))
traj_x(count_index,:) = tofloat(str_get_field(data,4,delim))
traj_y(count_index,:) = tofloat(str_get_field(data,5,delim))
traj_z(count_index,:) = tofloat(str_get_field(data,6,delim))
traj_u(count_index,:) = tofloat(str_get_field(data,7,delim))
traj_v(count_index,:) = tofloat(str_get_field(data,8,delim))
traj_w(count_index,:) = tofloat(str_get_field(data,9,delim))
traj_pres(count_index,:) = tofloat(str_get_field(data,10,delim))
traj_refl(count_index,:) = tofloat(str_get_field(data,13,delim))
; traj_thp(m,:) = tofloat(str_get_field(data,14,delim))
; traj_qv(m,:) = tofloat(str_get_field(data,15,delim))
; traj_qc(m,:) = tofloat(str_get_field(data,16,delim))
; traj_qr(m,:) = tofloat(str_get_field(data,17,delim))
; traj_qi(m,:) = tofloat(str_get_field(data,18,delim))
; traj_qs(m,:) = tofloat(str_get_field(data,19,delim))
; traj_qg(m,:) = tofloat(str_get_field(data,20,delim))
delete(lines)
delete(header)
delete(data)
delete(delim)
print("done reading in " + count_index + "/" + num_trajs + " files...")
count_index = count_index + 1
end do
print("done reading in all trajectory data!")
traj_lat at units = "degrees north"
traj_lon at units = "degrees east"
traj_z at units = "m"
print("all trajectory data read in")
;----------------------------------------------------------------------;
;----------------------------------------------------------------------;
;----------------------------------------------------------------------;
; Now read in reflectivity from the WRF output file. In this case,
; I'll use the reflectivity from the center of the times that the
; trajectories were computed over.
;
; Trajectories were computed from 0430 - 0435 UTC.
; Surface reflectivity is plotted at 0432:30 UTC.
;
; NOTE FOR SANITY: All variables related to the WRF output will be
; preceeded by "wrf". All trajectory variables will
; be proceeded by "traj".
;----------------------------------------------------------------------;
print("now reading in reflectivity from WRF NetCDF output file...")
nc_filepath = "/work/michael.coniglio//PECAN/wrf-regrid/mem30restart/wrfcart_20150706-043300.nc" ; 042800.nc"
nc_file = addfile(nc_filepath,"r")
;--------------------------------------------------
; Read in the variables...
;--------------------------------------------------
wrf_refl = nc_file->REFLECTIVITY ; (z,y,x) <-- reflectivity (dBZ)
wrf_lats = nc_file->LATITUDE ; (y,x)
wrf_lons = nc_file->LONGITUDE ; (y,x)
wrf_hgts = nc_file->HEIGHT ; (z) <-- height above ground level (m)
; wrf_latgrid(:,:) = wrf_user_getvar(nc_file,"lat",0) ; Model latitude
; wrf_longrid(:,:) = wrf_user_getvar(nc_file,"lon",0) ; Model longriditude
;--------------------------------------------------
; Assign the lat/lon grids to "wrf_refl".
;--------------------------------------------------
wrf_refl at lat2d = wrf_lats
wrf_refl at lon2d = wrf_lons
print("WRF variables read in... now plot!")
;----------------------------------------------------------------------;
;----------------------------------------------------------------------;
;----------------------------------------------------------------------;
; Now plot! Produce a 3-D plot with surface reflectivity overlaid with
; the trajectory track.
;----------------------------------------------------------------------;
plotdir = "/home/matt.flournoy/wrf_trajectories/plots/traj_back/"
plotfile = plotdir + "/traj_refl_back"
plotformat = "png"
colors = "MPL_Greys" ; <-- make reflectivity gray-scale so trajectories can be in color
xwks = gsn_open_wks(plotformat,plotfile)
gsn_define_colormap(xwks,colors)
res = True
res at gsnDraw = False
res at gsnFrame = False
;---------------------------------------------
; Map projection stuff.
;---------------------------------------------
res at mpProjection = "Mercator"
res at mpLimitMode = "Corners"
clat = 43.531399 ; center of the plot grid
clon = -97.875127 ; center of the plot grid
plot_width = 0.30 ; 0.20 ; degrees longitude (5 /8/30 for 1/3/15 km grid) NOTE: Keep the aspect ratio (width/height) 1.666!
plot_height = 0.18 ; 0.12 ; degrees latitude (2.5/4/18 for 1/3/15 km grid)
lat_ll = clat - 0.5*plot_height
lat_ur = clat + 0.5*plot_height
lon_ll = clon - 0.5*plot_width
lon_ur = clon + 0.5*plot_width
res at mpLeftCornerLatF = lat_ll
res at mpLeftCornerLonF = lon_ll
res at mpRightCornerLatF = lat_ur
res at mpRightCornerLonF = lon_ur
res at mpShapeMode = "FreeAspect"
res at mpOutlineOn = True
res at mpGridAndLimbOn = False
res at mpPerimOn = True
res at mpFillOn = False
res at mpPerimDrawOrder = "PostDraw"
res at mpGeophysicalLineColor = "Black"
res at mpNationalLineColor = "Black"
res at mpUSStateLineColor = "Black"
res at mpGridLineColor = "Black"
res at mpLimbLineColor = "Black"
res at mpPerimLineColor = "Black"
res at mpGeophysicalLineThicknessF = 1.0
res at mpGridLineThicknessF = 1.0
res at mpLimbLineThicknessF = 1.0
res at mpNationalLineThicknessF = 0.5
res at mpUSStateLineThicknessF = 1.0
res at mpDataResolution = "FinestResolution"
res at mpOutlineBoundarySets = "AllBoundaries"
res at mpUSStateLineThicknessF = 4.0
res at mpDataBaseVersion = "Ncarg4_1"
res at mpDataSetName = "Earth..2"
res at mpCountyLineThicknessF = 0.5
;---------------------------------------------
; Title and text.
;---------------------------------------------
res at tiMainString = "Lagrangian Trajectories ~C~~Z70~Member 30 (0423 - 0433 UTC)"
res at tiMainJust = "CenterLeft"
res at tiMainPosition = "Left"
;---------------------------------------------
; Plotting specifications and COLORS.
;---------------------------------------------
res at cnFillMode = "AreaFill"
res at cnLinesOn = False
res at cnFillOn = True
res at cnLineLabelsOn = False
res at cnFillOpacityF = 1.0 ; <-- no transparency
res at cnLevelSelectionMode = "ManualLevels"
res at cnMinLevelValF = 0.0
res at cnMaxLevelValF = 70.0
res at cnLevelSpacingF = 5.0
res at cnFillColors = (/5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80/) ; <-- got these numbers from the MPL_Greys site since that's the color bar I'm using
;---------------------------------------------
; Label bar resources (manual).
;---------------------------------------------
res at lbLabelBarOn = True
res at pmLabelBarSide = "Bottom"
res at lbAutoManage = False ; we control label bar
res at pmLabelBarDisplayMode = "Always" ; turns on label bar
res at lbOrientation = "Horizontal" ; ncl default is vertical
res at pmLabelBarSide = "Bottom" ; default is right
res at lbLabelStride = 2 ; skip every other label
res at pmLabelBarWidthF = 0.4 ; default is shorter
res at pmLabelBarHeightF = 0.1 ; default is taller
res at lbLabelFontHeightF = .018 ; default is HUGE
res at lbPerimOn = False ; default has box
;---------------------------------------------
; Plot the 8 m AGL reflectivity.
;---------------------------------------------
plot = gsn_csm_contour_map(xwks,wrf_refl(0,:,:),res)
;---------------------------------------------
; Now overlay it with the parcel trajectory,
; based on the lat/lon points.
;---------------------------------------------
cmap = "BlAqGrYeOrReVi200"
print("done plotting reflectivity, now overlaying the parcel swarm...")
min_z = 50.0 ; <-- 50 m
max_z = max(traj_z(:,:)) ; <-- max height of any of the trajectories (m)
print("Maximum height reached by any parcel is " + max_z)
levels = 197 ; <-- max amount allowed by color table I'm using (200 colors)
cnLevels = fspan(min_z,max_z,levels)
nvalues = dimsizes(traj_lat(0,:))-1 ; <-- just get the number of timesteps in each trajectory
dum_line = new( (/num_trajs,nvalues/) , graphic )
pres = True
pres at gsLineThicknessF = 2.0
do t = 0 , num_trajs - 1 ; <-- loop over the trajectories
traj_num = t + 1
do i = 0 , dimsizes(traj_lat(0,:)) - 2 ; <-- loop over the lat/lon points
pres at gsLineColor = GetFillColor(cnLevels,cmap,avg( (/traj_z(t,i),traj_z(t,i+1)/) ))
dum_line(t,i) = gsn_add_polyline(xwks,plot,(/traj_lon(t,i),traj_lon(t,i+1)/),(/traj_lat(t,i),traj_lat(t,i+1)/),pres)
end do
print("done plotting " + traj_num + "/" + num_trajs + " trajectories")
end do
;---------------------------------------------
; Draw the initial spot of the trajectory
; with a marker.
;---------------------------------------------
print("Drawing initial location of each trajectory")
first = True
first at gsMarkerSizeF = 3.0
first at gsMarkerColor = "black"
first at gsMarkerIndex = 16
dum_marker = new(num_trajs,graphic)
do n = 0 , num_trajs - 1
dum_marker(n) = gsn_add_polymarker(xwks,plot,traj_lon(n,0),traj_lat(n,0),first)
end do
;---------------------------------------------
; Now add the label bar...
;---------------------------------------------
xwks2 = gsn_open_wks(plotformat,plotfile)
gsn_define_colormap(xwks2,cmap)
print("lastly, adding the height color bar...")
res_lb = True
res_lb at lbPerimOn = False
res_lb at vpWidthF = 0.1
res_lb at vpHeightF = 0.5
res_lb at lbOrientation = "Vertical"
res_lb at lbFillColors = cmap
res_lb at lbLabelFontHeightF = 0.015
res_lb at lbMonoFillPattern = True ; Fill them all solid.
res_lb at lbMonoFillColor = False ; Fill with different colors .
res_lb at lbLabelAlignment = "InteriorEdges" ; Default is "BoxCenters".
res_lb at lbBoxLinesOn = False ; Turn off labelbar box lines
LabelStrings = (/"50","75","100","125","150","175","200"/)
print("error1")
gsn_labelbar_ndc(xwks2,dimsizes(cnLevels)+1,LabelStrings,0.83,0.78,res_lb)
print("error2")
draw(plot)
frame(xwks)
; frame(xwks2)
print("done plotting!")
end
More information about the ncl-talk
mailing list