[ncl-talk] Draw only outer boundaries of a shapefiles polygon

[Tabish Ansari]

Hi

I'm trying to map three distinct regions, two of which are a collection of
several districts, and create boundaries for them on a map. I have
shapefiles for each district but drawing polylines using all the member
districts would also draw internal boundaries which I don't want. I only
want to draw outer boundaries for each region. Please see the attached
plot.

My script is pasted here:




























































































































*begin  wks  = gsn_open_wks("png","pulse-districts") ; send graphics to PNG
file  res                     = True  res at gsnDraw             = False
; don't draw yet  res at gsnFrame            = False       ; don't advance
frame yet  res at gsnMaximize         = True       ; maximize plot in frame;
 res at gsnPaperOrientation = "portrait"  res at vpKeepAspect        = True
res at mpProjection        = "LambertConformal"    ; choose projection
res at mpLambertParallel1F = 20          ; first parallel
res at mpLambertParallel2F = 50          ; second parallel
res at mpLambertMeridianF  = 110         ; meridian  res at mpLimitMode         =
"Corners"   ; corner method of zoom  res at mpLeftCornerLatF  = 35
res at mpLeftCornerLonF  = 109.35  res at mpRightCornerLatF = 44.41
res at mpRightCornerLonF = 124.94  res at pmTickMarkDisplayMode = "Always"  ;
turn on tickmarks  res at tmXBMajorOutwardLengthF = 0.0
res at tmXTMajorOutwardLengthF = 0.0  res at tmYLMajorOutwardLengthF = 0.0
res at tmYRMajorOutwardLengthF = 0.0   res at tiMainFontHeightF  = 0.025
res at tiYAxisFontHeightF = 0.01  res at tiXAxisFontHeightF = 0.01
res at tmXBLabelFontHeightF = 0.01  plot = gsn_csm_map(wks,res)   ; Create
map, but don't draw it
yet.;*************************************************; Section to add
polygons to map.;*************************************************  f =
addfile("/home/tabish/Shapefiles/China/CHN_adm2.shp", "r")   ; Open
shapefile  ;; Read data off shapefile;  segments = f->segments  geometry =
f->geometry  segsDims = dimsizes(segments)  geomDims = dimsizes(geometry);;
Read global attributes  ;  geom_segIndex = f at geom_segIndex  geom_numSegs  =
f at geom_numSegs  segs_xyzIndex = f at segs_xyzIndex  segs_numPnts  =
f at segs_numPnts  lines       = new(segsDims(0),graphic)   ; Array to hold
polygons  numFeatures = geomDims(0)  name =  f->NAME_2  lon   = f->x  lat
= f->y  plres             = True       ; resources for polylines
plres at gsEdgesOn   = True       ; draw border around polygons
plres at gsEdgeColor = "grey"    ;  plres at gsFillColor      = "Navy";
 plres at gsFillOpacityF    = 0.5    segNum = 0  do i=0, numFeatures-1       ;
color assignment (probably a better way to do this?)     if
(name(i).eq."Beijing") then        plres at gsFillColor = "tomato"
 else if
(name(i).eq."Langfang").or.(name(i).eq."Baoding").or.(name(i).eq."Shijiazhuang").or.(name(i).eq."Xingtai").or.(name(i).eq."Handan").or.(name(i).eq."Tangshan").or.(name(i).eq."Tianjin").or.(name(i).eq."Cangzhou").or.(name(i).eq."Hengshui").or.(name(i).eq."Dezhou").or.(name(i).eq."Binzhou").or.(name(i).eq."Dongying").or.(name(i).eq."Liaocheng").or.(name(i).eq."Jinan").or.(name(i).eq."Zibo")then
      plres at gsFillColor = "limegreen"          else
plres at gsFillColor = "dodgerblue"     end if     end if     startSegment =
geometry(i, geom_segIndex)     numSegments  = geometry(i, geom_numSegs)
 do seg=startSegment, startSegment+numSegments-1        startPT =
segments(seg, segs_xyzIndex)        endPT = startPT + segments(seg,
segs_numPnts) - 1        lines(segNum) = gsn_add_polygon(wks, plot,
lon(startPT:endPT),  \
 lat(startPT:endPT), plres)        segNum = segNum + 1     end do
plres at gsFillColor      = "Grey"  end doprint("done making
polygons")delete(lat)delete(lon);---Attach some dummy text strings  txres
            = True  txres at txFontHeightF = 0.020  txres at txJust        =
"CenterCenter";STARTING TO MARK CITIES WHERE APEC CONTROLS WERE
IMPLEMENTEDlat=40.15lon=116.4074text = "Beijing"text_id1 =
gsn_add_text(wks,plot,text,lon,lat,txres)lat=38lon=116.2text = "Near
Neighbourhood"text_id8 =
gsn_add_text(wks,plot,text,lon,lat,txres)lat=41lon=112.2text = "Far
Neighbourhood"text_id8 = gsn_add_text(wks,plot,text,lon,lat,txres);
Maximize output in frame. This will draw everything and advance; the
frame.;  maximize_output(wks,False)end*

Thanks very much.

Cheers,

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/20190530/912a8efb/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pulse-districts.png
Type: image/png
Size: 195698 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190530/912a8efb/attachment.png>