For a project I am working on, I needed to figure out which elecion districts make up a New York City Council District. While this would be a simple lookup to conduct for one district, it would be very time consuming to do for all fifty one districts in New York City. Fortunately, the wonderful QGIS comes with a built-in Python console that makes solving these sorts of problems much more managable.
The method that I used involved working directly with the QGIS console. There is a way to access the associated QGIS library without using QGIS explicitly, but I wanted to be able to check on things as they progressed visually, so my approach involved using QGIS directly. Note also that I’m using QGIS version 2.2 Valmiera (which has a much improved Python console from 1.8 Lisboa, along with a revamped API).
After loading up your shapefiles into QGIS, accessing all of the layers is fairly easy:
From this point, we need to break each QGIS Vector layer down into discrete features. The easiest way to do that is this:
This gives an iterator object. At this point, getting the geographic intersection is pretty simple.
That approach, though, can be really slow if you have lots of feature elements to iterate through. This can be sped up significantly by using a spatial index. For details about the timing of the speed increase, check out this blog post by Nathan Woodrow.
We start by creating an empty
Now, we search for the intersection based on the index. This will return the same results, but will be much faster.
Now that we have a method to select out the IDs that are contained within any given parent feature, we have to add those IDs to the parent layer in order to make use of them. We are going to add our IDs as one long comma-separated string and parse that later. Adding a new field (or attribute) to a layer is fairly simple:
Now that we’ve added the new field, we have to add the data to it. This is also fairly straightforward.
Now, after we run this script from the console, we should see the list of districts populated in the proper field.
The full code that I ending up using is available at this gist.