Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I'm working with a large dataset, and I need to be able to combine up to 400,000 polygons. I have been using the QgsGeometry::combine(QgsGeometry* geometry) function and simply looping over all geometries in an interative fashion like this:

union_geoms = geometries[0]
for geometry in geometries:
  union_geoms = union_geoms.combine(geometry)

This works fine for most of my data, (the largest list of polygons I have gotten it to complete is ~30,000.) But I have two sets of ~200,000 and one of ~400,000 and I can't get this to complete in a timely fashion (I have let it run for days to no avail). So I started looking for a better method to combine these geometries, and I found geos::operation::geounion::CascadedPolygonUnion This algorithm, described here, is well suited to my application as the polygons of interest contain a lot of adjacency.

I found a couple references to implementations of this in Qgis, here in 2010 and here in 2013, but it doesn't appear to have made it into the API (unless I'm looking in the wrong place!)

So my question is how would I go about implementing this type of algorithm, or more simply how can I use the existing implementation in GEOS from my own (python) scripts? From the 2013 link,

you can use QgsGeometry::asGeos() in the plugin instead of wkt-export / wkt-import. For cascaded union, there is the method GEOSUnaryUnion(const GEOSGeometry* g1) in the c-api. The argument seems to be a geometry collection, which can be created with GEOSGeom_createCollection.

But there doesn't seem to by python bindings for the asGeos() method.

Anyone have suggestions as to how to proceed?

share|improve this question

1 Answer 1

up vote 3 down vote accepted

Shapely is one of the Geos Python bindings and has cascaded_union and unary_union implemented since versions 2.16 (GEOSCascadedUnion is deprecated since GEOS version 3.2.+ and GEOSUnaryUnion must be used instead: it can operate on different geometry types, not only polygons as is the case for the older cascaded unions).

Convert QGIS geometries to Shapely geometries in the Python console is very easy:

from shapely.wkb import loads
from shapely.ops import cascaded_union, unary_union
# transform QGIS geometries to a list of shapely geometries
geoms =[loads(feature.geometry().asWkb()) for feature in layer]
# cascaded union
cu= cascaded_union(geoms)
cu
<shapely.geometry.polygon.Polygon object at 0x1257d9750>
# unary union
un = unary_union(geoms)
un
<shapely.geometry.polygon.Polygon object at 0x1257d9c90>
# convert to Qgis geometry
geom = QgsGeometry.fromWkt(un.wkt)
geom
<qgis.core.QgsGeometry object at 0x1247c7dd0>

You can even do it without QGIS with modules like Fiona to read and write a Shapefile layer.

share|improve this answer
    
Can you explain the difference between cascaded_union and unary_union? –  Nelz11 Jun 19 '13 at 21:40
    
Also, how would I then use this shapely geometry to create a new feature? To I need to convert back to wkb and load this into a features geometry? –  Nelz11 Jun 19 '13 at 22:17
    
Nevermind on the converting back part...I see that in the last couple lines of your answer! –  Nelz11 Jun 19 '13 at 22:21
    
unary_union operates on any type of geometry, not just polygons. –  gene Jun 20 '13 at 7:23

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.