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 my array of longitude-latitude points?

UPDATE: I tested some more on how to plot the polygon. I removed the sort routine and just used the data in the order they occur in the file. This seems to improve the result, but as @tcaswell mentioned, the polygon shape still undercuts itself (see new plot). I am hoping that there could be a path/polygon routine that could solve my problem and sort of merge all shapes or paths within the boundaries of the polygon. Suggestions are very welcome.

UPDATE 2:

I have now a working version of my script that is based on suggestions by @Rutger Kassies and Roland Smith. I ended up reading the Shapefile using org which worked relatively well. It worked well for the standard lmes_64.shp file but when I used more detailed LME files where each LME could consist of several polygons this script broke down. I would have to find a way to merge the various polygons for identical LME names to make that work. I attach my script I ended up with in case anyone would take a look at it. I very much appreciate comments for how to improve this script or to make it more generic. This script creates the polygons and extracts data within the polygon region that I read from a netcdf file. The grid of the input file is -180 to 180 and -90 to 90.

import numpy as np
import math
from pylab import *
import matplotlib.patches as patches
import string, os, sys
import datetime, types
from netCDF4 import Dataset
import matplotlib.nxutils as nx
from mpl_toolkits.basemap import Basemap
import ogr
import matplotlib.path as mpath
import matplotlib.patches as patches


def getLMEpolygon(coordinatefile,mymap,index,first):

    ds = ogr.Open(coordinatefile)
    lyr = ds.GetLayer(0)
    numberOfPolygons=lyr.GetFeatureCount()

    if first is False:
        ft = lyr.GetFeature(index)
        print "Found polygon:",  ft.items()['LME_NAME']
        geom = ft.GetGeometryRef()

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

        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
            all_XY +=mymap(x,y)


        if len(all_XY)==0:
            all_XY=None
            mypoly=None
        else:
            mypoly=np.empty((len(all_XY[:][0]),2))
            mypoly[:,0]=all_XY[:][0]
            mypoly[:,1]=all_XY[:][3]
    else:
        print "Will extract data for %s polygons"%(numberOfPolygons)
        mypoly=None
    first=False
    return mypoly, first, numberOfPolygons


def openCMIP5file(CMIP5name,myvar,mymap):
    if os.path.exists(CMIP5name):
        myfile=Dataset(CMIP5name)
        print "Opened CMIP5 file: %s"%(CMIP5name)
    else:
        print "Could not find CMIP5 input file %s : abort"%(CMIP5name)
        sys.exit()
    mydata=np.squeeze(myfile.variables[myvar][-1,:,:]) - 273.15
    lonCMIP5=np.squeeze(myfile.variables["lon"][:])
    latCMIP5=np.squeeze(myfile.variables["lat"][:])

    lons,lats=np.meshgrid(lonCMIP5,latCMIP5)

    lons=lons.flatten()
    lats=lats.flatten()
    mygrid=np.empty((len(lats),2))
    mymapgrid=np.empty((len(lats),2))

    for i in xrange(len(lats)):
        mygrid[i,0]=lons[i]
        mygrid[i,1]=lats[i]
        X,Y=mymap(lons[i],lats[i])
        mymapgrid[i,0]=X
        mymapgrid[i,1]=Y

    return mydata, mygrid, mymapgrid

def drawMap(NUM_COLORS):

    ax = plt.subplot(111)
    cm = plt.get_cmap('RdBu')
    ax.set_color_cycle([cm(1.*j/NUM_COLORS) for j in range(NUM_COLORS)])
    mymap = Basemap(resolution='l',projection='robin',lon_0=0)

    mymap.drawcountries()

    mymap.drawcoastlines()
    mymap.fillcontinents(color='grey',lake_color='white')
    mymap.drawparallels(np.arange(-90.,120.,30.))
    mymap.drawmeridians(np.arange(0.,360.,60.))
    mymap.drawmapboundary(fill_color='white')

    return ax, mymap, cm

"""Edit the correct names below:"""

LMEcoordinatefile='ShapefileBoundaries/lmes_64.shp'
CMIP5file='tos_Omon_CCSM4_rcp85_r1i1p1_200601-210012_regrid.nc'

mydebug=False
doPoints=False
first=True


"""initialize the map:"""
mymap=None
mypolyXY, first, numberOfPolygons = getLMEpolygon(LMEcoordinatefile, mymap, 0,first)
NUM_COLORS=numberOfPolygons
ax, mymap, cm = drawMap(NUM_COLORS)


"""Get the CMIP5 data together with the grid"""
SST,mygrid, mymapgrid = openCMIP5file(CMIP5file,"tos",mymap)

"""For each LME of interest create a polygon of coordinates defining the boundaries"""
for counter in xrange(numberOfPolygons-1):

    mypolyXY,first,numberOfPolygons = getLMEpolygon(LMEcoordinatefile, mymap,counter,first)

    if mypolyXY != None:
        """Find the indices inside the grid that are within the polygon"""
        insideBoolean = plt.mlab.inside_poly(np.c_[mymapgrid[:,0],mymapgrid[:,1]],np.c_[mypolyXY[:,0],mypolyXY[:,1]])
        SST=SST.flatten()
        SST=np.ma.masked_where(SST>50,SST)

        mymapgrid=np.c_[mymapgrid[:,0],mymapgrid[:,1]]
        myaverageSST=np.mean(SST[insideBoolean])

        mycolor=cm(myaverageSST/SST.max())
        scaled_z = (myaverageSST - SST.min()) / SST.ptp()
        colors = plt.cm.coolwarm(scaled_z)

        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)

        if doPoints is True:

            for point in xrange(len(insideBoolean)):
                pointX=mymapgrid[insideBoolean[point],0]
                pointY=mymapgrid[insideBoolean[point],1]
                ax.scatter(pointX,pointY,8,color=colors)
                ax.hold(True)


if doPoints is True:
    colorbar()
print "Extracted average values for %s LMEs"%(numberOfPolygons)
plt.savefig('LMEs.png',dpi=300)
plt.show()

Final image attached. Thanks for all help.

Cheers, Trond

Answer 1

Having an array of points is not enough. You need to know the order of the points. Normally the points of a polygon are given sequentially. So you draw a line from the first point to the second, then from the second to the third et cetera.

If your list is not in sequential order, you need extra information to be able to make a sequential list.

A shapefile (see the documentation) contains a list of shapes, like a Null shape, Point, PolyLine, Polygon, with variants containing also the Z and M (measure) coordinates. So just dumping the points will not do. You have to divide them up in the different shapes and render the ones you are interested in. In this case probably a PolyLine or Polygon. See the link above for the data format for these shapes. Keep in mind that some parts of the file are big-endian, while others are little-endian. What a mess.

I would suggest using the struct module to parse the binary .shp file, because again according to the documentation, the points of a single Polygon are in order, and they form a closed chain (the last point is the same as the first).

Another thing you could try with your current list of coordinates is to start with a point, and then look for the same point furtheron in the list. Everything between those identical points should be one polygon. This is probably not foolproof, but see how far it gets you.

Answer 2

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.

Answer 3

There is a blog post here which talks about shapefiles and basemap: http://www.geophysique.be/2013/02/12/matplotlib-basemap-tutorial-10-shapefiles-unleached-continued/

If you're keen to try it out, cartopy might also be an option. Plotting data from a shapefile is designed to be fairly easy:

import matplotlib.pyplot as plt

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader

# pick a shapefile - Cartopy makes it easy to access Natural Earth
# shapefiles, but it could be anything
shapename = 'admin_1_states_provinces_lakes_shp'
states_shp = shpreader.natural_earth(resolution='110m',
                                      category='cultural', name=shapename)

.

# states_shp is just a filename to a shapefile
>>> print states_shp
/Users/pelson/.local/share/cartopy/shapefiles/natural_earth/cultural/110m_admin_1_states_provinces_lakes_shp.shp

.

# Create the mpl axes of a PlateCarree map
ax = plt.axes(projection=ccrs.PlateCarree())

# Read the shapes from the shapefile into a list of shapely geometries.
geoms = shpreader.Reader(states_shp).geometries()

# Add the shapes in the shapefile to the axes
ax.add_geometries(geoms, ccrs.PlateCarree(),
                  facecolor='coral', edgecolor='black')

plt.show()