[ncl-talk] plotting station data over WRF contour map: colored solid dots with black border?
Tabish Ansari
tabishumaransari at gmail.com
Mon Apr 30 12:35:51 MDT 2018
Hi
I want to plot surface station data on top of a WRF contour map. I want the
polymarkers to be circled with a black border - there's no polymarker that
can do this at the moment. So, I'm trying to plot gsmarker=4 type hollow
markers on top of solid coloured polymarkers (gsmarker=16) but unable to
achieve the desired result.
Here's my script:
*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";----------------------------------------------------------------------;
Procedure to attach colored markers to an NCL map, given 1D; data, and
lat/lon arrays, a set of levels, and a color
map.;----------------------------------------------------------------------procedure
add_obs_markers_to_map(wks,plot,levels,colormap,data_1d,lat_1d,lon_1d)local
mkres, nlevels, colors, nlevels, n, iibegin;---Set resources for the
markers mkres = True mkres at gsMarkerIndex = 16;; Based on the levels we
had for the contour plot and the colormap used,; get an array of colors
that span the colormap. This uses the same; algorithm that the contour
used, so the colors will be identical.; nlevels = dimsizes(levels) colors
= span_color_rgba(colormap,nlevels+1) ; Need one more color than number
of levels;; Loop through each level, gather all the data values in the
given level; range, and add colored markers to the existing map.; do
n=0,nlevels;---These "if" statements are how color contours are handled in
NCL. if(n.eq.0) then ii := ind(data_1d.lt.levels(n)) else
if(n.eq.nlevels) then ii := ind(data_1d.ge.levels(n-1)) else
ii := ind(data_1d.ge.levels(n-1).and.data_1d.lt.levels(n)) end if end
if if(any(ismissing(ii))) then continue end if;---Add the
markers mkres at gsMarkerColor = colors(n,:) ; colors is an N x 4
array plot@$unique_string("markers")$ =
gsn_add_polymarker(wks,plot,lon_1d(ii),lat_1d(ii),mkres) end
doend;----------------------------------------------------------------------;
Main
code;----------------------------------------------------------------------begin;---Open
WRF output file and read data DATADir =
"/data2/tabish/control-run-so4-ECMWF/"FILES = systemfunc (" ls -1 " +
DATADir + "wrfout_d03_2014-10* ")numFILES = dimsizes(FILES)s=1a =
addfile(FILES(0)+".nc","r"); THIS IS JUST FOR INITIALIZATION OF PM25 ARRAY,
IT WILL BE SUBSTRACTED FROM THE SUM IN THE ENDpm25init =
a->PM2_5_DRY(15,0,:,:)pm25sum = pm25inittimes =
wrf_user_getvar(a,"times",-1) ; get all times in the filentimes =
dimsizes(times)do it = 16,ntimes-1,1 ; TIME LOOP s=s+1
print("Working on time: " + times(it) ) pm25sum = pm25sum +
a->PM2_5_DRY(it,0,:,:)end do ; END OF TIME LOOPdo ifil =
1,numFILES-2 ; FILE LOOP a = addfile(FILES(ifil)+".nc","r")
; Open the next file times = wrf_user_getvar(a,"times",-1) ; get all
times in the file ntimes = dimsizes(times) do it =
0,ntimes-1,1 ; TIME LOOP s=s+1 print("Working on time: "
+ times(it) ) pm25sum = pm25sum + a->PM2_5_DRY(it,0,:,:) end do ; END
OF TIME LOOPend doa = addfile(FILES(numFILES-1)+".nc","r") ; Open the
next file times = wrf_user_getvar(a,"times",-1) ; get all times in the
file ntimes = dimsizes(times) do it = 0,16,1 ; TIME LOOP
s=s+1 print("Working on time: " + times(it) ) pm25sum = pm25sum +
a->PM2_5_DRY(it,0,:,:) end do ; END OF TIME LOOPpm25_avg=
pm25sum/s;pm25_avg at description = "NO2 conc in ppbv" lat2d =
wrf_user_getvar(a,"XLAT",it) ; latitude/longitude lon2d =
wrf_user_getvar(a,"XLONG",it) pm25_avg at lat2d = lat2d ; for plotting
pm25_avg at lon2d = lon2d;---Will use this for contours and filled markers
colormap = "BlAqGrYeOrReVi200" levels = ispan(0,150,15); levels =
(/0,10,20,30,40,50,60,70,80,90,100,125,150,175,200,225,250/) * 0.01 wks =
gsn_open_wks("x11","wrf_obs_d03_12-31octmean_pm25");---Common resources
shared by both contour plot and marker plot res =
True res at gsnMaximize = False res at gsnLeftString = ""
res at gsnRightString = "" res at mpDataBaseVersion = "MediumRes" ;
res at mpFillOn = False ; res at mpMinLatF =
min(pm25_avg at lat2d) ; res at mpMaxLatF = max(pm25_avg at lat2d) ;
res at mpMinLonF = min(pm25_avg at lon2d) ;
res at mpMaxLonF = max(pm25_avg at lon2d) ;
res at pmTickMarkDisplayMode = "Always" ; better map
tickmarks res at mpLimitMode = "Corners" ; choose range of
map res at mpLeftCornerLatF = 36.05806 res at mpLeftCornerLonF =
113.3356 res at mpRightCornerLatF = 42.4514 res at mpRightCornerLonF =
120.2583 res at mpProjection =
"LambertConformal" res at mpLambertParallel1F = 20 res at mpLambertParallel2F =
50 res at mpLambertMeridianF = 110;---Resources for filled contour plot
cnres = res cnres at cnFillOn = True
cnres at cnLevelSelectionMode = "ExplicitLevels" cnres at cnLevels =
levels cnres at cnFillPalette = colormap cnres at cnLinesOn =
False cnres at cnLineLabelsOn = False; cnres at lbLabelBarOn =
False ; will add in panel later cnres at pmTitleZone = 2 ;
move title down cnres at gsnAddCyclic = False
cnres at mpLimitMode = "Corners" ; choose range of map
cnres at mpLeftCornerLatF = 36.05806 cnres at mpLeftCornerLonF = 113.3356
cnres at mpRightCornerLatF = 42.4514 cnres at mpRightCornerLonF = 120.2583
cnres at mpProjection = "LambertConformal" cnres at mpLambertParallel1F
= 20 cnres at mpLambertParallel2F = 50 cnres at mpLambertMeridianF = 110
cnres at mpFillOn = True
cnres at mpOutlineDrawOrder = "PostDraw"
cnres at mpFillDrawOrder = "PreDraw" ;
cnres at mpOutlineBoundarySets = "National"
cnres at mpOutlineBoundarySets = "Allboundaries"
cnresmpNationalLineColor = "Black" cnresmpGeophysicalLineColor =
"Black" cnres at mpUSStateLineDashPattern = 2 cnres at mpOutlineOn
= True cnres at mpDataBaseVersion = "MediumRes"
cnres at mpDataSetName = "Earth..4" ; U.S. counties
cnres at mpGridAndLimbOn = True cnres at mpGridLineColor = "Black"
cnres at mpFillColors = (/"transparent","transparent","gray","transparent"/)
cnres at pmTickMarkDisplayMode = "Always" cnres at tiMainString = "WRF model
data" cnres at tiDeltaF = 3.3 cnres at tfDoNDCOverlay = True cnres at gsnDraw =
False ; don't draw cnres at gsnFrame =
False ; don't advance frameprint(res);---Create
smooth contour plot plot_wrf =
gsn_csm_contour_map(wks,pm25_avg,cnres);----------------------------------------------------------------------;
Recreate similar plot, but by using filled markers drawn over a
map.;----------------------------------------------------------------------
fname = "surfacedata_12-31_Octmean_improved.txt" lines =
asciiread(fname,-1,"string") lon = tofloat(str_get_field(lines(:),2," "))
lat = tofloat(str_get_field(lines(:),3," ")) pm25 =
tofloat(str_get_field(lines(:),4," "));---Create a map plot for adding
markers to. res at gsnDraw = False res at gsnFrame = False
res at tfPolyDrawOrder = "Draw" ; Necessary to make sure filled markers
drawn *under* map outlines res at tiMainString = "'Observational' data"
res at mpFillOn = True res at mpOutlineDrawOrder =
"PostDraw" res at mpFillDrawOrder = "PreDraw" ;
res at mpOutlineBoundarySets = "National"
res at mpOutlineBoundarySets = "Allboundaries"
res at mpUSStateLineDashPattern = 2 res at mpOutlineOn = True
res at mpDataBaseVersion = "MediumRes" res at mpDataSetName =
"Earth..4" ; U.S. counties ; res at mpGridAndLimbOn = True
res at mpFillColors = (/"transparent","transparent","gray","transparent"/)
res at pmTickMarkDisplayMode = "Always" res at tiMainString = "Surface data"
;plot_obs = gsn_csm_map(wks,res);---Add colored markers to the map
add_obs_markers_to_map(wks,plot_wrf,levels,colormap,pm25,lat,lon);*****************************************;
add black circles over obs data;*****************************************
polyres = True ; poly marker mods desired
polyres at gsMarkerIndex = 4 ; choose circle as polymarker
polyres at gsMarkerSizeF = 10.0 ; select size to avoid streaking
polyres at gsMarkerColor = (/"black"/) ; choose color
gsn_polymarker(wks,plot_wrf,lat,lon,polyres) ; draw polymarkers;---This
draws the map with the filled markers draw(plot_wrf) frame(wks)end*
Please let me know if that's possible.
Thanks
Tabish
Tabish U Ansari
PhD student, Lancaster Environment Center
Lancaster Univeristy
Bailrigg, Lancaster,
LA1 4YW, United Kingdom
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180430/5bd35771/attachment.html>
More information about the ncl-talk
mailing list