[ncl-talk] Shapefile does not attach to plot
zilore mumba
zmumba at yahoo.com
Mon Apr 22 12:26:07 MDT 2019
I hope someone can assist me where to put my shapefile drawing function so that it be attached to the plot. Currently I get my plot followed by the shapefile plotted separately.I have ready the solution by Maxine to this here [ncl-talk] I can't manage to create a plot with multiple lines on it, each line having different x and y values. Attribute assignment type issue, were he solved his proble by: either both "delete(res)" and "delete(wks)" before each plot
- or a new "res" (like "res3" for a 3rd plot) and a new "wks" (like "wks3" for a 3rd plot) forI have tried these suggestions, but really not clear how to do it.The wrf_pressureLevel1.ncl script I use is given below.
Help will be appreciated;:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::;;
load "./Shapefiles/shapefile_utils.ncl"
begin
;---Open WRF output file.
dir = "./"
filename = "wrfout_d01_2019-03-28_00:00:00"
a = addfile(dir + filename,"r")
; We generate plots, but what kind do we prefer?
type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"plt_PressureLevel1")
; Set some Basic Plot options
res = True
res at gsnDraw = False ; don't draw
res at gsnFrame = False ; don't advance frame
; res at MainTitle = "REAL-TIME WRF"
res at gsnMaximize = True ; Maximize plot in frame.
pltres = True
mpres = True ; Plot options desired.
mpres at mpGeophysicalLineColor = "Black" mpres at mpNationalLineColor = "Black"
mpres at mpInlandWaterFillColor = "Blue"
mpres at mpGridLineColor = "Black"
mpres at mpLimbLineColor = "Black"
mpres at mpPerimLineColor = "Black"
mpres at mpGeophysicalLineThicknessF = 2.0
mpres at mpGridLineThicknessF = 2.0
mpres at mpDataBaseVersion = "MediumRes" ; use finer database
mpres at mpDataSetName = "Earth..4"
mpres at mpOutlineBoundarySets = "AllBoundaries"
mpres at mpOutlineOn = True
mpres at mpNationalLineThicknessF = 3.0 ; turn on country boundaries
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; define shapefiles
shp_filename = ("Shapefiles/Zambia/ZMB_adm1.shp")
; set shapefile resources
shpres = True
shpres at gsLineThicknessF = 3.9 ; increase line thickness
shpres at gsLineColor = "Purple" ; line colorgsLineThicknessF
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What times and how many time steps are in the data set?
times = wrf_user_getvar(a,"times",-1) ; get all times in the file
ntimes = dimsizes(times) ; number of times in the file
; The specific pressure levels that we want the data interpolated to.
pressure_levels = (/ 850., 700., 500., 300./) ; pressure levels to plot
nlevels = dimsizes(pressure_levels) ; number of pressure levels
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
do it = 0,ntimes-1,2 ; TIME LOOP
print("Working on time: " + times(it) )
res at 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
u = wrf_user_getvar(a,"ua",it) ; u averaged to mass points
v = wrf_user_getvar(a,"va",it) ; v averaged to mass points
p = wrf_user_getvar(a, "pressure",it) ; pressure is our vertical coordinate
z = wrf_user_getvar(a, "z",it) ; grid point height
rh = wrf_user_getvar(a,"rh",it) ; relative humidity
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; printVarSummary(t)
wks = gsn_open_wks("x11","contour_wrf") ; Open X11 window
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
do level = 0,nlevels-1 ; LOOP OVER LEVELS
pressure = pressure_levels(level)
tc_plane = wrf_user_intrp3d(tc,p,"h",pressure,0.,False)
z_plane = wrf_user_intrp3d( z,p,"h",pressure,0.,False)
rh_plane = wrf_user_intrp3d(rh,p,"h",pressure,0.,False)
u_plane = wrf_user_intrp3d( u,p,"h",pressure,0.,False)
v_plane = wrf_user_intrp3d( v,p,"h",pressure,0.,False)
spd = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec
spd at description = "Wind Speed"
spd at units = "m/s"
u_plane = u_plane*1.94386 ; kts
v_plane = v_plane*1.94386 ; kts
u_plane at units = "kts"
v_plane at units = "kts"
res at gsnDraw = False ; don't draw
res at gsnFrame = False ; don't advance frame
; Plotting options for T
opts = res
opts at cnLineColor = "Red"
opts at ContourParameters = (/ 5.0 /)
opts at cnInfoLabelOrthogonalPosF = 0.07 ; offset second label information
opts at gsnContourLineThicknessesScale = 2.0
contour_tc = wrf_contour(a,wks,tc_plane,opts)
delete(opts)
; Plotting options for RH
opts = res
opts at cnFillOn = True
opts at pmLabelBarOrthogonalPosF = -0.1
opts at ContourParameters = (/ 10., 90., 10./)
opts at cnFillColors = (/"White","White","White", \
"White","Chartreuse","Green",\
"Green3","Green4", \
"ForestGreen","PaleGreen4"/)
contour_rh = wrf_contour(a,wks,rh_plane,opts)
delete(opts)
; Plotting options for Wind Speed
opts = res
opts at cnLineColor = "MediumSeaGreen"
opts at ContourParameters = (/ 10. /)
opts at cnInfoLabelOrthogonalPosF = 0.07 ; offset second label information
opts at gsnContourLineThicknessesScale = 3.0
contour_spd = wrf_contour(a,wks,spd,opts)
delete(opts)
; Plotting options for Wind Vectors
opts = res
opts at FieldTitle = "Wind" ; overwrite Field Title
opts at NumVectors = 47 ; wind barb density
vector = wrf_vector(a,wks,u_plane,v_plane,opts)
delete(opts)
; Plotting options for Geopotential Heigh
opts_z = res
opts_z at cnLineColor = "Blue"
opts_z at gsnContourLineThicknessesScale = 3.0
; MAKE PLOTS
if ( pressure .eq. 850 ) then ; plot temp, rh, height, wind barbs
opts_z at ContourParameters = (/ 20.0 /)
contour_height = wrf_contour(a,wks,z_plane,opts_z)
plot = wrf_map_overlays(a,wks,(/contour_rh,contour_tc,contour_height, \
vector/),pltres,mpres)
end if
if ( pressure .eq. 700 ) then ; plot temp, height, wind barbs
opts_z at ContourParameters = (/ 30.0 /)
contour_height = wrf_contour(a,wks, z_plane,opts_z)
plot = wrf_map_overlays(a,wks,(/contour_tc,contour_height, \
vector/),pltres,mpres)
end if
if ( pressure .eq. 500 ) then ; plot temp, height, wind barbs
opts_z at ContourParameters = (/ 60.0 /)
contour_height = wrf_contour(a,wks, z_plane,opts_z)
plot = wrf_map_overlays(a,wks,(/contour_tc,contour_height, \
vector/),pltres,mpres)
end if
if ( pressure .eq. 300 ) then ; plot windspeed, height, wind barbs
opts_z at ContourParameters = (/ 60.0 /)
contour_height = wrf_contour(a,wks, z_plane,opts_z)
plot = wrf_map_overlays(a,wks,(/contour_spd,contour_height, \
vector/),pltres,mpres)
end if
delete(opts_z)
shp= gsn_add_shapefile_polylines(wks,plot,shp_filename,shpres)
draw(plot)
frame(wks)
end do ; END OF LEVEL LOOP
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end do ; END OF TIME LOOP
end
|
|
| |
[ncl-talk] I can't manage to create a plot with multiple lines on it, ea...
|
|
|
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190422/f5d70e6e/attachment.html>
More information about the ncl-talk
mailing list