Finding polygons crossed by line using OGR?

Question

I'm trying to find all the polygons crossed by a single line (a GPS track). I'm using the OGR library (from python) for computing this, but it's currently a bit 'brute-force' (and slow). For every point of my track, I call the intersect method with all the polygons. The obvious optimization is to check only with adjacent polygons. But I guess this is a classical problem, with an already known solution (which I can't find...).

I would like to avoid using a dedicated database as I'm trying to write a standalone software (spatialite is an option if the DB is the way to go).

(FYI, the current source code is available here: https://github.com/dkm/airspace-checker )

Answer 1

For a Python solution, you may want to look at Shapely http://gispython.org/shapely/docs/1.2/ and RTree http://pypi.python.org/pypi/Rtree/

Rtree will help you create spatial indexes.

Answer 2

Instead of expansive intersect, you can perform pre-selection of polygons based on comparison of bounding boxes. In other words, find all polygons overlapped / adjacent to MBR of segments of your track. Then perform detailed test on the subset of polygons.

Answer 3

The proposals of mloskot and Nicklas to compare the bounding boxes are indeed correct.

If you are using shapefiles you could also consider calling this saga module: http://www.saga-gis.org/saga_modules_doc/shapes_transect/index.html

Answer 4

What a database like PostGIS does to speed this up is to first do an index, bounding box compare. It first finds all polygons that have bounding boxes interersecting with the bounding box of the line. The problem in your case might be that the linestring is long and will have a very big bounding box intersecting many polygons that is of no interest.

If the lines are very long you will probably also have to work with geodetic functions that is much more complex and slow than planar functions.

It might be quite complex to make things run smooth.

Why do you not want to rely on a database? That will not solve all your problems, but there is a lot of built in optimisations in PostGIS for instance. There you also have the geodetic calculations of intersection if you need it.

Update: I read your question again and realised that you are not using the linestring the trac forms but each vertex.

I think you are on the wrong trac ;)
Both because you are missing to check if the edge between the vertexpoints intersects the polygon and because you are moving the iteration between the vertex points to python instead of some C implementation which I think is much faster. Then you have that problem with indexes. To make things faster you will have to build and handle some sort of spatial index.

At the other hand, if you are doing this much of the work in your own code, why don't you do the intersection test too. That test is just a point in polygon test if you are dealing with the vertex points. Google for "point in polygon" and you will find several algorithms.

But, I would go for a database driven approach. That will give you the possibilities to use spatial indexes.

/Nicklas