Search code examples
pythongeojsonshapely

Delete polygons that have one or more side parts in common


I am trying to solve a particular case of comparison of polygons to others. I have five polygons distributed as in the figure below. The black polygon is the one with the largest area.

There may be other similar cases, the main rule is to remove the smallest polygons among all those that have one or more side portions in common.

The data for this case are in a GeoJson file as follows:

{"type":"FeatureCollection","features":[
    {"type":"Feature","properties":{"id":1},"geometry":{"type":"Polygon","coordinates":[[[3.4545135498046875,45.533288879467456],[3.4960556030273433,45.533288879467456],[3.4960556030273433,45.57055337226086],[3.4545135498046875,45.57055337226086],[3.4545135498046875,45.533288879467456]]]}},
    {"type":"Feature","properties":{"id":2},"geometry":{"type":"Polygon","coordinates":[[[3.4545135498046875,45.52917023833511],[3.4960556030273433,45.52917023833511],[3.4960556030273433,45.53891018749409],[3.4545135498046875,45.53891018749409],[3.4545135498046875,45.52917023833511]]]}},
    {"type":"Feature","properties":{"id":3},"geometry":{"type":"Polygon","coordinates":[[[3.4845542907714844,45.5298015824607],[3.5159683227539062,45.5298015824607],[3.5159683227539062,45.543388795387294],[3.4845542907714844,45.543388795387294],[3.4845542907714844,45.5298015824607]]]}},
    {"type":"Feature","properties":{"id":4},"geometry":{"type":"Polygon","coordinates":[[[3.465328216552734,45.542667432984864],[3.4735679626464844,45.542667432984864],[3.4735679626464844,45.5478369923404],[3.465328216552734,45.5478369923404],[3.465328216552734,45.542667432984864]]]}},
    {"type":"Feature","properties":{"id":5},"geometry":{"type":"Polygon","coordinates":[[[3.4545138850808144,45.56799974017372],[3.4588050842285156,45.56799974017372],[3.4588050842285156,45.57055290285386],[3.4545138850808144,45.57055290285386],[3.4545138850808144,45.56799974017372]]]}}]}

Is there a solution to delete only the two blue polygons(id 2 and 5)? In python.

By transforming the Polygons into LineString one could look if a Linestring is a portion of another Linestring ? But I don't see how to do it. Or maybe looking to see if the LineString of the black and blue polygons have more than two points in common ? But we can't convert a LineString into more than two points.

enter image description here


Solution

  • The following approach may work for you using shared_paths which correctly calls out the path overlap between polygons 1, 2 and 5:

    import json
    import shapely as sh
    import shapely.ops as ops
    import shapely.geometry as geo
    
    with open('./test.json') as f:
      features = json.load(f)['features']
    
    for f1 in features:
      for f2 in features:
        id1 = f1['properties']['id']
        id2 = f2['properties']['id']
        if int(id1) > int(id2):
          s1 = geo.shape(f1['geometry'])
          s2 = geo.shape(f2['geometry'])
          coll = ops.shared_paths(s1.boundary, s2.boundary)
          if not coll.is_empty:
            print(f"{id1} and {id2} have shared path")
            # update your feature collection etc
    

    I had to reduce the precision to 5 decimal places in your feature geometry for this to work as initially it only detects the overlap between polygon 1 and 2. The shared corner between polygon 1 and 5 is slightly out in your input FeatureCollection:

    {
      "type": "FeatureCollection",
      "features": [{
          "type": "Feature",
          "properties": {
            "id": 1
          },
          "geometry": {
            "type": "Polygon",
            "coordinates": [
              [
                [3.45451, 45.53328],
                [3.49605, 45.53328],
                [3.49605, 45.57055],
                [3.45451, 45.57055],
                [3.45451, 45.53328]
              ]
            ]
          }
        },
        {
          "type": "Feature",
          "properties": {
            "id": 2
          },
          "geometry": {
            "type": "Polygon",
            "coordinates": [
              [
                [3.45451, 45.52917],
                [3.49605, 45.52917],
                [3.49605, 45.53891],
                [3.45451, 45.53891],
                [3.45451, 45.52917]
              ]
            ]
          }
        },
        {
          "type": "Feature",
          "properties": {
            "id": 3
          },
          "geometry": {
            "type": "Polygon",
            "coordinates": [
              [
                [3.48455, 45.52980],
                [3.51596, 45.52980],
                [3.51596, 45.54338],
                [3.48455, 45.54338],
                [3.48455, 45.52980]
              ]
            ]
          }
        },
        {
          "type": "Feature",
          "properties": {
            "id": 4
          },
          "geometry": {
            "type": "Polygon",
            "coordinates": [
              [
                [3.465328, 45.54266],
                [3.473567, 45.54266],
                [3.473567, 45.54783],
                [3.465328, 45.54783],
                [3.465328, 45.54266]
              ]
            ]
          }
        },
        {
          "type": "Feature",
          "properties": {
            "id": 5
          },
          "geometry": {
            "type": "Polygon",
            "coordinates": [
              [
                [3.454513, 45.56799],
                [3.458805, 45.56799],
                [3.458805, 45.57055],
                [3.454513, 45.57055],
                [3.454513, 45.56799]
              ]
            ]
          }
        }
      ]
    }