[gdal-dev] Intersection of polygon layers

Francis Markham u2546226 at anu.edu.au
Fri Nov 5 17:54:32 EDT 2010

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,


On 6 November 2010 06:08, yosuke kimura <yosukesabai 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/7b5a4f461a4db319/
> Intersect of two shapefiles
> http://groups.google.com/group/gdal/browse_thread/thread/4b178c776aecf425/
> 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-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20101106/63aeb375/attachment.html

More information about the gdal-dev mailing list