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.<br><br>Roger<br>
--<br><br><div class="gmail_quote">On Jan 23, 2008 9:41 AM, Jose Luis Gomez Dans <<a href="mailto:josegomez@gmx.net">josegomez@gmx.net</a>> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Hi,<br><div class="Ih2E3d"><br><br>> I think that if you built GDAL with GEOS support, you can call on the GEOS<br>> methods from directly within your ogr objects. (Not sure I'm saying that<br>> right) There is an example here -
<br>> <a href="http://mapserver.gis.umn.edu/community/conferences/MUM3/workshop/python/geometriespdf" target="_blank">http://mapserver.gis.umn.edu/community/conferences/MUM3/workshop/python/geometriespdf</a><br>><br>
> Code in the example looks something like this:<br>><br>> #! /usr/bin/python<br>><br>> import ogr<br>><br>> r1 = {'minx': -5.0, 'miny': 0.0, 'maxx': 5.0, 'maxy':10.0}<br>
> r2 = {'minx': 0.0, 'miny': -5.0, 'maxx': 10.0, 'maxy':5.0}<br>><br>> template = 'POLYGON ((%(minx)f %(miny)f, %(minx)f %(maxy)f, %(maxx)f<br>> %(maxy)<br>> f, %(maxx)f %(miny)f, %(minx)f %(miny)f))'
<br>><br>> w1 = template % r1<br>> w2 = template % r2<br>><br>> g1 = ogr.CreateGeometryFromWkt(w1)<br>> g2 = ogr.CreateGeometryFromWkt(w2)<br>><br>> inter = g1.Intersection(g2) # <---------------GEOS Intersection command?
<br>><br>> print inter<br>><br><br></div>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:
<br><div class="Ih2E3d">import ogr<br>r1 = {'minx': -5.0 , 'miny': 0.0, 'maxx': 5.0, 'maxy':10.0}<br><br></div><div class="Ih2E3d">template = 'POLYGON ((%(minx)f %(miny)f, %(minx)f %(maxy)f, %(maxx)f %(maxy) f, %(maxx)f %(miny)f, %(minx)f %(miny)f))'
<br><br>w1 = template % r1<br><br></div>g1 = ogr.CreateGeometryFromWkt(w1)<br>s = ogr.Open (<fname>)<br>L = s.GetLayer(0) #Or layer whatever<br>L.SetSpatialFilter ( g1 )<br>#If it is a square region,<br>#L.SetSpatialFilterRect(r1['minx'], r1['miny'], r1['maxx'], r1['maxy'])
<br><br><br>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.
<br><br>Jose<br><font color="#888888">--<br>GMX FreeMail: 1 GB Postfach, 5 E-Mail-Adressen, 10 Free SMS.<br>Alle Infos und kostenlose Anmeldung: <a href="http://www.gmx.net/de/go/freemail" target="_blank">http://www.gmx.net/de/go/freemail
</a><br></font></blockquote></div><br>