I'm trying to perform lots of unions in GeoPandas. I've got around fifty sets of between 7-12 GeoDataFrames of the same place at different times that I'm trying to union into one layer per set. Right now, about 20% of the unions are failing due to non-noded intersection exceptions. I'm not sure why snapping the layers doesn't solve this issue. The union looks like this:

# Snap the add layer to the base
flat = flatten(base)
add.geometry  = [snap(g, flat, 0.0001) for g in add.geometry]
base.geometry = [snap(g, flat, 0.0001) for g in base.geometry]

# Ensure the geometries are still 
add  = ensureValid(add)
base = ensureValid(base)

# Union the add layer to the base
base= gpd.overlay(add, base, "union")

Here are the helper functions:

def flatten(gdf):

    # Try performing a unary union
        flat = unary_union(gdf)


        # Try performing a unary union with buffers
            flat = unary_union(gdf.buffer(0.00001).buffer(-0.00001))

        # Break the MultiPolygons down into LineStrings and do the unary union on those
            lines = []
            for f in mapping(ensureValid(gdf).geometry)["features"]:
                 coords = json.loads(json.dumps(f["geometry"]["coordinates"]))
                 cpairs = []
                 for a in range(len(coords)):
                     for b in range(len(coords[a][0])):
            flat = unary_union(lines)

    # Return the result
    return flat


def ensureValid(gdf, geom_type=MultiPolygon):

    # See if there are invalid geometries
    if len(gdf[~gdf.is_valid]):

        # Use a zero buffer to fix problematic intersections
        gdf.geometry = gdf.geometry.apply(lambda g: g.buffer(0.0) if not g.is_valid else g)

        # Remove remaining invalid geometries because I have no idea how to fix them
        gdf = gdf.loc[gdf.is_valid & ~gdf.is_empty]

    # Ensure that every geometry has the correct type
    gdf.geometry = gdf.geometry.apply(lambda g: geom_type([g]) if type(g) != geom_type else g)

    return gdf

Relevant imports:

import geopandas as gpd
import json, shapely
from shapely.geometry import mapping, LineString
from shapely.geometry.multipolygon import MultiPolygon
from shapely.ops import snap, unary_union


    base = gpd.overlay(add, base, "union")

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 371, in overlay
    result = _overlay_union(df1, df2)

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 299, in _overlay_union
    dfsym = _overlay_symmetric_diff(df1, df2)

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 273, in _overlay_symmetric_diff
    dfdiff1 = _overlay_difference(df1, df2)

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 260, in _overlay_difference
    [geom] + list(df2.geometry.iloc[neighbours]))

  File "C:\ProgramData\Anaconda3\lib\site-packages\geopandas\tools\overlay.py", line 259, in <lambda>
    new = reduce(lambda x, y: x.difference(y).buffer(0),

  File "C:\ProgramData\Anaconda3\lib\site-packages\shapely\geometry\base.py", line 616, in difference
    return geom_factory(self.impl['difference'](self, other))

  File "C:\ProgramData\Anaconda3\lib\site-packages\shapely\topology.py", line 70, in __call__
    self._check_topology(err, this, other)

  File "C:\ProgramData\Anaconda3\lib\site-packages\shapely\topology.py", line 39, in _check_topology
    raise err

TopologicalError: This operation could not be performed. Reason: unknown

What is wrong here?

  • have you tried makeValid? – Elio Diaz Dec 9 '19 at 19:24
  • @ElioDiaz Is that possible to do from inside my Python script? – Kung Fu Howie Dec 9 '19 at 19:55
  • reading on makeValid I found it is included (I've used it) in pyqgis dir(QgsVectorLayer); it is a function from postgis, also included in R (I've also used that one); I thought it would be available for geopandas or shapely, there's a function on the web, github.com/ftwillms/makevalid which I haven't tried – Elio Diaz Dec 9 '19 at 21:06
  • @ElioDiaz I've now tried makeValid and it hasn't solved the issue. The problem doesn't appear to be with either gdf before the overlay, but rather from a non-noded intersection appearing during the Difference portion of the union. I posted the traceback in the OP. – Kung Fu Howie Dec 14 '19 at 11:59

I figured it out! There were two issues.

  1. The geometries I was trying to put in were not all valid. Huge thanks to @ElioDiaz for finding a version of ST_MakeValid that I could easily include in my script. While the buffer(0) trick works sometimes, the makevalid git repository is a lot more effective. I still had to drop a bunch of invalid geometries, but I think that it's fine for my purposes.

  2. Some of the coordinates created after the unions were not rounded to the correct precision. Although I had rounded them to 1m (0.00001) before, any points that were created during the union were not rounded. Rounding the coordinates again solved the issue.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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