Merge pull request #425 from GIScience/bug/#424_incorrect_geometry_mu…

…ltipolygon

Incorrect geometry for Multipolygons with shells which share a single point

CHANGELOG.md

@@ -17,6 +17,7 @@ Changelog
 
 ### bugfixes
 
+* fix building valid geometries for multipolygons with shells which share a single point ([#424])
 * change geometry filters to be based on full (unclipped) geometries ([#433])
 * make sure area computation never returns negative results (instead zero is returned for the invalid geometries which previously resulted in negative values) ([#438])
 
@@ -33,6 +34,7 @@ Changelog
 
 [#260]: https://github.com/GIScience/oshdb/issues/260
 [#419]: https://github.com/GIScience/oshdb/pull/419
+[#424]: https://github.com/GIScience/oshdb/pull/424
 [#433]: https://github.com/GIScience/oshdb/issues/433
 [#438]: https://github.com/GIScience/oshdb/pull/438
 [#441]: https://github.com/GIScience/oshdb/pull/441
@@ -45,6 +47,7 @@ Changelog
 [#459]: https://github.com/GIScience/oshdb/pull/459
 
 
+
 ## 0.7.2
 
 ### bugfixes

oshdb-util/src/main/java/org/heigit/ohsome/oshdb/util/geometry/OSHDBGeometryBuilder.java

@@ -251,16 +251,21 @@ private static Geometry getMultiPolygonGeometry(
         timestamp
     );
 
-    // construct rings from polygons
-    List<LinearRing> outerRings = OSHDBGeometryBuilder.buildRings(outerLines).stream()
-        .map(ring -> geometryFactory.createLinearRing(
-            ring.stream().map(node -> new Coordinate(node.getLongitude(), node.getLatitude()))
-                .toArray(Coordinate[]::new)))
-        .collect(Collectors.toList());
+    // construct inner and outer rings
+    List<LinkedList<OSMNode>> outerRingsNodes = OSHDBGeometryBuilder.buildRings(outerLines);
     List<LinkedList<OSMNode>> innerRingsNodes = OSHDBGeometryBuilder.buildRings(innerLines);
+    // check if there are any pinched off sections in outer rings
+    splitPinchedRings(outerRingsNodes, innerRingsNodes, geometryFactory);
     // check if there are any touching inner/outer rings, merge any
     mergeTouchingRings(innerRingsNodes);
     // create JTS rings for non-degenerate rings only
+
+    List<LinearRing> outerRings = outerRingsNodes.stream()
+        .filter(ring -> ring.size() >= LinearRing.MINIMUM_VALID_SIZE)
+        .map(ring -> geometryFactory.createLinearRing(
+            ring.stream().map(node -> new Coordinate(node.getLongitude(), node.getLatitude()))
+                .toArray(Coordinate[]::new)))
+        .collect(Collectors.toList());
     List<LinearRing> innerRings = innerRingsNodes.stream()
         .filter(ring -> ring.size() >= LinearRing.MINIMUM_VALID_SIZE)
         .map(ring -> geometryFactory.createLinearRing(
@@ -332,7 +337,7 @@ private static Stream<Stream<OSMNode>> waysToLines(
    *   Touching rings are defined as rings which share at least one segment (a segment is formed by
    *   two consecutive ring nodes, regardless of their order). An example is:
    *   [r1 = (A,B,C,D,E,F,A); r2 = (X,Y,B,C,D,E,X)].
-   *   The result would be: [r1 = (B,A,F,E,X,Y,B)] "or any equivalent representation of this ring
+   *   The result would be: [r1 = (B,A,F,E,X,Y,B)] "or any equivalent representation of this ring"
    * </p>
    *
    * <pre>
@@ -343,7 +348,7 @@ private static Stream<Stream<OSMNode>> waysToLines(
    * A----B--Y       A----B--Y
    * </pre>
    *
-   * @param ringsNodes a collection of node-lists, each forming a ring closed linestring
+   * @param ringsNodes a collection of node-lists, each forming a ring (i.e. a closed linestring)
    */
   private static void mergeTouchingRings(Collection<LinkedList<OSMNode>> ringsNodes) {
     // ringSegments will hold a reference of which ring a particular segment is part of.
@@ -400,6 +405,128 @@ private static void mergeTouchingRings(Collection<LinkedList<OSMNode>> ringsNode
     }
   }
 
+
+
+  /**
+   * Search and split self-intersecting/pinched/figure-8 rings.
+   *
+   * <p>Attention: modifies the input data, such that there are no more figure-8 rings.</p>
+   *
+   * <p>
+   *   A pinched ring forms a figure-8 configuration where the ring touches itself in a single
+   *   point. An example is: [r = (A,B,C,D,E,F,C,G,A)].
+   *   The result would be: [r1 = (C,D,E,G,C); r2 = (A,B,C,G,A)].
+   * </p>
+   *
+   * <pre>
+   *  A--B
+   *  |  |
+   *  G--C--D
+   *     |  |
+   *     F--E
+   * </pre>
+   *
+   * @param ringsNodes a collection of node-lists, each forming a ring (i.e. a closed linestring)
+   * @param holeRingsNodes a collection where holes formed by "upended" figure-8's should be stored
+   * @param geometryFactory a {@link GeometryFactory} object to create temporary features
+   */
+  private static void splitPinchedRings(
+      Collection<LinkedList<OSMNode>> ringsNodes,
+      Collection<LinkedList<OSMNode>> holeRingsNodes,
+      GeometryFactory geometryFactory
+  ) {
+    Map<Long, Integer> nodeIds = new HashMap<>();
+    Collection<LinkedList<OSMNode>> additionalRings = new LinkedList<>();
+    for (LinkedList<OSMNode> ringNodes : ringsNodes) {
+      var splitRings = splitPinchedRing(ringNodes, nodeIds);
+      if (splitRings != null) {
+        // if self-intersection(s) were found, we need to check whether these are next to or
+        // overlapping each other. to do this, we convert the rings to polygon geometries first
+        splitRings.add(new LinkedList<>(ringNodes));
+        ringNodes.clear();
+        var splitRingsGeoms = splitRings.stream()
+            .map(ring -> {
+              if (ring.size() >= LinearRing.MINIMUM_VALID_SIZE) {
+                return geometryFactory.createPolygon(ring.stream()
+                    .map(node -> new Coordinate(node.getLongitude(), node.getLatitude()))
+                    .toArray(Coordinate[]::new));
+              } else {
+                return geometryFactory.createPolygon();
+              }
+            })
+            .collect(Collectors.toList());
+        // determine which of the rings is "coveredBy" how many of the others
+        var nestingNumbers = Collections.nCopies(splitRingsGeoms.size(), 0)
+            .toArray(new Integer [] {});
+        for (var i = 0; i < splitRingsGeoms.size(); i++) {
+          for (var j = 0; j < splitRingsGeoms.size(); j++) {
+            if (i == j) {
+              continue;
+            }
+            if (splitRingsGeoms.get(i).coveredBy(splitRingsGeoms.get(j))) {
+              nestingNumbers[i]++;
+            }
+          }
+        }
+        // sort result into (additional) rings and holes
+        for (var i = 0; i < splitRingsGeoms.size(); i++) {
+          if (nestingNumbers[i] % 2 == 0) {
+            additionalRings.add(splitRings.get(i));
+          } else {
+            holeRingsNodes.add(splitRings.get(i));
+          }
+        }
+      }
+    }
+    ringsNodes.addAll(additionalRings);
+  }
+
+  /**
+   * Search and split pinched (figure-8) rings.
+   *
+   * @return null if no self-intersection is found,
+   *         otherwise a collection containing additional split-off rings
+   */
+  private static List<LinkedList<OSMNode>> splitPinchedRing(
+      LinkedList<OSMNode> ringNodes,
+      Map<Long, Integer> nodeIds
+  ) {
+    List<LinkedList<OSMNode>> result = null;
+    boolean wasSplittable;
+    do {
+      wasSplittable = false;
+      nodeIds.clear();
+      var currentNodePos = 0;
+      for (OSMNode ringNode : ringNodes) {
+        long nodeId = ringNode.getId();
+        if (nodeIds.containsKey(nodeId)) {
+          // split off ring between previous and current ring position
+          int nodePos = nodeIds.get(nodeId);
+          final var additionalRing = new LinkedList<>(ringNodes.subList(nodePos, currentNodePos + 1));
+          final var remainingRing = new LinkedList<OSMNode>();
+          remainingRing.addAll(ringNodes.subList(0, nodePos));
+          remainingRing.addAll(ringNodes.subList(currentNodePos, ringNodes.size()));
+          wasSplittable = true;
+          // add to results
+          ringNodes.clear();
+          ringNodes.addAll(remainingRing);
+          if (result == null) {
+            result = new ArrayList<>();
+          }
+          result.add(additionalRing);
+          break;
+        }
+        if (currentNodePos > 0) {
+          // don't memorize start node, since it is always repeated at the end of the ring
+          nodeIds.put(nodeId, currentNodePos);
+        }
+        currentNodePos++;
+      }
+      // repeat until the ring doesn't have any more self intersections
+    } while (wasSplittable);
+    return result;
+  }
+
   /**
    * Cut a ring at the given segment.
    *
@@ -533,26 +660,33 @@ private static List<LinkedList<OSMNode>> buildRings(Stream<Stream<OSMNode>> line
             what.removeFirst();
             current.addAll(what);
             waysIterator.remove();
+            lastId = current.getLast().getId();
             joinable = true;
           } else if (firstId == what.getLast().getId()) {
             // start of partial ring matches end of current line
             what.removeLast();
             current.addAll(0, what);
             waysIterator.remove();
+            firstId = current.getFirst().getId();
             joinable = true;
           } else if (lastId == what.getLast().getId()) {
             // end of partial ring matches end of current line
             what.removeLast();
             current.addAll(Lists.reverse(what));
             waysIterator.remove();
+            lastId = current.getLast().getId();
             joinable = true;
           } else if (firstId == what.getFirst().getId()) {
             // start of partial ring matches start of current line
             what.removeFirst();
             current.addAll(0, Lists.reverse(what));
             waysIterator.remove();
+            firstId = current.getFirst().getId();
             joinable = true;
           }
+          if (firstId == lastId) {
+            break;
+          }
         }
         // joinable==false for invalid geometries (dangling way, unclosed ring)
       } while (joinable);

oshdb-util/src/test/java/org/heigit/ohsome/oshdb/util/geometry/OSHDBGeometryBuilderTest.java

@@ -6,9 +6,13 @@
 import static org.junit.jupiter.api.Assertions.assertFalse;
 import static org.junit.jupiter.api.Assertions.assertTrue;
 
+import java.util.function.Consumer;
 import org.heigit.ohsome.oshdb.OSHDBBoundingBox;
 import org.heigit.ohsome.oshdb.OSHDBTimestamp;
+import org.heigit.ohsome.oshdb.osm.OSM;
 import org.heigit.ohsome.oshdb.osm.OSMEntity;
+import org.heigit.ohsome.oshdb.osm.OSMMember;
+import org.heigit.ohsome.oshdb.osm.OSMRelation;
 import org.heigit.ohsome.oshdb.util.geometry.helpers.FakeTagInterpreterAreaAlways;
 import org.heigit.ohsome.oshdb.util.geometry.helpers.FakeTagInterpreterAreaMultipolygonAllOuters;
 import org.heigit.ohsome.oshdb.util.geometry.helpers.FakeTagInterpreterAreaNever;
@@ -18,6 +22,7 @@
 import org.locationtech.jts.geom.Coordinate;
 import org.locationtech.jts.geom.Envelope;
 import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.MultiPolygon;
 import org.locationtech.jts.geom.Point;
 import org.locationtech.jts.geom.Polygon;
 import org.locationtech.jts.io.ParseException;
@@ -29,8 +34,11 @@
 class OSHDBGeometryBuilderTest extends OSHDBGeometryTest {
   private static final double DELTA = 1E-6;
 
-  OSHDBGeometryBuilderTest() {
-    super("./src/test/resources/geometryBuilder.osh");
+  public OSHDBGeometryBuilderTest() {
+    super(
+        "./src/test/resources/geometryBuilder.osh",
+        "./src/test/resources/relations/multipolygonShellsShareNode.osm"
+    );
   }
 
   @Test
@@ -209,4 +217,111 @@ void testBoundingBoxGetGeometry() {
       new Coordinate(0, 0)};
     assertArrayEquals(test, geometry.getCoordinates());
   }
+
+  /** Test building of multipolygons with self-touching outer rings.
+   *
+   * <p>Example:
+   * <pre>
+   * lat
+   *  ^
+   *  |
+   *  4 a--d
+   *  3 |  |
+   *  2 b--c--e
+   *  1    |  |
+   *  0    g--f
+   *    0 1 2 3456 -> lon
+   *
+   * A: (a,b,c)\
+   *            one shell (a,b,c,d,a)
+   * B: (c,d,a)/
+   *
+   * C: (c,e,f)\
+   *            one shell (c,e,f,g,c)
+   * D: (f,g,c)/
+   * </pre></p>
+   *
+   * <p>Expected:
+   * <pre>
+   * - MULTIPOLYGON (((0 4, 0 2, 3 2, 3 4, 0 4)), ((3 2, 6 2, 6 0, 3 0, 3 2)))
+   * - MULTIPOLYGON (((3 2, 6 2, 6 0, 3 0, 3 2)), ((0 4, 0 2, 3 2, 3 4, 0 4)))
+   * - or similar
+   * </pre></p>
+   */
+  @Test
+  public void testMultipolygonShellsShareNode() {
+    // single figure-8
+    var members = testData.relations().get(100L).get(0).getMembers();
+    getMultipolygonSharedNodeCheck(geom -> {
+      assertTrue(geom instanceof MultiPolygon);
+      assertEquals(2, geom.getNumGeometries());
+    }).accept(members);
+    // all permutations of relation members should lead to the same result
+    members = testData.relations().get(101L).get(0).getMembers();
+    checkAllMemberPermutations(members.length, members, getMultipolygonSharedNodeCheck(geom -> {
+      assertTrue(geom instanceof MultiPolygon);
+      assertEquals(2, geom.getNumGeometries());
+    }));
+    // double loop
+    members = testData.relations().get(102L).get(0).getMembers();
+    getMultipolygonSharedNodeCheck(geom -> {
+      assertTrue(geom instanceof MultiPolygon);
+      assertEquals(3, geom.getNumGeometries());
+    }).accept(members);
+    // second loop forming whole
+    members = testData.relations().get(103L).get(0).getMembers();
+    getMultipolygonSharedNodeCheck(geom -> {
+      assertTrue(geom instanceof Polygon);
+      assertEquals(1, ((Polygon) geom).getNumInteriorRing());
+    }).accept(members);
+    // multiple holes
+    members = testData.relations().get(104L).get(0).getMembers();
+    getMultipolygonSharedNodeCheck(geom -> {
+      assertTrue(geom instanceof Polygon);
+      assertEquals(2, ((Polygon) geom).getNumInteriorRing());
+    }).accept(members);
+    // hole in hole = additional outer
+    members = testData.relations().get(105L).get(0).getMembers();
+    checkAllMemberPermutations(members.length, members, getMultipolygonSharedNodeCheck(geom -> {
+      assertTrue(geom instanceof MultiPolygon);
+      assertEquals(2, geom.getNumGeometries());
+      assertEquals(1,
+          ((Polygon) geom.getGeometryN(0)).getNumInteriorRing()
+              + ((Polygon) geom.getGeometryN(1)).getNumInteriorRing()
+      );
+    }));
+  }
+
+  private Consumer<OSMMember[]> getMultipolygonSharedNodeCheck(Consumer<Geometry> tester) {
+    var areaDecider = new FakeTagInterpreterAreaMultipolygonAllOuters();
+    var timestamp = TimestampParser.toOSHDBTimestamp("2001-01-01");
+    return relMembers -> {
+      var relation = OSM.relation(1, 1, timestamp.getEpochSecond(), 0, 0, null, relMembers);
+      var geom = OSHDBGeometryBuilder.getGeometry(relation, timestamp, areaDecider);
+      assertTrue(geom.isValid());
+      tester.accept(geom);
+    };
+  }
+
+  private static <T> void checkAllMemberPermutations(int n, T[] elements, Consumer<T[]> consumer) {
+    if (n == 1) {
+      consumer.accept(elements);
+    } else {
+      for (int i = 0; i < n - 1; i++) {
+        checkAllMemberPermutations(n - 1, elements, consumer);
+        if (n % 2 == 0) {
+          swap(elements, i, n - 1);
+        } else {
+          swap(elements, 0, n - 1);
+        }
+      }
+      checkAllMemberPermutations(n - 1, elements, consumer);
+    }
+  }
+
+  private static <T> void swap(T[] elements, int a, int b) {
+    T tmp = elements[a];
+    elements[a] = elements[b];
+    elements[b] = tmp;
+  }
 }

oshdb-util/src/test/java/org/heigit/ohsome/oshdb/util/geometry/OSHDBGeometryTest.java

@@ -21,8 +21,10 @@ public abstract class OSHDBGeometryTest {
   protected final ListMultimap<Long, OSMRelation> relations;
   protected final TagInterpreter areaDecider;
 
-  protected OSHDBGeometryTest(String testdata) {
-    testData.add(testdata);
+  protected OSHDBGeometryTest(String... testdata) {
+    for (var data : testdata) {
+      testData.add(data);
+    }
     nodes = testData.nodes();
     ways = testData.ways();
     relations = testData.relations();