[gdal-dev] Re: Intersection of polygon layers

[yosuke kimura]

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 <u2546... at anu.edu.au> 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 <yosukesa... at gmail.com> 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
> > gdal-... at lists.osgeo.org
> >http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-... at lists.osgeo.orghttp://lists.osgeo.org/mailman/listinfo/gdal-dev