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 static int tmsXYToPosition(int z, int x, int y) {
     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 static TileExtents computeFromWorldBounds(int maxzoom, Envelope worldBoun
 
       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 @@ private void renderLineOrPolygon(FeatureCollector.Feature feature, Geometry inpu
         continue;
       }
 
+      double buffer = feature.getBufferPixelsAtZoom(z) / 256;
+      TileExtents.ForZoom extents = config.bounds().tileExtents().getForZoom(z);
+
       // 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);
-
+      Geometry scaled = AffineTransformation.scaleInstance(scale, scale).transform(input);
+      TiledGeometry sliced;
+      Geometry geom = DouglasPeuckerSimplifier.simplify(scaled, 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);
+      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.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 @@ private static int wrapInt(int value, int max) {
    * @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 static TiledGeometry sliceIntoTiles(Geometry scaledGeom, double minSize,
    * @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 static CoveredTiles getCoveredTiles(Geometry scaledGeom, int zoom, TileEx
    * @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 Map<TileCoord, List<List<CoordinateSequence>>> getTileData() {
    *                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 @@ private IntObjectMap<List<MutableCoordinateSequence>> sliceX(CoordinateSequence
    * 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 @@ record SkippedSegment(Direction side, int lo, int hi) {}
 
           // 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 @@ private void removeFilledRange(int x, IntRangeSet yRange) {
     }
   }
 
+  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 @@ private static <T> List<T> orEmpty(List<T> in) {
   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 @@ void testBadRelation() throws Exception {
         }, 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")) {
-          features.polygon("layer")
-            .setZoomRange(12, 14)
-            .setBufferPixels(4);
+        if (in.hasTag("natural", "sand")) {
+          features.polygon("test")
+            .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 @@
 
 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 @@ void testPoint() {
   }
 
   @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 @@ void testMultiPoint() {
   }
 
   @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 @@ void testLine() {
   }
 
   @Test
-  void testMultiLine() {
+  void testMultiLine() throws GeometryException {
     var tiledGeom = TiledGeometry.getCoveredTiles(TestUtils.newMultiLineString(
       TestUtils.newLineString(
         0.5, 0.5,
@@ -81,7 +88,7 @@ void testMultiLine() {
   }
 
   @Test
-  void testPolygon() {
+  void testPolygon() throws GeometryException {
     var tiledGeom =
       TiledGeometry.getCoveredTiles(TestUtils.newPolygon(
         TestUtils.rectangleCoordList(25.5, 27.5),
@@ -103,7 +110,7 @@ void testPolygon() {
   }
 
   @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 @@ void testMultiPolygon() {
   }
 
   @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 @@ void testGeometryCollection() {
       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));
+  }
 }