# Best way to find the polygons crossed by a line

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 )

-

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/

-
I've already seen shapely, but didn't know about rtree! Thanks! – Marc Dec 22 '10 at 22:04
Ok, I've been able to (finally) test rtree. It works fine as I can reduce the test set from 1096 to 19 :). I'm thinking of using shapely, but I have to best understand what it implies (shapely does not handle projections, I don't want to loose anything in the process). – Marc Jan 11 '11 at 8:51

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.

-

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

-
I think I'll check solutions with BBoxes. My data is not random: tracks are paraglider/hangglider flights and the polygons are the airspaces. I'll have a look to see if I can come up with a simple filter on the bbox. I'm not directly using shapefiles but I've written a script that read OpenAIR format ("standard" airspace descriptors) to OGR (from which I can directly exports to shp) – Marc Dec 22 '10 at 14:24

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

-
I was thinking the same thing on the linestring. I'm not as familiar with GIS & python, but there should be a way to build spatial indexes in memory (I know various there are various .net options for doing this). This might be another good question for gis.se. – Jay Cummins Dec 22 '10 at 13:15
Thanks for the answer. I'm not using a DB for the setup simplicity. But if avoiding DB implies too much code complexity, then I will change my mind :). As for the "spatial indexes", I will have to google a bit to understand exactly what it is. – Marc Dec 22 '10 at 14:20
As for "I read your question again and realised that you are not using the linestring the trac forms but each vertex." I guess I could also test for intersection from the linestring object, but then I will need to extract the intersecting part. But that may be faster, you are absolutely right ! – Marc Dec 22 '10 at 14:27
About spatial indexes you should google gist and multidimensional indexes. The idea is to build a multidimmensional index of the bounding boxes. In db-evironment the planner then decides if it is worth the effort to search the index first to find the intersecting bounding boxes before doing the real intersection test. – Nicklas Avén Dec 22 '10 at 14:35