4

I'm using Python, Shapely and Fiona to process some shapefiles.

Is there an easy way to count the number of points in a polygon/multipolygon in Shapely and Fiona?

8

1) With Fiona, you don't need shapely to count the number of points in a polygon/multipolygon. Simply use the resulting GeoJSON format (= a Python dictionary).

Polygon simple:

>>> import fiona
>>> shape = fiona.open("simplePoly.shp")
>>> # first feature
>>> feature = shape.next()
>>> geom = feature['geometry']
>>> print geom
{'type': 'Polygon', 'coordinates': [[(1.0, 1.0), (1.0, 7.0), (7.0, 7.0), (7.0, 1.0), (1.0, 1.0)]]}
>>> print geom['coordinates']
[[(1.0, 1.0), (1.0, 7.0), (7.0, 7.0), (7.0, 1.0), (1.0, 1.0)]]
>>> print len(geom['coordinates'])
1 # -> only one polygon
>>> print "number of points =", len(geom['coordinates'][0])
number of points = 5 

or simply:

>>> shape = fiona.open("simplePoly.shp")
>>> print "number of points =", len(shape.next()['geometry']['coordinates'][0]
number of points = 5

Polygon with one hole

>>> shape = fiona.open("Poly_hole.shp")
>>> geom = shape.next()['geometry']['coordinates']
[[(1.0, 1.0), (1.0, 7.0), (7.0, 7.0), (7.0, 1.0), (1.0, 1.0)], [(2.0, 3.0), (4.0, 3.0), (4.0, 5.0), (2.0, 5.0), (2.0, 3.0)]]
>>> print len(geom)
2 # -> polygon and 1 hole
>>> print "number of exterior points =", len(geom)[0])
number of exterior points =5
>>> print "number of interior points =", len(geom)[1])
number of interior points =5 

MultiPolygon

>>> shape = fiona.open("multiPoly.shp")
>>> geom = shape.next()['geometry']['coordinates']
>>> print geom
{'type': 'MultiPolygon', 'coordinates': [[[(1.0, 1.0), (1.0, 7.0), (7.0, 7.0), (7.0, 1.0), (1.0, 1.0)]], [[(2.0, 3.0), (4.0, 3.0), (4.0, 5.0), (2.0, 5.0), (2.0, 3.0)]]]}
>>> print "number of points =", sum([len(poly[0]) for poly in geom])
number of points =10

So the solution is

if shape.next()['geometry']['type'] == 'Polygon':
      if len(shape.next()['geometry']) == 1: # -> polygon simple
      if len(shape.next()['geometry']) > 1: # -> polygon with holes
if shape.next()['geometry']['type'] == 'MultiPolygon': # -> Multipolygon

2) if you want to use the solution of againstflow with Fiona and shapely, use the shapely shape function (Python Geo Interface):

>>> from shapely.geometry import shape
>>> shape = fiona.open("simplePoly.shp")
>>> geom_shapely = shape(shape.next())
>>> print type(geom_shapely)
<class 'shapely.geometry.polygon.Polygon'>
>>> print geom_shapely
POLYGON ((1 1, 1 7, 7 7, 7 1, 1 1))
6

You can get polygon exterior and interior points coordinates this way:

def extract_poly_coords(geom):
    if geom.type == 'Polygon':
        exterior_coords = geom.exterior.coords[:]
        interior_coords = []
        for interior in geom.interiors:
            interior_coords += interior.coords[:]
    elif geom.type == 'MultiPolygon':
        exterior_coords = []
        interior_coords = []
        for part in geom:
            epc = extract_poly_coords(part)  # Recursive call
            exterior_coords += epc['exterior_coords']
            interior_coords += epc['interior_coords']
    else:
        raise ValueError('Unhandled geometry type: ' + repr(geom.type))
    return {'exterior_coords': exterior_coords,
            'interior_coords': interior_coords}

Then you can count for example all exterior coords:

polyPoints = len(exterior_coords);

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.