Efficiently identifying overlapping polygons with Python

Question

I'm parsing the NOAA severe weather alert feed and using it as a datasource for a weather alerting system.

In the feed, I occasionally find a polygon entry for an alert. Unfortunately, counties that fall under these polygons are not always identified in the FIPS/UGC codes section of the alert.

In cases where I come across a polygon entry, I need to interate over all my county polygons to make sure all counties potentially affected by weather are included in my output.

Here's an example of a polygon entry from today:

<cap:polygon>41.38,-97.93 41.40,-97.81 41.45,-97.69 41.40,-97.69 41.36,-97.92 41.38,-97.93</cap:polygon>

With 3200+ US County polygons that I have stored, I need to find the most efficient way of iterating over them and comparing the county polygon with the alert polygon to see if they overlap. (I'm not sure if I'm looking for overlap, intersection, or union?)

Does anyone have suggestions on the best way, using Python, to compare a polygon against a list of polygons to identify areas that overlap?

Answer 1

First, welcome to the site!

Without knowing what other GIS software you've got available I'd suggest looking at a spatially enabled database, such as PostGIS to store your data, and accessing the data with GeoAlchemy in Python (I like SpatailLite as a great lightweight alternative, but access in Python is a bit painful as you have to recompile the sqlite3 library from scratch).

The database will automatically build an RTree index which stores extents of objects. It then queries the extents of objects to see if they intersect before querying the actual geometries to significantly increase the speed of the query.

If you don't want to use a database you can build stand-alone RTree indexes in Python with the Rtree module. Construct an RTree from the bounds of your US County polygons and then compare that to the bounds of your incoming polygons to narrow down your query region. Then do a spatial intersection of the counties with your actual geometry to get whether or not they intersect. If you're having trouble constructing geometries I'd also strongly recommend looking at the Shapely library.