Counting sides or edges of polygon using open source GIS?

Question

I wondered if it is possible to determine the number of sides / edges (not vertices) that a polygon has. Imagine the following case. There is a shapefile with 10 polygons.

Is there any way to determine for each polygon, how many sides there are and then write these into the attribute table?

I would be happy for any solution using open source tools.

In addition, I aim to use that in for describing voronoi polyons.

Answer 1

With Python and Fiona, Polygons and MutiPolygons (multi-parts) are different geometries:

1) multi-parts geometries

 import fiona
 shape = fiona.open("polygons.shp")
 # shapefile schema
 print c.schema
 {'geometry': 'Polygon', 'properties': OrderedDict([(u'id', 'int:10')])}
 # first feature
 first = shape.next()
 print first (GeoJSON format)
 {'geometry': {'type': 'Polygon', 'coordinates': [[(244697.45179524383, 1000369.2307574936), (244827.15493968062, 1000373.0455558595), (244933.96929392271, 1000353.9715640305), (244933.96929392271, 1000353.9715640305), (244930.15449555693, 1000147.9724522779), (244697.45179524383, 1000159.4168473752), (244697.45179524383, 1000369.2307574936)]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 1)])}


 multi = fiona.open("multipolygons.shp")
 print c.schema
 {'geometry': 'Polygon', 'properties': OrderedDict([(u'id', 'int:10')])}
 # first feature
 first = shape.next()
 print first
 {'geometry': {'type': 'MultiPolygon', 'coordinates': [[[(244697.45179524383, 1000369.2307574936), (244827.15493968062, 1000373.0455558595), (244933.96929392271, 1000353.9715640305), (244933.96929392271, 1000353.9715640305), (244930.15449555693, 1000147.9724522779), (244697.45179524383, 1000159.4168473752), (244697.45179524383, 1000369.2307574936)]], [[(246082.22360202507, 1000453.1563215409), (246139.44557751188, 1000460.7859182726), (246189.03795626713, 1000403.5639427857), (246189.03795626713, 1000403.5639427857), (246086.03840039085, 1000132.7132588148), (245990.66844124615, 1000205.1944277647), (246082.22360202507, 1000453.1563215409)]]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 1)])}

2) With simple Polygons, the coordinates are

coord = first['geometry']['coordinates']
[[(244697.45179524383, 1000369.2307574936), (244827.15493968062, 1000373.0455558595), (244933.96929392271, 1000353.9715640305), (244933.96929392271, 1000353.9715640305), (244930.15449555693, 1000147.9724522779), (244697.45179524383, 1000159.4168473752), (244697.45179524383, 1000369.2307574936)]
print len(coord)
1 # -> one polygon

And the equivalent LinearRing coordinates:

linearR = coor[0] # or coord = first['geometry']['coordinates'][0]
print linearR
[(244697.45179524383, 1000369.2307574936), (244827.15493968062, 1000373.0455558595), (244933.96929392271, 1000353.9715640305), (244933.96929392271, 1000353.9715640305), (244930.15449555693, 1000147.9724522779), (244697.45179524383, 1000159.4168473752), (244697.45179524383, 1000369.2307574936)]

So

nb_vertices = len(lineaR) #(x,y) points
print nb_vertices
7 
nb_edges = nb_vertices -1
print nb_edges
6

3) With MultiPolygons, use a for loop

4) If the Polygons have holes, use Shapely with Fiona

holes = fiona.open("poly_holes.shp")
from shapely.geometry import shape
# First feature
first = holes.next()
print first['geometry']['coordinates']
# conversion to shapely geometry
shape = first['geometry']
[[(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)], [(5.0, 5.0), (6.0, 5.0), (6.0, 6.0), (5.0, 6.0), (5.0, 5.0)]]
# exterior coordinates
print list(shape.exterior.coords) # = LinearRing
[(1.0, 1.0), (1.0, 7.0), (7.0, 7.0), (7.0, 1.0), (1.0, 1.0)]
edges_ext= len(list(shape.exterior.coords))
print edges_ext
5

Same with interior coordinates (look at the Shapely manual)

Answer 2

In QGIS if your polygons do not have any holes or multi-parts:

l = iface.activeLayer()
for f in l.getFeatures():
    print f['NAME']
    print 'no. edges: %d' %(len(f.geometry().asPolygon()[0])-1)

replace 'NAME' with some identifier in your layer attribute table.

Concerning writing to the attribute table check the instructions in the PyQGIS Cookbook - Modifying Vector Layers.

Answer 3

I finally found also an R-solution for my problem. The rgeos package offers the get.pts function which allows to extract the number of vertices. As we have learned from the first comments, the number of edges is equal to the number of vertices - 1. Programming a solution for my task in R seems now much easier.