Taking Union of several geometries in GEOS Python?

Question

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?

Answer 1

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.

Answer 2

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:
        multi.AddGeometry(g.geometry())
    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.

Answer 3

I had a similar problem and solved it like this, taking all the previous answers into account:

import os
from osgeo import ogr

def createDS(ds_name, ds_format, geom_type, srs, overwrite=False):
    drv = ogr.GetDriverByName(ds_format)
    if os.path.exists(ds_name) and overwrite is True:
        deleteDS(ds_name)
    ds = drv.CreateDataSource(ds_name)
    lyr_name = os.path.splitext(os.path.basename(ds_name))[0]
    lyr = ds.CreateLayer(lyr_name, srs, geom_type)
    return ds, lyr

def dissolve(input, output, multipoly=False, overwrite=False):
    ds = ogr.Open(input)
    lyr = ds.GetLayer()
    out_ds, out_lyr = createDS(output, ds.GetDriver().GetName(), lyr.GetGeomType(), lyr.GetSpatialRef(), overwrite)
    defn = out_lyr.GetLayerDefn()
    multi = ogr.Geometry(ogr.wkbMultiPolygon)
    for feat in lyr:
        if feat.geometry():
            feat.geometry().CloseRings() # this copies the first point to the end
            wkt = feat.geometry().ExportToWkt()
            multi.AddGeometryDirectly(ogr.CreateGeometryFromWkt(wkt))
    union = multi.UnionCascaded()
    if multipoly is False:
        for geom in union:
            poly = ogr.CreateGeometryFromWkb(geom.ExportToWkb())
            feat = ogr.Feature(defn)
            feat.SetGeometry(poly)
            out_lyr.CreateFeature(feat)
    else:
        out_feat = ogr.Feature(defn)
        out_feat.SetGeometry(union)
        out_lyr.CreateFeature(out_feat)
        out_ds.Destroy()
    ds.Destroy()
    return True

UnionCascaded requires MultiPolygon as a geometry type, which is why I implemented the option to re-create single polyons. You could also use ogr2ogr from the command line with option -explodecollections:

ogr2ogr -f "ESRI Shapefile" -explodecollections dissolved.shp input.shp -dialect sqlite -sql "select ST_union(Geometry) from input"

Answer 4

I have found that the Union() function will return None if it is passed a polygon, which looks good otherwise, but is actually invalid because it is missing the closing point at the end of the polygon.

Make sure to duplicate the first point's coordinates as the last point of the polygon.