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.


Where is Western Sydney?

‘Western Sydney’ is a pretty politically and culturally charged term in Sydney, see:

But, there is no real definition of what Western Sydney is. Seeing as it is essentially an opinon, I figured asking people is as good a way as any to define it.

After seeing this sweet post by  and , I realised I could use the  use the tools generously uploaded by Nick Martinelli here to allow people to submit their definition of Western Sydney and then combine them.

The response was pretty good, and once I had removed all the continent sized cock and balls drawings, submissions for the Sydney in Nova Scotia etc.  I had 119 definitions of Western Sydney to compare.

What was surprisingly difficult was finding a way to tally up how many definitions covered a given location. QGIS does not appear to have an integrated tool to achieve this. What I ended up doing was building an arbitrary hexagonal grid and then counting the number of the definitions that fell within each of the cells in this grid. To do that I had to build a python script using fiona and shapely, I’ll chuck the details of that code up later, it could help out anyone else who needs to count the number of polygons in another set of polygons.

The results are shown below. Each category represents the proportion of submitted definitions that agree a cell is in Western Sydney.

People seem to agree the eastern extent is somewhere between Lidcombe and Strathfield and the Western extent is somewhere between Penrith and the escarpment of the Blue Mountains. Compared to the northern and southern extent there is SIGNIFICANTLY less variation in opinion.

In general the spread of definitions to the north seems narrower than the south. The northern extent of 75% consensus area seems to follow the western train line almost exactly.

I also found it intersting that roughly equal numbers of  people consider Campbelltown and parts of Epping to both be Western Sydney.

Edit: I have updated this post, removing references to the website people entered their definitions on that I was no longer using.