All, I think I found what I should have been using. pick a feature from one layer, use its geometry as an argument of SetSpatialFilter method of the other layer. Seems to be yielding same answer quicker...
Francis, Thank you for suggestion. I looked at the Shapely documentation http://gispython.org/shapely/docs/1.0/manual.html . It wasn't obvious to me advantage of the library over ogr, except that it looks more pythonic. I even didn't find mention of layer vs. feature, and to solve my problem it is a step back to me? On Nov 5, 4:54 pm, Francis Markham <[email protected]> wrote: > If you're using Python, I'd recommend checking out the shapely package > (http://pypi.python.org/pypi/Shapely): > > You can load your shapefile using ogr, and then import the geometry to > shapely where you can do the geometry predicate calculations: > > >>> import ogr > >>> from shapely.wkb import loads > >>> source = ogr.Open("/tmp/world_borders.shp") > >>> borders = source.GetLayerByName("world_borders") > >>> feature = borders.GetNextFeature() > >>> loads(feature.GetGeometryRef().ExportToWkb()) > > Good luck, > > Francis > > On 6 November 2010 06:08, yosuke kimura <[email protected]> wrote: > > > All, > > > I am trying to calculate intersection of two polygon layers with > > polygons. I searched threads, found following thread, and what I got > > was that OGR/GEOS provide me to intersection of two features, but I > > have to write something in order to perform this operation for two > > layers with multiple features. > > > Open source vector geoprocessing libraries? > >http://groups.google.com/group/gdal/browse_thread/thread/7b5a4f461a4d... > > > Intersect of two shapefiles > >http://groups.google.com/group/gdal/browse_thread/thread/4b178c776aec... > > > My questions are: > > > 1. Did I miss something? Is there development of library which takes > > two polygon layers to calculate intersection or other basic geometry > > operations that i find in OGR? > > > 2. Assuming no to above, I got started writing code using python > > binding. Not having too much experience, I codes to go through all > > possible pair (complete bipartite), test with intersect, then run > > intersection if it returns true. I am wondering there is faster way > > to do this. I am thinking something like get envelope polygon of each > > feature first, and then test intersection of them and then perform > > intersect operation only when first test passes (i assume that it is > > easier to test intersect if there are only handful of vertices as > > oppose to hundreds of vertices in my case). > > > Right now I am trying with layer of hundreds of thousands of polygons > > vs. thousands of polygon, and it seems to be taking 10s of hours to > > finish in my environment. > > > Thank you, > > > Yosuke Kimura > > _______________________________________________ > > gdal-dev mailing list > > [email protected] > >http://lists.osgeo.org/mailman/listinfo/gdal-dev > > > > _______________________________________________ > gdal-dev mailing list > [email protected]http://lists.osgeo.org/mailman/listinfo/gdal-dev _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
