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

[estrada.adam at gmail.com]

All,
    How important is it to create the index for the intersection
function I am trying to perform.  We are talking over 15 million LIDAR
points already loaded in to my PostgreSQL DB for a single dataset.  I
read something that said that each set of 1 million points takes an
hour to process when creating the index.  I have 60 datasets that are
similar in size...Arghhh...I suppose that I will try it with and then
without the indexes and report back.

AWE

On Jan 29, 1:59 am, Joe Healy <j... at omc-international.com.au> wrote:
> 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/pyth...
>
> >> 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
>
> _______________________________________________
> gdal-dev mailing list
> gdal-... at lists.osgeo.orghttp://lists.osgeo.org/mailman/listinfo/gdal-dev- Hide quoted text -
>
> - Show quoted text -