[ncl-talk] plotting a polygon over an existing plot
Alison Bridger
alison.bridger at sjsu.edu
Tue Jul 31 17:02:26 MDT 2018
Ooops - wtest.ncl attached!
Another way to state my problem is by looking at the sample code under
gsn_add_polygon
The quantities dum1, dum2, and dum3 are defined (left of equals sign) but I
don't see them again in the code, so how are they used???
On Tue, Jul 31, 2018 at 2:06 PM Alison Bridger <alison.bridger at sjsu.edu>
wrote:
> Although my NCL skills are gradually improving, I am still stumped by
> simple tasks!
>
> The attached code reads in WRF output, extracts rainfall, sums it over the
> run, and contour plots the total each hour. I got it from somebody else
> here!
>
> I want to superimpose ("overplot") small columns at a few locations
> showing obs versus WRF rain totals. So imagine - say - small red and blue
> columns side-by-side somewhere on this plot (e.g., SFO) showing obs and WRF
> rain totals.
>
> I have spent 2 hours (!) trying to morph the sample code in polyg_4.ncl
> (from the NCL pages) into my code (wtest.ncl, attached), but I get nothing
> visual and no useful diagnostics.
>
> I fell this is related to things like "frame" and "draw" and "plot" and so
> forth, but no amount of trying seems to make me understand how these relate
> to each other!
>
> The added code starts near the end with "xpts = ...".
>
> Thanks for any pointers!
>
> Alison
>
> PS can supply wrfout file if needed
>
> --
>
> Alison F.C. Bridger
> Professor & Chair
>
> Department of Meteorology and Climate Science
>
> San Jose State University tel 408.924.5206
> One Washington Square fax 408.924.5191
> San Jose, CA 95192-0104
>
> email: Alison.Bridger at sjsu.edu
>
> *Global CO2 levels...410 ppm and still rising*
>
>
> www.sjsu.edu/meteorology <http://www.met.sjsu.edu>
>
>
--
Alison F.C. Bridger
Professor & Chair
Department of Meteorology and Climate Science
San Jose State University tel 408.924.5206
One Washington Square fax 408.924.5191
San Jose, CA 95192-0104
email: Alison.Bridger at sjsu.edu
*Global CO2 levels...410 ppm and still rising*
www.sjsu.edu/meteorology <http://www.met.sjsu.edu>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180731/7b3ceab6/attachment.html>
-------------- next part --------------
; Example script to produce plots for a WRF real-data run,
; with the ARW coordinate dynamics option.
; Edited by moi to provide a closeup on the Bay Area!
; We only plot sigma(rain)
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
; a = addfile("../wrfout_d01_2000-01-24_12:00:00.nc","r")
; domain 01
a = addfile("wrfout_d02_2016-10-28_00:00:00_45lcu0mp8d02","r")
; a = addfile("/home/008143031/sandbox/wrfout/v39/wrfout215gcu5mp815s_161028gfs","r")
; domain 02
; a = addfile("wrfout_d01_2016-10-28_00:00:00.nc","r")
; We generate plots, but what kind do we prefer?
type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"Precip_10-30-16-45levels-cu1mp8-dom02")
; Set some basic resources
res = True
res at MainTitle = "REAL-TIME WRF"
pltres = True
mpres = True
mpres at mpGeophysicalLineColor = "Black"
mpres at mpNationalLineColor = "Black"
mpres at mpUSStateLineColor = "Black"
mpres at mpGridLineColor = "Black"
mpres at mpLimbLineColor = "Black"
mpres at mpPerimLineColor = "Black"
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What times and how many time steps are in the data set?
FirstTime = True
times = wrf_user_getvar(a,"times",-1) ; get all times in the file
ntimes = dimsizes(times) ; number of times in the file
; print("ntimes:" + ntimes )
;i;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; do it = 0,ntimes-1,2 ; TIME LOOP
; do it = 0,ntimes-1 ; TIME LOOP
; do it = 7,31 ; TIME LOOP
; do it = 0,31 ; TIME LOOP
do it = 0,2 ; TIME LOOP
print("Working on time: " + times(it) )
if (FirstTime) then ; Save some times for tracking tendencies
times_sav = times(it)
end if
res at TimeLabel = times(it) ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
slp = wrf_user_getvar(a,"slp",it) ; slp
wrf_smooth_2d( slp, 3 ) ; smooth slp
; Prep for zoom map, per Frank
lat = a->XLAT(it,:,:)
lon = a->XLONG(it,:,:)
dims = dimsizes(lat)
nlat = dims(0)
nlon = dims(1)
latS = 36.83
latN = 38.83
lonW = -123.58
lonE = -121.48
; latS = 35.
; latN = 40.
; lonW = -126.
; lonE = -120.
ji = region_ind(lat, lon, latS, latN, lonW, lonE)
jStrt = ji(0)
jLast = ji(1)
iStrt = ji(2)
iLast = ji(3)
latz = lat(jStrt:jLast , iStrt:iLast)
lonz = lon(jStrt:jLast , iStrt:iLast)
slpz = slp(jStrt:jLast , iStrt:iLast)
dimz = dimsizes(latz)
nlatz = dimz(0)
nlonz = dimz(1)
; data gathering and printing
latSJ = 37.3639
lonSJ = -121.9289
latOK = 37.7126
lonOK = -122.2197
latSF = 37.6213
lonSF = -122.3790
latSC = 36.9741
lonSC = -122.0308
llres = True
llres at ReturnInt = True
locij = wrf_user_ll_to_ij(a,lonSJ,latSJ,llres)
locij = locij - 1
locXSJ = locij(0)
locYSJ = locij(1)
locij = wrf_user_ll_to_ij(a,lonOK,latOK,llres)
locij = locij - 1
locXOK = locij(0)
locYOK = locij(1)
locij = wrf_user_ll_to_ij(a,lonSF,latSF,llres)
locij = locij - 1
locXSF = locij(0)
locYSF = locij(1)
locij = wrf_user_ll_to_ij(a,lonSC,latSC,llres)
locij = locij - 1
locXSC = locij(0)
locYSC = locij(1)
; print("San Jose X grid point:" + locXSJ )
; print("San Jose Y grid point:" + locYSJ )
; print("Oakland X grid point:" + locXOK )
; print("Oakland Y grid point:" + locYOK )
; print("SFO X grid point:" + locXSF )
; print("SFO Y grid point:" + locYSF )
; print("Santa Cruz X grid point:" + locXSC )
; print("Santa Cruz Y grid point:" + locYSC )
; Get non-convective, convective and total precipitation
; Calculate tendency values
rain_exp = wrf_user_getvar(a,"RAINNC",it)
rain_con = wrf_user_getvar(a,"RAINC",it)
; print(rain_con(locYSJ,locXSJ))
rain_tot = rain_exp + rain_con
; grab sigma(rain) after 7h = midnight
; then subtract to get midnight-midnight rain
; rain_tot7 = wrf_user_getvar(a,"RAINNC",7)+wrf_user_getvar(a,"RAINC",7)
; rain_tot = rain_tot - rain_tot7
rain_tot at description = "Total Precipitation"
if( FirstTime ) then
if ( it .eq. 0 ) then
rain_exp_save = rain_exp
rain_con_save = rain_con
rain_tot_save = rain_tot
else
rain_exp_save = wrf_user_getvar(a,"RAINNC",it-1)
rain_con_save = wrf_user_getvar(a,"RAINC",it-1)
rain_tot_save = rain_exp_save + rain_con_save
FirstTime = False
times_sav = times(it-1)
end if
end if
rain_exp_tend = rain_exp - rain_exp_save
rain_con_tend = rain_con - rain_con_save
rain_tot_tend = rain_tot - rain_tot_save
rain_exp_tend at description = "Explicit Precipitation Tendency"
rain_con_tend at description = "Param Precipitation Tendency"
rain_tot_tend at description = "Precipitation Tendency"
; Bookkeeping, just to allow the tendency at the next time step
rain_exp_save = rain_exp
rain_con_save = rain_con
rain_tot_save = rain_tot
rain_totz = rain_tot(jStrt:jLast , iStrt:iLast)
; convert to bloody inches!
rain_totz = rain_totz / 25.4
; grab totals at SJ and Oakland only
rtzSJ = rain_tot_save(locYSJ,locXSJ) / 25.4
; print("San Jose sigma rain (in): " + rtzSJ)
rtzOK = rain_tot_save(locYOK,locXOK) / 25.4
; print("Oakland sigma rain (in): " + rtzOK)
rtzSF = rain_tot_save(locYSF,locXSF) / 25.4
; print("SFO sigma rain (in): " + rtzSF)
rtzSC = rain_tot_save(locYSC,locXSC) / 25.4
; print("Santa Cruz sigma rain (in): " + rtzSC)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if( .not. FirstTime ) then ; We will skip the first time
; Plotting options for Sea Level Pressure
; opts_psl = res
; opts_psl at ContourParameters = (/ 900., 1100., 2. /)
; opts_psl at cnLineColor = "Blue"
; opts_psl at cnInfoLabelOn = False
; opts_psl at cnLineLabelFontHeightF = 0.01
; opts_psl at cnLineLabelPerimOn = False
; opts_psl at gsnContourLineThicknessesScale = 1.5
; opts_psl at gsnMaximize = False
; contour_psl = wrf_contour(a,wks,slp,opts_psl)
; delete(opts_psl)
; Plotting options for Precipitation
opts_r = res
opts_r at UnitLabel = "in"
opts_r at cnLevelSelectionMode = "ExplicitLevels"
; millimeters
; opts_r at cnLevels = (/ .1, .2, .4, .8, 1.6, 3.2, 6.4, \
; 12.8, 25.6, 51.2, 102.4/)
; opts_r at cnFillColors = (/"White","White","DarkOliveGreen1", \
; "DarkOliveGreen3","Chartreuse", \
; "Chartreuse3","Green","ForestGreen", \
; "Yellow","Orange","Red","Violet"/)
; inches
opts_r at cnLevels = (/ .05, .1, .2, .3, .4, .5, .6, \
.75, 1., 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, \
2.75, 3.0, 3.25, 3.5 /)
opts_r at cnFillColors = (/"White","DarkOliveGreen1","DarkOliveGreen2", \
"DarkOliveGreen3","Chartreuse", \
"Chartreuse3","Green","ForestGreen", \
"Yellow","Orange","Red","Violet", \
"CadetBlue1","Blue","Blue4","Blueviolet", \
"sienna","sienna3","sienna4","tan4"/)
; opts_r at cnLevels = (/ .05, .1, .2, .3, .4, .5, .6, \
; .75, 1., 1.25, 1.5/)
; opts_r at cnFillColors = (/"White","DarkOliveGreen1","DarkOliveGreen2", \
; "DarkOliveGreen3","Chartreuse", \
; "Chartreuse3","Green","ForestGreen", \
; "Yellow","Orange","Red","Violet"/)
opts_r at cnInfoLabelOn = False
opts_r at cnConstFLabelOn = False
opts_r at cnFillOn = True
; Zooming
mpres at mpLimitMode = "Corners"
mpres at mpLeftCornerLatF = latS
mpres at mpLeftCornerLonF = lonW
mpres at mpRightCornerLatF = latN
mpres at mpRightCornerLonF = lonE
; mpres at mpLeftCornerLatF = latz(0,0)
; mpres at mpRightCornerLonF = lonz(nlatz-1,nlonz-1)
; Total Precipitation (color fill)
contour_tot = wrf_contour(a,wks, rain_totz, opts_r)
; add a random box
xpts = (/ 37.0, 37.2, 37.2, 37.0, 37.0/)
ypts = (/-123.0,-123.0,-122.8,-122.8,-123.0/)
;
resp = True
resp at gsLineColor = "red"
resp at gsLineThicknessF = 2.0
resp at gsLineLabelString = "test"
dum = new(4,graphic)
do i=0,3
dum(i) = gsn_add_polyline(wks,contour_tot,xpts(i:i+1),ypts(i:i+1),resp)
end do
;
; dum2 = gsn_add_polygon(wks,contour_tot,xpts,ypts,resp)
; MAKE PLOTS
; Total Precipitation
plot = wrf_map_overlays(a,wks,contour_tot,pltres,mpres)
draw(plot)
; frame(wks)
;
end if ; END IF FOR SKIPPING FIRST TIME
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
times_sav = times(it)
FirstTime = False
end do ; END OF TIME LOOP
end
More information about the ncl-talk
mailing list