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

[Martel, Christian]

Frank,

Thanks. This is the hint I was looking for.

I read in the doc about SetSpatialFilter:
"Currently this test is may be inaccurately implemented, but it is
guaranteed that all features who's envelope (as returned by
OGRGeometry::getEnvelope()) overlaps the envelope of the spatial filter
will be returned. This can result in more shapes being returned that
should strictly be the case."

It's ok because I only want a bbox comparison. Do you plan to update the
function to do an intersection test with the actual geometry ? If so I
should try something else to ensure my code will always return the same
result.

Thanks,
Christian

-----Original Message-----
From: Frank Warmerdam [mailto:warmerdam at pobox.com] 
Sent: 6 juin 2007 14:11
To: Martel, Christian
Cc: gdal-dev at lists.maptools.org
Subject: Re: [Gdal-dev] spatial intersection of ogr geometry and
shapefile features

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