2

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:

Example dataset

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.

Wrongly dissolved polygons

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).

1
  • Please decide whether it is Geopandas, PyQGIS or another spatial library that you wish to make the focus of your question.
    – PolyGeo
    Commented Feb 16, 2023 at 7:33

1 Answer 1

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:

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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