Clip tileset to a polygonal boundary

tiles/src/main/java/com/protomaps/basemap/Basemap.java

@@ -16,16 +16,20 @@
 import com.protomaps.basemap.layers.Roads;
 import com.protomaps.basemap.layers.Transit;
 import com.protomaps.basemap.layers.Water;
+import com.protomaps.basemap.postprocess.Clip;
 import com.protomaps.basemap.text.FontRegistry;
+import java.io.IOException;
+import java.nio.file.Files;
 import java.nio.file.Path;
+import java.nio.file.Paths;
 import java.util.HashMap;
 import java.util.List;
 import java.util.Map;
 
 
 public class Basemap extends ForwardingProfile {
 
-  public Basemap(NaturalEarthDb naturalEarthDb, QrankDb qrankDb) {
+  public Basemap(NaturalEarthDb naturalEarthDb, QrankDb qrankDb, Clip clip) {
 
     var admin = new Boundaries();
     registerHandler(admin);
@@ -72,6 +76,10 @@ public Basemap(NaturalEarthDb naturalEarthDb, QrankDb qrankDb) {
     registerSourceHandler("osm", earth::processOsm);
     registerSourceHandler("osm_land", earth::processPreparedOsm);
     registerSourceHandler("ne", earth::processNe);
+
+    if (clip != null) {
+      registerHandler(clip);
+    }
   }
 
   @Override
@@ -132,14 +140,14 @@ static void run(Arguments args) {
     String area = args.getString("area", "geofabrik area to download", "monaco");
 
     var planetiler = Planetiler.create(args)
-      .addNaturalEarthSource("ne", nePath, neUrl)
+            .addNaturalEarthSource("ne", nePath, neUrl)
       .addOsmSource("osm", Path.of("data", "sources", area + ".osm.pbf"), "geofabrik:" + area)
-      .addShapefileSource("osm_water", sourcesDir.resolve("water-polygons-split-3857.zip"),
-        "https://osmdata.openstreetmap.de/download/water-polygons-split-3857.zip")
-      .addShapefileSource("osm_land", sourcesDir.resolve("land-polygons-split-3857.zip"),
-        "https://osmdata.openstreetmap.de/download/land-polygons-split-3857.zip")
-      .addGeoPackageSource("landcover", sourcesDir.resolve("daylight-landcover.gpkg"),
-        "https://r2-public.protomaps.com/datasets/daylight-landcover.gpkg");
+          .addShapefileSource("osm_water", sourcesDir.resolve("water-polygons-split-3857.zip"),
+            "https://osmdata.openstreetmap.de/download/water-polygons-split-3857.zip")
+          .addShapefileSource("osm_land", sourcesDir.resolve("land-polygons-split-3857.zip"),
+            "https://osmdata.openstreetmap.de/download/land-polygons-split-3857.zip")
+          .addGeoPackageSource("landcover", sourcesDir.resolve("daylight-landcover.gpkg"),
+            "https://r2-public.protomaps.com/datasets/daylight-landcover.gpkg");
 
     Path pgfEncodingZip = sourcesDir.resolve("pgf-encoding.zip");
     Downloader.create(planetiler.config()).add("ne", neUrl, nePath)
@@ -155,9 +163,17 @@ static void run(Arguments args) {
     FontRegistry fontRegistry = FontRegistry.getInstance();
     fontRegistry.setZipFilePath(pgfEncodingZip.toString());
 
+
+
+    Clip clip = null;
+    var clipArg = args.getString("clip","File path to GeoJSON Polygon or MultiPolygon geometry to clip tileset.");
+    if (!clipArg.isEmpty()) {
+      clip = Clip.fromGeoJSONFile(args.getStats(), clipArg);
+    }
+
     fontRegistry.loadFontBundle("NotoSansDevanagari-Regular", "1", "Devanagari");
 
-    planetiler.setProfile(new Basemap(naturalEarthDb, qrankDb)).setOutput(Path.of(area + ".pmtiles"))
+    planetiler.setProfile(new Basemap(naturalEarthDb, qrankDb, clip)).setOutput(Path.of(area + ".pmtiles"))
       .run();
   }
 }

tiles/src/main/java/com/protomaps/basemap/GeoJSON.java

@@ -0,0 +1,55 @@
+package com.protomaps.basemap;
+
+import com.fasterxml.jackson.databind.node.ArrayNode;
+import com.onthegomap.planetiler.geo.GeoUtils;
+import org.locationtech.jts.geom.Coordinate;
+import com.fasterxml.jackson.databind.JsonNode;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.LinearRing;
+import org.locationtech.jts.geom.Polygon;
+
+import java.util.ArrayList;
+import java.util.List;
+
+public class GeoJSON {
+  private static Coordinate[] parseCoordinates(ArrayNode coordinateArray) {
+    Coordinate[] coordinates = new Coordinate[coordinateArray.size()];
+    for (int i = 0; i < coordinateArray.size(); i++) {
+      ArrayNode coordinate = (ArrayNode) coordinateArray.get(i);
+      double x = coordinate.get(0).asDouble();
+      double y = coordinate.get(1).asDouble();
+      coordinates[i] = new Coordinate(x, y);
+    }
+    return coordinates;
+  }
+
+  private static Polygon coordsToPolygon(JsonNode coords) {
+    ArrayNode outerRingNode = (ArrayNode) coords.get(0);
+    Coordinate[] outerRingCoordinates = parseCoordinates(outerRingNode);
+    LinearRing outerRing = GeoUtils.JTS_FACTORY.createLinearRing(outerRingCoordinates);
+
+    LinearRing[] innerRings = new LinearRing[coords.size() - 1];
+    for (int j = 1; j < coords.size(); j++) {
+      ArrayNode innerRingNode = (ArrayNode) coords.get(j);
+      Coordinate[] innerRingCoordinates = parseCoordinates(innerRingNode);
+      innerRings[j - 1] = GeoUtils.JTS_FACTORY.createLinearRing(innerRingCoordinates);
+    }
+    return GeoUtils.JTS_FACTORY.createPolygon(outerRing, innerRings);
+  }
+
+  // return a Polygon or MultiPolygon from a GeoJSON geometry object.
+  public static Geometry parseGeometry(JsonNode jsonGeometry) {
+    var coords = jsonGeometry.get("coordinates");
+    if (jsonGeometry.get("type").asText().equals("Polygon")) {
+      return coordsToPolygon(coords);
+    } else if (jsonGeometry.get("type").asText().equals("MultiPolygon")) {
+      List<Polygon> polygons = new ArrayList<>();
+      for (var polygonCoords : coords) {
+        polygons.add(coordsToPolygon(polygonCoords));
+      }
+      return GeoUtils.createMultiPolygon(polygons);
+    } else {
+      throw new IllegalArgumentException();
+    }
+  }
+}

tiles/src/main/java/com/protomaps/basemap/postprocess/Clip.java

@@ -0,0 +1,155 @@
+package com.protomaps.basemap.postprocess;
+
+import static com.onthegomap.planetiler.geo.GeoUtils.WORLD_BOUNDS;
+import static com.onthegomap.planetiler.geo.GeoUtils.latLonToWorldCoords;
+import static com.onthegomap.planetiler.render.TiledGeometry.getCoveredTiles;
+import static com.onthegomap.planetiler.render.TiledGeometry.sliceIntoTiles;
+
+import com.fasterxml.jackson.databind.JsonNode;
+import com.fasterxml.jackson.databind.ObjectMapper;
+import com.onthegomap.planetiler.ForwardingProfile;
+import com.onthegomap.planetiler.VectorTile;
+import com.onthegomap.planetiler.geo.*;
+import com.onthegomap.planetiler.render.TiledGeometry;
+import java.io.IOException;
+import java.nio.file.Files;
+import java.nio.file.Paths;
+import java.util.*;
+
+import com.onthegomap.planetiler.stats.Stats;
+import com.protomaps.basemap.GeoJSON;
+import org.locationtech.jts.geom.*;
+import org.locationtech.jts.geom.util.AffineTransformation;
+import org.locationtech.jts.operation.overlayng.OverlayNG;
+import org.locationtech.jts.operation.overlayng.OverlayNGRobust;
+
+public class Clip implements ForwardingProfile.TilePostProcessor {
+  private final Map<Integer, Map<TileCoord, List<List<CoordinateSequence>>>> tiledGeometries;
+  private final Map<Integer, TiledGeometry.CoveredTiles> coverings;
+  private final Stats stats;
+
+  public Clip(Stats stats, Geometry input) {
+    this.stats = stats;
+    var clipGeometry = latLonToWorldCoords(input).buffer(0.00001);
+    tiledGeometries = new HashMap<>();
+    coverings = new HashMap<>();
+    try {
+      for (var i = 0; i <= 15; i++) {
+        var extents = TileExtents.computeFromWorldBounds(i, WORLD_BOUNDS);
+        double scale = 1 << i;
+        Geometry scaled = AffineTransformation.scaleInstance(scale, scale).transform(clipGeometry);
+//        var simplified = DouglasPeuckerSimplifier.simplify(scaled, 0.25/256);
+        this.tiledGeometries.put(i, sliceIntoTiles(scaled, 0, 0.015625, i, extents.getForZoom(i)).getTileData());
+        this.coverings.put(i, getCoveredTiles(scaled, i, extents.getForZoom(i)));
+      }
+    } catch (GeometryException e) {
+      throw new RuntimeException("Error clipping");
+    }
+  }
+
+  public static Clip fromGeoJSONFile(Stats stats, String filename) {
+    try {
+      return fromGeoJSON(stats, Files.readAllBytes(Paths.get(filename)));
+    } catch (IOException e) {
+      throw new IllegalArgumentException("Could not open clip file");
+    }
+  }
+
+  public static Clip fromGeoJSON(Stats stats, byte[] bytes) {
+    try {
+      ObjectMapper mapper = new ObjectMapper();
+      JsonNode geoJson = mapper.readTree(bytes);
+      return new Clip(stats, GeoJSON.parseGeometry(geoJson));
+    } catch (IOException e) {
+      throw new IllegalArgumentException("Clip GeoJSON is invalid");
+    }
+  }
+
+  // Copied from elsewhere in planetiler
+  private static Polygon reassemblePolygon(List<CoordinateSequence> group) throws GeometryException {
+    try {
+      LinearRing first = GeoUtils.JTS_FACTORY.createLinearRing(group.getFirst());
+      LinearRing[] rest = new LinearRing[group.size() - 1];
+      for (int j = 1; j < group.size(); j++) {
+        CoordinateSequence seq = group.get(j);
+        CoordinateSequences.reverse(seq);
+        rest[j - 1] = GeoUtils.JTS_FACTORY.createLinearRing(seq);
+      }
+      return GeoUtils.JTS_FACTORY.createPolygon(first, rest);
+    } catch (IllegalArgumentException e) {
+      throw new GeometryException("reassemble_polygon_failed", "Could not build polygon", e);
+    }
+  }
+
+  // Copied from elsewhere in Planetiler
+  static Geometry reassemblePolygons(List<List<CoordinateSequence>> groups) throws GeometryException {
+    int numGeoms = groups.size();
+    if (numGeoms == 1) {
+      return reassemblePolygon(groups.getFirst());
+    } else {
+      Polygon[] polygons = new Polygon[numGeoms];
+      for (int i = 0; i < numGeoms; i++) {
+        polygons[i] = reassemblePolygon(groups.get(i));
+      }
+      return GeoUtils.JTS_FACTORY.createMultiPolygon(polygons);
+    }
+  }
+
+  @Override
+  public Map<String, List<VectorTile.Feature>> postProcessTile(TileCoord tileCoord,
+    Map<String, List<VectorTile.Feature>> map) throws GeometryException {
+    if (this.coverings.containsKey(tileCoord.z()) &&
+      this.coverings.get(tileCoord.z()).test(tileCoord.x(), tileCoord.y())) {
+      if (this.tiledGeometries.containsKey(tileCoord.z()) && this.tiledGeometries.get(tileCoord.z()).containsKey(tileCoord)) {
+        List<List<CoordinateSequence>> coords = tiledGeometries.get(tileCoord.z()).get(tileCoord);
+        var clipGeometry = reassemblePolygons(coords);
+        var clipGeometry2 = GeoUtils.fixPolygon(clipGeometry);
+        clipGeometry2.reverse();
+        Map<String, List<VectorTile.Feature>> output = new HashMap<>();
+
+        for (Map.Entry<String, List<VectorTile.Feature>> layer : map.entrySet()) {
+          List<VectorTile.Feature> clippedFeatures = new ArrayList<>();
+          for (var feature : layer.getValue()) {
+            try {
+              var newGeom = OverlayNGRobust.overlay(feature.geometry().decode(), clipGeometry2, OverlayNG.INTERSECTION);
+              if (!newGeom.isEmpty() && newGeom.getNumGeometries() > 0) {
+                if (newGeom instanceof Polygonal) {
+                  newGeom = GeoUtils.snapAndFixPolygon(newGeom, stats, "clip");
+                  newGeom = newGeom.reverse();
+                  if (!newGeom.isEmpty() && newGeom.getNumGeometries() > 0) {
+                    if (newGeom instanceof GeometryCollection) {
+                      for (int i = 0; i < newGeom.getNumGeometries(); i++) {
+                        // geometrycollection
+                        clippedFeatures.add(feature.copyWithNewGeometry(newGeom.getGeometryN(i)));
+                      }
+                    } else {
+                      // a multipolygon/polygon
+                      clippedFeatures.add(feature.copyWithNewGeometry(newGeom));
+                    }
+                  }
+                } else {
+                  if (!newGeom.isEmpty() && newGeom.getNumGeometries() > 0) {
+                    if (newGeom instanceof GeometryCollection) {
+                      for (int i = 0; i < newGeom.getNumGeometries(); i++) {
+                        clippedFeatures.add(feature.copyWithNewGeometry(newGeom.getGeometryN(i)));
+                      }
+                    } else {
+                      clippedFeatures.add(feature.copyWithNewGeometry(newGeom));
+                    }
+                  }
+                }
+              }
+            } catch (GeometryException e) {
+              System.err.println("Could not clip geometry");
+            }
+          }
+
+          output.put(layer.getKey(), clippedFeatures);
+        }
+        return output;
+      }
+      return map;
+    }
+    return Map.of();
+  }
+}

tiles/src/test/java/com/protomaps/basemap/GeoJSONTest.java

@@ -0,0 +1,4 @@
+package com.protomaps.basemap;
+
+public class GeoJSONTest {
+}

tiles/src/test/java/com/protomaps/basemap/layers/LayerTest.java

@@ -24,7 +24,7 @@ abstract class LayerTest {
     List.of(new NaturalEarthDb.NeAdmin1StateProvince("California", "US-CA", "Q2", 5.0, 8.0)),
     List.of(new NaturalEarthDb.NePopulatedPlace("San Francisco", "Q3", 9.0, 2))
   );
-  final Basemap profile = new Basemap(naturalEarthDb, null);
+  final Basemap profile = new Basemap(naturalEarthDb, null, null);
 
   static void assertFeatures(int zoom, List<Map<String, Object>> expected, Iterable<FeatureCollector.Feature> actual) {
     var expectedList = expected.stream().toList();

tiles/src/test/java/com/protomaps/basemap/postprocess/ClipTest.java

@@ -0,0 +1,4 @@
+package com.protomaps.basemap.postprocess;
+
+public class ClipTest {
+}