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###'
print
count = 0
output_file = open('GridCount.csv', 'a')
Multi = MultiPolygon([shape(Y['geometry']) for Y in fiona.open('yourpolygons.shp')])
print Multi
for X in fiona.collection('yourgrid.shp'):
intsct = 0
count += 1
print 'ID', count
output_file.write(str(count) + &amp;amp;quot;&amp;amp;amp;&amp;amp;quot;)
grid_shape = shapely.geometry.shape(X['geometry'])
print grid_shape
for Y in Multi:
if Y.intersects(grid_shape):
intsct += 1
else:
pass
print 'INTERSECTS', intsct
output_file.write(str(intsct) + &amp;amp;quot;\n&amp;amp;quot;)
output_file.close()
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.

### Like this:

Like Loading...

*Related*