Count the number of points in a (multi)polygon in shapely?

Question

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?

Answer 1

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))

Answer 2

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);

Answer 3

Flatten geometry and count coordinates

I'm a fan of this flatten function, so here's the shapely geometry equivalent, and a coords length calculator.

def geometry_flatten(geom):
  if hasattr(geom, 'geoms'):  # Multi<Type> / GeometryCollection
    for g in geom.geoms:
      yield from geometry_flatten(g)
  elif hasattr(geom, 'interiors'):  # Polygon
    yield geom.exterior
    yield from geom.interiors
  else:  # Point / LineString
    yield geom
def geometry_length(geom):
  return sum(len(g.coords) for g in geometry_flatten(geom))

Example Usage

df = gpd.read_file(gpd.datasets.get_path('nybb'))
g = df.loc[0, 'geometry']
geometry_length(g)  # 8991
df['geometry'].apply(geometry_length)  # 8991, 29219, 22986, 6362, 8505
geometry_length(df.unary_union)  # 75545

similar geopandas question
geopandas question on coordinates