Dissolving overlapping multipolygon geometries into new multipolygons and keeping non-overlapping features intact

Question

I want to merge (multi)polygon geometries any time they overlap partly or fully. For example the Dissolve function in QGIS does not fit my case, since I also want to:

1. Include all parts of the multipolygons that do not overlap with another geometry in the newly dissolved feature's geometry.
2. Keep any multipolygon geometries that do not overlap at all as a single feature.

I'll illustrate my point with a toy geojson dataset:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [[0,0], [0,1], [1,1], [1,0], [0,0]]
          ],
          [
            [[1,1], [1,2], [2,2], [2,1], [1,1]]
          ]
        ]
      },
      "properties": {
        "name": "A"
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [[4,2], [4,3], [5,3], [5,2], [4,2]]
          ],
          [
            [[2.5,2], [2,3], [3.5,3], [3,2], [2.5,2]]
          ]
        ]
      },
      "properties": {
        "name": "B"
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [[2,4], [2,5], [3,5], [3,4], [2,4]]
          ],
          [
            [[2,2], [2,3], [3,3], [3,2], [2,2]]
          ]
        ]
      },
      "properties": {
        "name": "C"
      }
    },
        {
      "type": "Feature",
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [[2.5,4.1], [2.1,4.75], [2.75,4.75], [2.9,4.25], [2.5,4.1]]
          ]
        ]
      },
      "properties": {
        "name": "D"
      }
        }
  ]
}

Which looks like:

In this case, B and C overlap and D is within C. A does not overlap with the others.

Using the method proposed here:

import geopandas as gpd

with open('data.geojson', 'r') as f:
    data = json.load(f)
    
gdf = gpd.GeoDataFrame.from_features(data['features'])

gpd.GeoDataFrame(
  geometry=[gdf.unary_union]).explode(
  index_parts=False).reset_index(
  drop=True).plot(categorical=True, cmap='Set1', alpha=0.66)

This dissolves the overlaps correctly, but it also breaks up the geometrically separate parts of the multipolygon into single polygons. Using Dissolve with Keep disjoint features separate in QGIS has the same result.

The expected output of the dissolving process in the example would be a layer with two features, one consisting of multipolygon A and the other the newly merged geometries from B-C-D. Attributes do not need to be retained.

How should I move forward with this, using for example the Geopandas script above as a basis?

If it matters, the actual data I would be using this for is polygons from the World Database of Protected Areas (WDPA).

Answer 1

I think the following code gives the result you mentioned as expected output. The inline comments explain what is happening:

import json
from pathlib import Path

import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd

url = "https://raw.githubusercontent.com/theroggy/pysnippets/07d778d3ac149b3ba2be22735d0a54bf3d52a6d7/pysnippets/stackoverflow_questions/2023-02-15_dissolve_intersecting.geojson"
gdf = gpd.read_file(url)

# Calculate intersections within the layer
intersection_gdf = gdf.overlay(gdf, how="intersection", keep_geom_type=True)
intersection_gdf = intersection_gdf.loc[
    intersection_gdf.name_1 != intersection_gdf.name_2
]

# The features to dissolve are the intersecting ones, excluding the self-intersections
to_dissolve_gdf = gdf.loc[
    gdf.name.isin(intersection_gdf.name_1) | gdf.name.isin(intersection_gdf.name_2)
]

# Other features should not be dissolved
no_dissolve_gdf = gdf.loc[~gdf.index.isin(to_dissolve_gdf.index)]

# Compile + plot the result
result_gdf = pd.concat([to_dissolve_gdf.dissolve(), no_dissolve_gdf])
result_gdf.plot(categorical=True, cmap="Set1", alpha=0.66)
plt.show()

Result: