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

[Joe Healy]

We use postgis for this sort of work. It is orders of magnitude faster 
than ArcGIS.

 From memory the largest we have dealt with so far was in the 10's of 
millions of points. It handles this very well.

Depending on the size of the dataset, we find loading the data using the 
copy command in psql, creating indexes and then querying is fast enough 
that people are happy to reload the data etc. each time they perform the 
join. For the larger data sets, loading the data and building the 
indexes once is a big gain.

Hope this helps,

Joe


Jose Luis Gomez Dans 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
>