[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