[gdal-dev] OGR Select Features by Location...

[Jose Luis Gomez Dans]

Hi,


> I think that if you built GDAL with GEOS support, you can call on the GEOS
> methods from directly within your ogr objects.  (Not sure I'm saying that
> right)  There is an example here -
> http://mapserver.gis.umn.edu/community/conferences/MUM3/workshop/python/geometriespdf
> 
> Code in the example looks something like this:
> 
> #! /usr/bin/python
> 
> import ogr
> 
> r1 = {'minx': -5.0, 'miny': 0.0, 'maxx': 5.0, 'maxy':10.0}
> r2 = {'minx': 0.0, 'miny': -5.0, 'maxx': 10.0, 'maxy':5.0}
> 
> template = 'POLYGON ((%(minx)f %(miny)f, %(minx)f %(maxy)f, %(maxx)f
> %(maxy)
> f, %(maxx)f %(miny)f, %(minx)f %(miny)f))'
> 
> w1 = template % r1
> w2 = template % r2
> 
> g1 = ogr.CreateGeometryFromWkt(w1)
> g2 = ogr.CreateGeometryFromWkt(w2)
> 
> inter = g1.Intersection(g2)   # <---------------GEOS Intersection command?
> 
> print inter
> 

While this will work, I doubt it will cope with multi-million points. Say that you know the extent of your polygon (min and max x and y coordinates), you could extract the points within that region "easily" with OGR and python:
import ogr
r1 = {'minx': -5.0 , 'miny': 0.0, 'maxx': 5.0, 'maxy':10.0}

template = 'POLYGON ((%(minx)f %(miny)f, %(minx)f %(maxy)f, %(maxx)f %(maxy) f, %(maxx)f %(miny)f, %(minx)f %(miny)f))'

w1 = template % r1

g1 = ogr.CreateGeometryFromWkt(w1)
s = ogr.Open (<fname>)
L = s.GetLayer(0) #Or layer whatever
L.SetSpatialFilter ( g1 )
#If it is a square region,
#L.SetSpatialFilterRect(r1['minx'], r1['miny'], r1['maxx'], r1['maxy'])


Clearly, you can use any geometry, not just a square. You can have a shapefile with your polygon geometry, read the first feature, and use that as your SpatialFilter geometry, but I don't know how fast this is. If you decide to go this way, it would be interesting to know how well it performs.

Jose
-- 
GMX FreeMail: 1 GB Postfach, 5 E-Mail-Adressen, 10 Free SMS.
Alle Infos und kostenlose Anmeldung: http://www.gmx.net/de/go/freemail