Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am using the GEOS/OGR library in Python. I have files full of geometries that overlap, and I want to convert these to files of geometries that don't overlap.

The way I have been approaching this is to just read in all of the geometries from a file, and then looping over each geometry and taking the Union() of that with all of the previous geometries. So something like:

def dedup(geometries):
   """Return a geometry that is the union of all geometries."""
   if not geometries:  return None
   current_union = geometries[0].Clone()
   for g in geometries:
      current_union = current_union.Union(g).Clone()
   return current_union

Most of the time this appears to work, but occasionally Union() will return None. My debugging reveals that under these cases:

  • IsValid() returns true for both g and current_union
  • Intersects() return true between g and current_union
  • GetGeometryCount() for current_union typically returns a number over 100.
  • I suspect the error might be some kind of memory issue because it seems like the error "moves" depending on context. Ie if I say run this on files 1-100, it might crash on file 20, but if I say run it on files 20-100, it might crash on file 40. I never use Destroy() on any geometries. Task Manager does not show me running out of memory though.

Anyway, maybe there is an easier way to do this? Or else does anyone have any idea why it is crashing here?

Thanks much in advance!

share|improve this question

No need to clone any geometries. Also, iterate over geometries[1:] to save one Union().

If osgeo.ogr has something like Shapely's cascaded_union use that instead of accumulating unions in a loop.

share|improve this answer
There is a UnionCascaded function apparently but I don't think it works. There's no documentation for it in the SWIG wrapper and it throws an exception if you pass any arguments into it. I am starting to think there might be a bug in the GEOS Union() function as well. – Ben Jun 3 '12 at 23:01
I've seen reports of union bugs in the GEOS tracker. Find what your version is and look in… to see whether any issue has been closed since or if it looks like one of the open ones. – sgillies Jun 4 '12 at 14:38

(This is an old question, but still is the top GIS.SE result for a search of "Python Union", so I am offering an answer. Hopefully it will be of some use.)

Python osgeo inlcudes a UnionCascaded method that implements the underlying GEOS CascadedUnion method. It is a little strange in that while Union takes two geometries (in form geom1.Union(geom2)), UnionCascaded operates on a geometry collection. According to GEOS documentation, this method is more efficient than performing a union piecewise on individual geometries.

I don't know exactly what you are passing into the function in geometries. If it is a geometry collection, possibly geometries.UnionCascaded() will work. Bugs discussed in comments (by OP Ben and user sgillies) could have been fixed in the last three years. Union will work to combine geometries of different dimension (e.g. polygons and lines), and possibly UnionCascaded will as well, but I'm having a hard time understanding how that would be useful or what you would do with the resulting object. So, as an example, I'm going to assume that geometries is a layer with features of the same dimension. I think the easiest way to operationalize what you want to do is to collect your geometries into a multi and then call UnionCascaded on that object. For now, let's assume we are working with polygons.

def dedup(geometries):
    """Return a geometry that is the union of all geometries."""
    if not geometries:  return None
    multi  = ogr.Geometry(ogr.wkbMultiPolygon)
    for g in geometries:
    return multi.UnionCascaded()

Possibly (I haven't tested) if you replace ogr.wkbMultiPolygon with ogr.wkbUnknown, this would work for lines or points. But possibly the geometry type of the output would also be unknown, and that may make it difficult to work with.

share|improve this answer

Your Answer


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.