Detect certain kinds of invalid polygons when slicing into tiles (#499)

planetiler-core/src/main/java/com/onthegomap/planetiler/geo/TileCoord.java

@@ -167,4 +167,8 @@ public record TileCoord(int encoded, int x, int y, int z) implements Comparable<
     int dim = 1 << z;
     return x * dim + (dim - 1 - y);
   }
+
+  public TileCoord parent() {
+    return ofXYZ(x / 2, y / 2, z - 1);
+  }
 }

planetiler-core/src/main/java/com/onthegomap/planetiler/geo/TileExtents.java

@@ -52,7 +52,12 @@ public class TileExtents implements Predicate<TileCoord> {
 
       if (mercator != null) {
         Geometry scaled = AffineTransformation.scaleInstance(1 << zoom, 1 << zoom).transform(mercator);
-        TiledGeometry.CoveredTiles covered = TiledGeometry.getCoveredTiles(scaled, zoom, forZoom);
+        TiledGeometry.CoveredTiles covered;
+        try {
+          covered = TiledGeometry.getCoveredTiles(scaled, zoom, forZoom);
+        } catch (GeometryException e) {
+          throw new IllegalArgumentException("Invalid geometry: " + scaled);
+        }
         forZoom = forZoom.withShape(covered);
         LOGGER.info("prepareShapeForZoom z{} {}", zoom, covered);
       }

planetiler-core/src/main/java/com/onthegomap/planetiler/render/FeatureRenderer.java

@@ -189,15 +189,28 @@ public class FeatureRenderer implements Consumer<FeatureCollector.Feature>, Clos
         continue;
       }
 
-      // TODO potential optimization: iteratively simplify z+1 to get z instead of starting with original geom each time
-      // simplify only takes 4-5 minutes of wall time when generating the planet though, so not a big deal
-      Geometry geom = AffineTransformation.scaleInstance(scale, scale).transform(input);
-      geom = DouglasPeuckerSimplifier.simplify(geom, tolerance);
-
-      List<List<CoordinateSequence>> groups = GeometryCoordinateSequences.extractGroups(geom, minSize);
       double buffer = feature.getBufferPixelsAtZoom(z) / 256;
       TileExtents.ForZoom extents = config.bounds().tileExtents().getForZoom(z);
-      TiledGeometry sliced = TiledGeometry.sliceIntoTiles(groups, buffer, area, z, extents);
+
+      // TODO potential optimization: iteratively simplify z+1 to get z instead of starting with original geom each time
+      // simplify only takes 4-5 minutes of wall time when generating the planet though, so not a big deal
+      Geometry scaled = AffineTransformation.scaleInstance(scale, scale).transform(input);
+      TiledGeometry sliced;
+      Geometry geom = DouglasPeuckerSimplifier.simplify(scaled, tolerance);
+      List<List<CoordinateSequence>> groups = GeometryCoordinateSequences.extractGroups(geom, minSize);
+      try {
+        sliced = TiledGeometry.sliceIntoTiles(groups, buffer, area, z, extents);
+      } catch (GeometryException e) {
+        try {
+          geom = GeoUtils.fixPolygon(geom);
+          groups = GeometryCoordinateSequences.extractGroups(geom, minSize);
+          sliced = TiledGeometry.sliceIntoTiles(groups, buffer, area, z, extents);
+        } catch (GeometryException ex) {
+          ex.log(stats, "slice_line_or_polygon", "Error slicing feature at z" + z + ": " + feature);
+          // omit from this zoom level, but maybe the next will be better
+          continue;
+        }
+      }
       Map<String, Object> attrs = feature.getAttrsAtZoom(sliced.zoomLevel());
       if (numPointsAttr != null) {
         // if profile wants the original number of points that the simplified but untiled geometry started with

planetiler-core/src/main/java/com/onthegomap/planetiler/render/TiledGeometry.java

@@ -23,6 +23,7 @@ import com.carrotsearch.hppc.cursors.IntObjectCursor;
 import com.onthegomap.planetiler.collection.Hppc;
 import com.onthegomap.planetiler.collection.IntRangeSet;
 import com.onthegomap.planetiler.geo.GeoUtils;
+import com.onthegomap.planetiler.geo.GeometryException;
 import com.onthegomap.planetiler.geo.MutableCoordinateSequence;
 import com.onthegomap.planetiler.geo.TileCoord;
 import com.onthegomap.planetiler.geo.TileExtents;
@@ -130,7 +131,7 @@ public class TiledGeometry {
    * @return each tile this feature touches, and the points that appear on each
    */
   public static TiledGeometry sliceIntoTiles(Geometry scaledGeom, double minSize, double buffer, int z,
-    TileExtents.ForZoom extents) {
+    TileExtents.ForZoom extents) throws GeometryException {
 
     if (scaledGeom.isEmpty()) {
       // ignore
@@ -161,7 +162,8 @@ public class TiledGeometry {
    * @param extents    The tile extents for this zoom level.
    * @return A {@link CoveredTiles} instance for the tiles that are covered by this geometry.
    */
-  public static CoveredTiles getCoveredTiles(Geometry scaledGeom, int zoom, TileExtents.ForZoom extents) {
+  public static CoveredTiles getCoveredTiles(Geometry scaledGeom, int zoom, TileExtents.ForZoom extents)
+    throws GeometryException {
     if (scaledGeom.isEmpty()) {
       return new CoveredTiles(new RoaringBitmap(), zoom);
     } else if (scaledGeom instanceof Puntal || scaledGeom instanceof Polygonal || scaledGeom instanceof Lineal) {
@@ -190,9 +192,10 @@ public class TiledGeometry {
    * @param z       zoom level
    * @param extents range of tile coordinates within the bounds of the map to generate
    * @return each tile this feature touches, and the points that appear on each
+   * @throws GeometryException for a polygon that is invalid in a way that interferes with clipping
    */
   static TiledGeometry sliceIntoTiles(List<List<CoordinateSequence>> groups, double buffer, boolean area, int z,
-    TileExtents.ForZoom extents) {
+    TileExtents.ForZoom extents) throws GeometryException {
     TiledGeometry result = new TiledGeometry(extents, buffer, z, area);
     EnumSet<Direction> wrapResult = result.sliceWorldCopy(groups, 0);
     if (wrapResult.contains(Direction.RIGHT)) {
@@ -323,8 +326,10 @@ public class TiledGeometry {
    *                content that wraps too far west)
    * @return {@link Direction#LEFT} if there is more content to the west and {@link Direction#RIGHT} if there is more
    *         content to the east.
+   * @throws GeometryException for a polygon that is invalid in a way that interferes with clipping
    */
-  private EnumSet<Direction> sliceWorldCopy(List<List<CoordinateSequence>> groups, int xOffset) {
+  private EnumSet<Direction> sliceWorldCopy(List<List<CoordinateSequence>> groups, int xOffset)
+    throws GeometryException {
     EnumSet<Direction> overflow = EnumSet.noneOf(Direction.class);
     for (List<CoordinateSequence> group : groups) {
       Map<TileCoord, List<CoordinateSequence>> inProgressShapes = new HashMap<>();
@@ -494,7 +499,7 @@ public class TiledGeometry {
    * polygon.
    */
   private IntRangeSet sliceY(CoordinateSequence stripeSegment, int x, boolean outer,
-    Map<TileCoord, List<CoordinateSequence>> inProgressShapes) {
+    Map<TileCoord, List<CoordinateSequence>> inProgressShapes) throws GeometryException {
     if (stripeSegment.size() == 0) {
       return null;
     }
@@ -583,6 +588,13 @@ public class TiledGeometry {
 
           // if this is tile is inside a fill from an outer tile, infer that fill here
           if (area && !outer && toAddTo.isEmpty()) {
+            // since we process outer shells before holes, if a hole is the first thing to intersect
+            // a tile then it must be inside a filled tile from the outer shell. If that's not the case
+            // then the geometry is invalid, so throw an exception so the caller can decide how to handle,
+            // for example fix the polygon then try again.
+            if (!isFilled(x, y)) {
+              throw new GeometryException("bad_polygon_fill", x + ", " + y + " is not filled!");
+            }
             toAddTo.add(fill(buffer));
           }
           toAddTo.add(slice);
@@ -692,6 +704,17 @@ public class TiledGeometry {
     }
   }
 
+  private boolean isFilled(int x, int y) {
+    if (filledRanges == null) {
+      return false;
+    }
+    var filledCol = filledRanges.get(x);
+    if (filledCol == null) {
+      return false;
+    }
+    return filledCol.contains(y);
+  }
+
   private enum Direction {
     RIGHT,
     LEFT

planetiler-core/src/test/java/com/onthegomap/planetiler/PlanetilerTests.java

@@ -1597,6 +1597,24 @@ class PlanetilerTests {
   void testBadRelation() throws Exception {
     // this threw an exception in OsmMultipolygon.build
     OsmXml osmInfo = TestUtils.readOsmXml("bad_spain_relation.xml");
+    List<OsmElement> elements = convertToOsmElements(osmInfo);
+
+    var results = runWithOsmElements(
+      Map.of("threads", "1"),
+      elements,
+      (in, features) -> {
+        if (in.hasTag("landuse", "forest")) {
+          features.polygon("layer")
+            .setZoomRange(12, 14)
+            .setBufferPixels(4);
+        }
+      }
+    );
+
+    assertEquals(11, results.tiles.size());
+  }
+
+  private static List<OsmElement> convertToOsmElements(OsmXml osmInfo) {
     List<OsmElement> elements = new ArrayList<>();
     for (var node : orEmpty(osmInfo.nodes())) {
       elements.add(new OsmElement.Node(node.id(), node.lat(), node.lon()));
@@ -1626,20 +1644,42 @@ class PlanetilerTests {
         }, member.ref(), member.role()));
       }
     }
+    return elements;
+  }
+
+  @Test
+  void testIssue496BaseballMultipolygon() throws Exception {
+    // this generated a polygon that covered an entire z11 tile where the buffer intersected the baseball field
+    OsmXml osmInfo = TestUtils.readOsmXml("issue_496_baseball_multipolygon.xml");
+    List<OsmElement> elements = convertToOsmElements(osmInfo);
 
     var results = runWithOsmElements(
       Map.of("threads", "1"),
       elements,
       (in, features) -> {
-        if (in.hasTag("landuse", "forest")) {
+        if (in.hasTag("natural", "sand")) {
-          features.polygon("layer")
+          features.polygon("test")
-            .setZoomRange(12, 14)
+            .setBufferPixels(4)
-            .setBufferPixels(4);
+            .setPixelTolerance(0.5)
+            .setMinPixelSize(0.1)
+            .setAttr("id", in.id());
         }
       }
     );
 
-    assertEquals(11, results.tiles.size());
+    double areaAtZ14 = 20;
+
+    for (var entry : results.tiles().entrySet()) {
+      var tile = entry.getKey();
+      for (var feature : entry.getValue()) {
+        var geom = feature.geometry().geom();
+        double area = geom.getArea();
+        double expectedMaxArea = areaAtZ14 / (1 << (14 - tile.z()));
+        assertTrue(area < expectedMaxArea, "tile=" + tile + " area=" + area + " geom=" + geom);
+      }
+    }
+
+    assertEquals(8, results.tiles.size());
   }
 
   @ParameterizedTest

planetiler-core/src/test/java/com/onthegomap/planetiler/render/TiledGeometryTest.java

@@ -2,21 +2,28 @@ package com.onthegomap.planetiler.render;
 
 import static org.junit.jupiter.api.Assertions.assertEquals;
 import static org.junit.jupiter.api.Assertions.assertFalse;
+import static org.junit.jupiter.api.Assertions.assertThrows;
 import static org.junit.jupiter.api.Assertions.assertTrue;
 
 import com.onthegomap.planetiler.TestUtils;
+import com.onthegomap.planetiler.geo.GeometryException;
+import com.onthegomap.planetiler.geo.MutableCoordinateSequence;
 import com.onthegomap.planetiler.geo.TileCoord;
 import com.onthegomap.planetiler.geo.TileExtents;
 import java.util.List;
 import java.util.Set;
 import java.util.stream.Collectors;
 import org.junit.jupiter.api.Test;
+import org.junit.jupiter.params.ParameterizedTest;
+import org.junit.jupiter.params.provider.ValueSource;
+import org.locationtech.jts.geom.CoordinateSequence;
+import org.locationtech.jts.geom.util.AffineTransformation;
 
 class TiledGeometryTest {
   private static final int Z14_TILES = 1 << 14;
 
   @Test
-  void testPoint() {
+  void testPoint() throws GeometryException {
     var tiledGeom = TiledGeometry.getCoveredTiles(TestUtils.newPoint(0.5, 0.5), 14,
       new TileExtents.ForZoom(14, 0, 0, Z14_TILES, Z14_TILES, null));
     assertTrue(tiledGeom.test(0, 0));
@@ -36,7 +43,7 @@ class TiledGeometryTest {
   }
 
   @Test
-  void testMultiPoint() {
+  void testMultiPoint() throws GeometryException {
     var tiledGeom = TiledGeometry.getCoveredTiles(TestUtils.newMultiPoint(
       TestUtils.newPoint(0.5, 0.5),
       TestUtils.newPoint(2.5, 1.5)
@@ -49,7 +56,7 @@ class TiledGeometryTest {
   }
 
   @Test
-  void testLine() {
+  void testLine() throws GeometryException {
     var tiledGeom = TiledGeometry.getCoveredTiles(TestUtils.newLineString(
       0.5, 0.5,
       1.5, 0.5
@@ -61,7 +68,7 @@ class TiledGeometryTest {
   }
 
   @Test
-  void testMultiLine() {
+  void testMultiLine() throws GeometryException {
     var tiledGeom = TiledGeometry.getCoveredTiles(TestUtils.newMultiLineString(
       TestUtils.newLineString(
         0.5, 0.5,
@@ -81,7 +88,7 @@ class TiledGeometryTest {
   }
 
   @Test
-  void testPolygon() {
+  void testPolygon() throws GeometryException {
     var tiledGeom =
       TiledGeometry.getCoveredTiles(TestUtils.newPolygon(
         TestUtils.rectangleCoordList(25.5, 27.5),
@@ -103,7 +110,7 @@ class TiledGeometryTest {
   }
 
   @Test
-  void testMultiPolygon() {
+  void testMultiPolygon() throws GeometryException {
     var tiledGeom = TiledGeometry.getCoveredTiles(TestUtils.newMultiPolygon(
       TestUtils.rectangle(25.5, 26.5),
       TestUtils.rectangle(30.1, 30.9)
@@ -120,14 +127,14 @@ class TiledGeometryTest {
   }
 
   @Test
-  void testEmpty() {
+  void testEmpty() throws GeometryException {
     var tiledGeom = TiledGeometry.getCoveredTiles(TestUtils.newGeometryCollection(), 14,
       new TileExtents.ForZoom(14, 0, 0, Z14_TILES, Z14_TILES, null));
     assertEquals(Set.of(), tiledGeom.stream().collect(Collectors.toSet()));
   }
 
   @Test
-  void testGeometryCollection() {
+  void testGeometryCollection() throws GeometryException {
     var tiledGeom = TiledGeometry.getCoveredTiles(TestUtils.newGeometryCollection(
       TestUtils.rectangle(0.1, 0.9),
       TestUtils.newPoint(1.5, 1.5),
@@ -146,4 +153,35 @@ class TiledGeometryTest {
       TileCoord.ofXYZ(4, 10, 14)
     ), tiledGeom.stream().collect(Collectors.toSet()));
   }
+
+  private void rotate(CoordinateSequence coordinateSequence, double x, double y, int degrees) {
+    var transformation = AffineTransformation.rotationInstance(Math.toRadians(degrees), x, y);
+    for (int i = 0; i < coordinateSequence.size(); i++) {
+      transformation.transform(coordinateSequence, i);
+    }
+  }
+
+  @ParameterizedTest
+  @ValueSource(ints = {0, 90, 180, 270})
+  void testOnlyHoleTouchesOtherCellBottom(int degrees) {
+    // hole falls outside shell, and hole touches neighboring tile
+    // make sure that we don't emit that other tile at all
+    MutableCoordinateSequence outer = new MutableCoordinateSequence();
+    outer.addPoint(1.5, 1.5);
+    outer.addPoint(1.6, 1.5);
+    outer.addPoint(1.5, 1.6);
+    outer.closeRing();
+    MutableCoordinateSequence inner = new MutableCoordinateSequence();
+    inner.addPoint(1.4, 1.8);
+    inner.addPoint(1.6, 1.8);
+    inner.addPoint(1.5, 2);
+    inner.closeRing();
+    rotate(outer, 1.5, 1.5, degrees);
+    rotate(inner, 1.5, 1.5, degrees);
+
+    List<List<CoordinateSequence>> coordinateSequences = List.of(List.of(outer, inner));
+    var extent = new TileExtents.ForZoom(11, 0, 0, 1 << 11, 1 << 11, null);
+    assertThrows(GeometryException.class,
+      () -> TiledGeometry.sliceIntoTiles(coordinateSequences, 0.1, true, 11, extent));
+  }
 }

planetiler-core/src/test/resources/issue_496_baseball_multipolygon.xml

@@ -0,0 +1,98 @@
+<osm version="0.6" generator="CGImap 0.8.8 (2922573 spike-06.openstreetmap.org)"
+  copyright="OpenStreetMap and contributors" attribution="http://www.openstreetmap.org/copyright"
+  license="http://opendatacommons.org/licenses/odbl/1-0/">
+  <node id="6797432136" lat="40.9820785" lon="-74.1108750"/>
+  <node id="6797432137" lat="40.9822061" lon="-74.1107262"/>
+  <node id="6797432138" lat="40.9822191" lon="-74.1107365"/>
+  <node id="6797432139" lat="40.9822597" lon="-74.1106672"/>
+  <node id="6797432140" lat="40.9823093" lon="-74.1106122"/>
+  <node id="6797432141" lat="40.9822972" lon="-74.1105398"/>
+  <node id="6797432142" lat="40.9822719" lon="-74.1104781"/>
+  <node id="6797432143" lat="40.9822273" lon="-74.1104110"/>
+  <node id="6797432144" lat="40.9821798" lon="-74.1103641"/>
+  <node id="6797432145" lat="40.9821129" lon="-74.1103386"/>
+  <node id="6797432146" lat="40.9820400" lon="-74.1103319"/>
+  <node id="6797432147" lat="40.9819823" lon="-74.1103413"/>
+  <node id="6797432148" lat="40.9819287" lon="-74.1103641"/>
+  <node id="6797432149" lat="40.9818750" lon="-74.1104244"/>
+  <node id="6797432150" lat="40.9818295" lon="-74.1104915"/>
+  <node id="6797432151" lat="40.9818001" lon="-74.1105572"/>
+  <node id="6797432152" lat="40.9818851" lon="-74.1106913"/>
+  <node id="6797432153" lat="40.9819013" lon="-74.1107007"/>
+  <node id="6797432154" lat="40.9820147" lon="-74.1108697"/>
+  <node id="6797432155" lat="40.9820238" lon="-74.1108549"/>
+  <node id="6797432156" lat="40.9820372" lon="-74.1108494"/>
+  <node id="6797432157" lat="40.9820491" lon="-74.1108449"/>
+  <node id="6797432158" lat="40.9820661" lon="-74.1108567"/>
+  <node id="6797432159" lat="40.9820648" lon="-74.1108530"/>
+  <node id="6797432160" lat="40.9819074" lon="-74.1106779"/>
+  <node id="6797432161" lat="40.9819165" lon="-74.1106591"/>
+  <node id="6797432162" lat="40.9819204" lon="-74.1106372"/>
+  <node id="6797432163" lat="40.9819186" lon="-74.1106149"/>
+  <node id="6797432164" lat="40.9820421" lon="-74.1104727"/>
+  <node id="6797432165" lat="40.9820562" lon="-74.1104781"/>
+  <node id="6797432166" lat="40.9820693" lon="-74.1104791"/>
+  <node id="6797432167" lat="40.9820823" lon="-74.1104763"/>
+  <node id="6797432168" lat="40.9821961" lon="-74.1106444"/>
+  <node id="6797432169" lat="40.9821899" lon="-74.1106658"/>
+  <node id="6797432170" lat="40.9821990" lon="-74.1107026"/>
+  <node id="9152172374" lat="40.9821914" lon="-74.1106855"/>
+  <node id="9152172375" lat="40.9820262" lon="-74.1108496"/>
+  <node id="9152172376" lat="40.9820406" lon="-74.1108444"/>
+  <node id="9152172377" lat="40.9820521" lon="-74.1108498"/>
+  <node id="9152172378" lat="40.9820573" lon="-74.1108478"/>
+
+  <way id="724791255">
+    <nd ref="6797432136"/>
+    <nd ref="6797432137"/>
+    <nd ref="6797432138"/>
+    <nd ref="6797432139"/>
+    <nd ref="6797432140"/>
+    <nd ref="6797432141"/>
+    <nd ref="6797432142"/>
+    <nd ref="6797432143"/>
+    <nd ref="6797432144"/>
+    <nd ref="6797432145"/>
+    <nd ref="6797432146"/>
+    <nd ref="6797432147"/>
+    <nd ref="6797432148"/>
+    <nd ref="6797432149"/>
+    <nd ref="6797432150"/>
+    <nd ref="6797432151"/>
+    <nd ref="6797432152"/>
+    <nd ref="6797432153"/>
+    <nd ref="6797432154"/>
+    <nd ref="6797432155"/>
+    <nd ref="6797432156"/>
+    <nd ref="9152172377"/>
+    <nd ref="6797432158"/>
+    <nd ref="6797432136"/>
+  </way>
+  <way id="724791256">
+    <nd ref="6797432159"/>
+    <nd ref="9152172378"/>
+    <nd ref="6797432157"/>
+    <nd ref="9152172376"/>
+    <nd ref="9152172375"/>
+    <nd ref="6797432160"/>
+    <nd ref="6797432161"/>
+    <nd ref="6797432162"/>
+    <nd ref="6797432163"/>
+    <nd ref="6797432164"/>
+    <nd ref="6797432165"/>
+    <nd ref="6797432166"/>
+    <nd ref="6797432167"/>
+    <nd ref="6797432168"/>
+    <nd ref="6797432169"/>
+    <nd ref="9152172374"/>
+    <nd ref="6797432170"/>
+    <nd ref="6797432159"/>
+    <tag k="landuse" v="grass"/>
+  </way>
+  <relation id="10041968">
+    <member type="way" ref="724791255" role="outer"/>
+    <member type="way" ref="724791256" role="inner"/>
+    <tag k="natural" v="sand"/>
+    <tag k="type" v="multipolygon"/>
+  </relation>
+</osm>