Create closed polygon from boundary points

Question

I have an array of longitude-latitude points that defines the boundaries of an area. I would like to create a polygon based on these points and plot the polygon on a map and fill it. Currently, my polygon seems to consist of many patches that connect all the points, but the order of the points is not correct and when I try to fill the polygon I get a weird looking area (see attached). I sort my longitude-latitude points (mypolyXY array) according to the center of the polygon, but my guess is that this is not correct: ``` cent=(np.sum([p[0] for p in mypolyXY])/len(mypolyXY),np.sum([p[1] for p in mypolyXY])/len(mypolyXY)) # sort by polar angle mypolyXY.sort(key=lambda p: math.atan2(p[1]-cent[1],p[0]-cent[0])) ``` I plot the point locations (black circles) and my polygons (colored patches) using ``` scatter([p[0] for p in mypolyXY],[p[1] for p in mypolyXY],2) p = Polygon(mypolyXY,facecolor=colors,edgecolor='none') ax.add_artist(p) ``` My question is: how can I close my polygon based on m…

Accepted Answer

I recommend using the original Shapefile, which is in a format appropriate for storing polygons. As an alternative to OGR you could use Shapely, or export the polygon to Wkt etc.

```
import ogr
import matplotlib.path as mpath
import matplotlib.patches as patches
import matplotlib.pyplot as plt

ds = ogr.Open('lmes_64.shp')
lyr = ds.GetLayer(0)
ft = lyr.GetFeature(38)
geom = ft.GetGeometryRef()
ds = None

codes = []
all_x = []
all_y = []

if (geom.GetGeometryType() == ogr.wkbPolygon):
  for i in range(geom.GetGeometryCount()):

r = geom.GetGeometryRef(i)
    x = [r.GetX(j) for j in range(r.GetPointCount())]
    y = [r.GetY(j) for j in range(r.GetPointCount())]

codes += [mpath.Path.MOVETO] + (len(x)-1)*[mpath.Path.LINETO]
    all_x += x
    all_y += y

if (geom.GetGeometryType() == ogr.wkbMultiPolygon):
  codes = []
  for i in range(geom.GetGeometryCount()):
    # Read ring geometry and create path
    r = geom.GetGeometryRef(i)
    for part in r:
      x = [part.GetX(j) for j in range(part.GetPointCount())]
      y = [part.GetY(j) for j in range(part.GetPointCount())]
      # skip boundary between individual rings
      codes += [mpath.Path.MOVETO] + (len(x)-1)*[mpath.Path.LINETO]
      all_x += x
      all_y += y

carib_path = mpath.Path(np.column_stack((all_x,all_y)), codes)    
carib_patch = patches.PathPatch(carib_path, facecolor='orange', lw=2)

poly1 = patches.Polygon([[-80,20],[-75,20],[-75,15],[-80,15],[-80,20]], zorder=5, fc='none', lw=3)
poly2 = patches.Polygon([[-65,25],[-60,25],[-60,20],[-65,20],[-65,25]], zorder=5, fc='none', lw=3)

fig, ax = plt.subplots(1,1)

for poly in [poly1, poly2]:
    if carib_path.intersects_path(poly.get_path()):
        poly.set_edgecolor('g')
    else:
        poly.set_edgecolor('r')

ax.add_patch(poly)

ax.add_patch(carib_patch)
ax.autoscale_view()
```

Also checkout Fiona (wrapper for OGR) if you want really easy Shapefile handling.

Create closed polygon from boundary points

Solution

Problem