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

[Roger André]

Hmm, performance is a good point.  I thought that SetSpatialFilter only
returned the intersection of the 2 envelopes though?  Will that accurately
capture only points that are within an odd-shaped polygon.

Roger
--

On Jan 23, 2008 9:41 AM, Jose Luis Gomez Dans <josegomez at gmx.net> wrote:

> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20080123/f2527636/attachment.html