[Gdal-dev] spatial intersection of ogr geometry and shapefile features

[Frank Warmerdam]

Martel, Christian wrote:
> Hi,
> 
> I have an OGR geometry and a shapefile of all worldwide countries. I
> need to find which countries intersect my geometry (bbox of country is
> fine). I'm working in Python with FWTools. I found some hints in the
> list but I need more help.
> 
> What should I try ?
> - Create a qix index file and work with shapefile ? If so, how to
> compare geometries with DataSource.ExecuteSQL(...) command ? How to use
> the index for that comparison ?
> - In a loop, build ogr geometry of each country and compare it with the
> current geometry (not very efficient I think) ?
> 
> I have experience with gdal/ogr with python in fwtool and I know how to
> write/read shapefiles.

Christian,

Before reading the features from the shapefile layer of countries,
first call SetSpatialFilter() passing in the geometry you want to
intersect with.  The spatial filter will be used to restrict the set
of features (countries) you get back.

If speed is important (hard to imagine with a list of countries) then
you could also build a spatial index on the country shapefile but that
is necessary to make things work.

The test script at:
   http://trac.osgeo.org/gdal/browser/trunk/autotest/ogr/ogr_basic_test.py#L78

should be a useful example-in-action of spatial filtering from python.

If you want to do a more details geometric test after preliminary spatial
filtering you might want to review the GEOS based examples at:

   http://trac.osgeo.org/gdal/browser/trunk/autotest/ogr/ogr_geos.py

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org