Clip tileset to a polygonal boundary

CHANGELOG.md

@@ -1,3 +1,7 @@
+Tiles v4.2.0
+------
+- add `--clip` option to tile generation which clips the entire tileset by a polygon or multipolygon. [#51]
+
 Styles v4.4.0
 ------
 - Improve boundary appearances at low zooms

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

@@ -16,16 +16,18 @@
 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.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 +74,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
@@ -86,7 +92,7 @@ public String description() {
 
   @Override
   public String version() {
-    return "4.1.0";
+    return "4.2.0";
   }
 
   @Override
@@ -155,9 +161,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(), planetiler.config().minzoom(), planetiler.config().maxzoom(), true,
+          Paths.get(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/postprocess/Clip.java

@@ -0,0 +1,159 @@
+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.onthegomap.planetiler.ForwardingProfile;
+import com.onthegomap.planetiler.Planetiler;
+import com.onthegomap.planetiler.VectorTile;
+import com.onthegomap.planetiler.geo.*;
+import com.onthegomap.planetiler.reader.FileFormatException;
+import com.onthegomap.planetiler.reader.geojson.GeoJson;
+import com.onthegomap.planetiler.render.TiledGeometry;
+import com.onthegomap.planetiler.stats.Stats;
+import java.nio.file.Path;
+import java.util.*;
+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>>>> boundaryTilesByZoom;
+  private final Map<Integer, TiledGeometry.CoveredTiles> coveredTilesByZoom;
+  private final Stats stats;
+
+  static final double DEFAULT_BUFFER = 4.0 / 256.0;
+
+  // A TilePostProcessor that clips all layers to a given geometry.
+  // the geometry must be in world coordinates ( world from 0 to 1 )
+  public Clip(Stats stats, int minzoom, int maxzoom, boolean doBuffer, Geometry input) {
+    this.stats = stats;
+    double bufferAmount = 0;
+    if (doBuffer) {
+      var envelope = input.getEnvelope().getEnvelopeInternal();
+      bufferAmount = Math.max(envelope.getWidth(), envelope.getHeight()) * DEFAULT_BUFFER;
+    }
+    var clipGeometry = input.buffer(bufferAmount);
+    boundaryTilesByZoom = new HashMap<>();
+    coveredTilesByZoom = new HashMap<>();
+    try {
+      for (var i = minzoom; i <= maxzoom; i++) {
+        var extents = TileExtents.computeFromWorldBounds(i, WORLD_BOUNDS);
+        double scale = 1 << i;
+        Geometry scaled = AffineTransformation.scaleInstance(scale, scale).transform(clipGeometry);
+        this.boundaryTilesByZoom.put(i,
+          sliceIntoTiles(scaled, 0, DEFAULT_BUFFER, i, extents.getForZoom(i)).getTileData());
+        this.coveredTilesByZoom.put(i, getCoveredTiles(scaled, i, extents.getForZoom(i)));
+      }
+    } catch (GeometryException e) {
+      throw new Planetiler.PlanetilerException("Error clipping", e);
+    }
+  }
+
+  public static Clip fromGeoJSONFile(Stats stats, int minzoom, int maxzoom, boolean doBuffer, Path path) {
+    var g = GeoJson.from(path);
+    if (g.count() == 0) {
+      throw new FileFormatException("Empty clipping geometry");
+    }
+    var feature = g.iterator().next();
+    return new Clip(stats, minzoom, maxzoom, doBuffer, latLonToWorldCoords(feature.geometry()));
+  }
+
+  // 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);
+    }
+  }
+
+  private boolean nonDegenerateGeometry(Geometry geom) {
+    return !geom.isEmpty() && geom.getNumGeometries() > 0;
+  }
+
+  private Geometry fixGeometry(Geometry geom) throws GeometryException {
+    if (geom instanceof Polygonal) {
+      geom = GeoUtils.snapAndFixPolygon(geom, stats, "clip");
+      return geom.reverse();
+    }
+    return geom;
+  }
+
+  private void addToFeatures(List<VectorTile.Feature> features, VectorTile.Feature feature, Geometry geom) {
+    if (nonDegenerateGeometry(geom)) {
+      if (geom instanceof GeometryCollection) {
+        for (int i = 0; i < geom.getNumGeometries(); i++) {
+          features.add(feature.copyWithNewGeometry(geom.getGeometryN(i)));
+        }
+      } else {
+        features.add(feature.copyWithNewGeometry(geom));
+      }
+    }
+  }
+
+  @Override
+  public Map<String, List<VectorTile.Feature>> postProcessTile(TileCoord tile,
+    Map<String, List<VectorTile.Feature>> layers) throws GeometryException {
+
+    var inCovering =
+      this.coveredTilesByZoom.containsKey(tile.z()) && this.coveredTilesByZoom.get(tile.z()).test(tile.x(), tile.y());
+
+    if (!inCovering)
+      return Map.of();
+
+    var inBoundary =
+      this.boundaryTilesByZoom.containsKey(tile.z()) && this.boundaryTilesByZoom.get(tile.z()).containsKey(tile);
+
+    if (!inBoundary)
+      return layers;
+
+    List<List<CoordinateSequence>> coords = boundaryTilesByZoom.get(tile.z()).get(tile);
+    var clippingPoly = reassemblePolygons(coords);
+    clippingPoly = GeoUtils.fixPolygon(clippingPoly);
+    clippingPoly.reverse();
+    Map<String, List<VectorTile.Feature>> output = new HashMap<>();
+
+    for (Map.Entry<String, List<VectorTile.Feature>> layer : layers.entrySet()) {
+      List<VectorTile.Feature> clippedFeatures = new ArrayList<>();
+      for (var feature : layer.getValue()) {
+        try {
+          var clippedGeom =
+            OverlayNGRobust.overlay(feature.geometry().decode(), clippingPoly, OverlayNG.INTERSECTION);
+          if (nonDegenerateGeometry(clippedGeom)) {
+            addToFeatures(clippedFeatures, feature, fixGeometry(clippedGeom));
+          }
+        } catch (GeometryException e) {
+          e.log(stats, "clip", "Failed to clip geometry");
+        }
+      }
+      if (!clippedFeatures.isEmpty())
+        output.put(layer.getKey(), clippedFeatures);
+    }
+    return output;
+  }
+}

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,150 @@
+package com.protomaps.basemap.postprocess;
+
+import static com.onthegomap.planetiler.TestUtils.newLineString;
+import static com.onthegomap.planetiler.TestUtils.newPolygon;
+import static org.junit.jupiter.api.Assertions.*;
+
+import com.onthegomap.planetiler.VectorTile;
+import com.onthegomap.planetiler.geo.GeometryException;
+import com.onthegomap.planetiler.geo.TileCoord;
+import com.onthegomap.planetiler.reader.FileFormatException;
+import com.onthegomap.planetiler.stats.Stats;
+import java.nio.file.Path;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Map;
+import org.junit.jupiter.api.Test;
+
+class ClipTest {
+  private final Stats stats = Stats.inMemory();
+
+  @Test
+  void testLoadGeoJSON() {
+    Path cwd = Path.of("").toAbsolutePath();
+    Path pathFromRoot = Path.of("tiles", "src", "test", "resources", "clip.geojson");
+    var clip = Clip.fromGeoJSONFile(stats, 0, 0, false, cwd.resolveSibling(pathFromRoot));
+    assertNotNull(clip);
+  }
+
+  @Test
+  void testLoadNonJSON() {
+    Path cwd = Path.of("").toAbsolutePath();
+    Path pathFromRoot = Path.of("tiles", "src", "test", "resources", "empty.geojson");
+    Path path = cwd.resolveSibling(pathFromRoot);
+    assertThrows(FileFormatException.class, () -> {
+      Clip.fromGeoJSONFile(stats, 0, 0, false, path);
+    });
+  }
+
+  @Test
+  void testClipLine() throws GeometryException {
+    List<VectorTile.Feature> unclipped = new ArrayList<>();
+    unclipped.add(new VectorTile.Feature("layer", 1,
+      // a horizontal line in the across the middle of the 0,0,0 tile.
+      VectorTile.encodeGeometry(newLineString(0, 128, 256, 128)),
+      Map.of("foo", "bar")
+    ));
+
+    // a rectangle that is 50% of the earths width, centered at null island.
+    var n = new Clip(stats, 0, 0, false, newPolygon(0.25, 0.25, 0.75, 0.25, 0.75, 0.75, 0.25, 0.75, 0.25, 0.25));
+    var clipped = n.postProcessTile(TileCoord.ofXYZ(0, 0, 0), Map.of("layer", unclipped));
+
+    assertEquals(1, clipped.size());
+    assertEquals(1, clipped.get("layer").size());
+    assertEquals(newLineString(64, 128, 192, 128), clipped.get("layer").getFirst().geometry().decode());
+  }
+
+  @Test
+  void testClipLineMulti() throws GeometryException {
+    List<VectorTile.Feature> unclipped = new ArrayList<>();
+    unclipped.add(new VectorTile.Feature("layer", 1,
+      // a V shape that enters and leaves the clipping square
+      VectorTile.encodeGeometry(newLineString(32, 128, 128, 224, 224, 128)),
+      Map.of("foo", "bar")
+    ));
+
+    // a rectangle that is 50% of the earths width, centered at null island.
+    var n = new Clip(stats, 0, 0, false, newPolygon(0.25, 0.25, 0.75, 0.25, 0.75, 0.75, 0.25, 0.75, 0.25, 0.25));
+    var clipped = n.postProcessTile(TileCoord.ofXYZ(0, 0, 0), Map.of("layer", unclipped));
+
+    assertEquals(1, clipped.size());
+    assertEquals(2, clipped.get("layer").size());
+    assertEquals(newLineString(64, 160, 96, 192), clipped.get("layer").get(0).geometry().decode());
+    assertEquals(newLineString(160, 192, 192, 160), clipped.get("layer").get(1).geometry().decode());
+  }
+
+  @Test
+  void testClipPolygon() throws GeometryException {
+    List<VectorTile.Feature> unclipped = new ArrayList<>();
+    unclipped.add(new VectorTile.Feature("layer", 1,
+      VectorTile.encodeGeometry(newPolygon(32, 160, 96, 160, 96, 224, 32, 224, 32, 160)),
+      Map.of("foo", "bar")
+    ));
+
+    // a rectangle that is 50% of the earths width, centered at null island.
+    var n = new Clip(stats, 0, 0, false, newPolygon(0.25, 0.25, 0.75, 0.25, 0.75, 0.75, 0.25, 0.75, 0.25, 0.25));
+    var clipped = n.postProcessTile(TileCoord.ofXYZ(0, 0, 0), Map.of("layer", unclipped));
+
+    assertEquals(1, clipped.size());
+    assertEquals(1, clipped.get("layer").size());
+    assertEquals(newPolygon(64, 160, 96, 160, 96, 192, 64, 192, 64, 160),
+      clipped.get("layer").getFirst().geometry().decode());
+  }
+
+  @Test
+  void testClipBelowMinZoom() throws GeometryException {
+    List<VectorTile.Feature> unclipped = new ArrayList<>();
+    unclipped.add(new VectorTile.Feature("layer", 1,
+      VectorTile.encodeGeometry(newLineString(0, 128, 256, 128)),
+      Map.of("foo", "bar")
+    ));
+
+    var n = new Clip(stats, 1, 1, false, newPolygon(0.25, 0.25, 0.75, 0.25, 0.75, 0.75, 0.25, 0.75, 0.25, 0.25));
+    var clipped = n.postProcessTile(TileCoord.ofXYZ(0, 0, 0), Map.of("layer", unclipped));
+    assertEquals(0, clipped.size());
+  }
+
+  @Test
+  void testClipWhollyOutside() throws GeometryException {
+    List<VectorTile.Feature> unclipped = new ArrayList<>();
+    unclipped.add(new VectorTile.Feature("layer", 1,
+      VectorTile.encodeGeometry(newLineString(0, 1, 5, 1)),
+      Map.of("foo", "bar")
+    ));
+
+    var n = new Clip(stats, 0, 0, false, newPolygon(0.25, 0.25, 0.75, 0.25, 0.75, 0.75, 0.25, 0.75, 0.25, 0.25));
+    var clipped = n.postProcessTile(TileCoord.ofXYZ(0, 0, 0), Map.of("layer", unclipped));
+    assertEquals(0, clipped.size());
+  }
+
+  @Test
+  void testClipInInterior() throws GeometryException {
+    List<VectorTile.Feature> unclipped = new ArrayList<>();
+    unclipped.add(new VectorTile.Feature("layer", 1,
+      VectorTile.encodeGeometry(newLineString(0, 1, 5, 1)),
+      Map.of("foo", "bar")
+    ));
+
+    var n = new Clip(stats, 0, 3, false, newPolygon(0.25, 0.25, 0.75, 0.25, 0.75, 0.75, 0.25, 0.75, 0.25, 0.25));
+    var clipped = n.postProcessTile(TileCoord.ofXYZ(3, 3, 3), Map.of("layer", unclipped));
+    assertEquals(1, clipped.size());
+    assertEquals(1, clipped.get("layer").size());
+  }
+
+  @Test
+  void testClipLineBuffer() throws GeometryException {
+    List<VectorTile.Feature> unclipped = new ArrayList<>();
+    unclipped.add(new VectorTile.Feature("layer", 1,
+      VectorTile.encodeGeometry(newLineString(0, 128, 256, 128)),
+      Map.of("foo", "bar")
+    ));
+
+    // a rectangle that is 50% of the earths width, centered at null island.
+    var n = new Clip(stats, 0, 0, true, newPolygon(0.25, 0.25, 0.75, 0.25, 0.75, 0.75, 0.25, 0.75, 0.25, 0.25));
+    var clipped = n.postProcessTile(TileCoord.ofXYZ(0, 0, 0), Map.of("layer", unclipped));
+
+    assertEquals(1, clipped.size());
+    assertEquals(1, clipped.get("layer").size());
+    assertEquals(newLineString(62, 128, 194, 128), clipped.get("layer").getFirst().geometry().decode());
+  }
+}

tiles/src/test/resources/clip.geojson

@@ -0,0 +1 @@
+{"type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":[[[7.417252411949448,43.73567091721708],[7.42905253797835,43.73567091721708],[7.42905253797835,43.7275042719568],[7.417252411949448,43.7275042719568],[7.417252411949448,43.73567091721708]]]}}

tiles/src/test/resources/empty.geojson

@@ -0,0 +1 @@
+{}