Counting Polygons in Polygons

For a recent project I had a shapefile that consisted of a number of polygons, I wanted to be able to see how many polygons over lapped at a given location.

To do this I decided to build another shapefile consisting of an arbitrary grid and then see how many of the polygons intersect the individual grid cells.

So all I needed next was a way to count the number of individual polygons from one shapefile that intersect the individual polygons of another.

There seems to be a fairly robust way to do this sort of thing in ArcGIS, but it does not appear that there is a way to do it in QGIS natively, so I decided to build a python script using fiona and shapely.

The solution I came up with is below:

import shapely
import pprint
import csv

from shapely.geometry import Point
from shapely.geometry import MultiPolygon
from shapely.geometry import shape

import fiona

print '###START###'
count = 0
output_file = open('GridCount.csv', 'a')

Multi = MultiPolygon([shape(Y['geometry']) for Y in'yourpolygons.shp')])
print Multi
for X in fiona.collection('yourgrid.shp'):
    intsct = 0
    count += 1
    print 'ID', count
    output_file.write(str(count) + "&")

    grid_shape = shapely.geometry.shape(X['geometry'])
    print grid_shape
    for Y in Multi:
        if Y.intersects(grid_shape):
            intsct += 1
    print 'INTERSECTS', intsct
    output_file.write(str(intsct) + "\n")    


print '###END###'

This output a csv file with the first column being an identifier for the polygon in the grid shapefile and the second column the number of intersections.

Then I just copied and pasted the data from this new csv into the shapefiles dbf.

It probably wouldn’t be too hard to improve the code so it automatically edits the shapefile itself, but I was too excited to see the results of my project so I didn’t get around to doing it.

Couple last things I’ll point out, I had problems if I imported fiona before shapely and, this code uses intersection as the topology test, so any polygon that intersects a grid cell is included in the results. Shapely supports a number of other topology tests so I imagine it would be pretty simple to change that test for other topological relationships.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s