diff --git a/CHANGELOG.md b/CHANGELOG.md index 1e92b60e..dd583729 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/tiles/src/main/java/com/protomaps/basemap/Basemap.java b/tiles/src/main/java/com/protomaps/basemap/Basemap.java index 5c5f50f6..21aa3dcb 100644 --- a/tiles/src/main/java/com/protomaps/basemap/Basemap.java +++ b/tiles/src/main/java/com/protomaps/basemap/Basemap.java @@ -16,8 +16,10 @@ 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; @@ -25,7 +27,7 @@ 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(); } } diff --git a/tiles/src/main/java/com/protomaps/basemap/postprocess/Clip.java b/tiles/src/main/java/com/protomaps/basemap/postprocess/Clip.java new file mode 100644 index 00000000..ea80daa2 --- /dev/null +++ b/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>>> boundaryTilesByZoom; + private final Map 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 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> 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 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> postProcessTile(TileCoord tile, + Map> 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> coords = boundaryTilesByZoom.get(tile.z()).get(tile); + var clippingPoly = reassemblePolygons(coords); + clippingPoly = GeoUtils.fixPolygon(clippingPoly); + clippingPoly.reverse(); + Map> output = new HashMap<>(); + + for (Map.Entry> layer : layers.entrySet()) { + List 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; + } +} diff --git a/tiles/src/test/java/com/protomaps/basemap/layers/LayerTest.java b/tiles/src/test/java/com/protomaps/basemap/layers/LayerTest.java index 9e2ec47d..2caf9cf7 100644 --- a/tiles/src/test/java/com/protomaps/basemap/layers/LayerTest.java +++ b/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> expected, Iterable actual) { var expectedList = expected.stream().toList(); diff --git a/tiles/src/test/java/com/protomaps/basemap/postprocess/ClipTest.java b/tiles/src/test/java/com/protomaps/basemap/postprocess/ClipTest.java new file mode 100644 index 00000000..12995a08 --- /dev/null +++ b/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 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 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 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 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 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 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 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()); + } +} diff --git a/tiles/src/test/resources/clip.geojson b/tiles/src/test/resources/clip.geojson new file mode 100644 index 00000000..015e181a --- /dev/null +++ b/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]]]}} diff --git a/tiles/src/test/resources/empty.geojson b/tiles/src/test/resources/empty.geojson new file mode 100644 index 00000000..0967ef42 --- /dev/null +++ b/tiles/src/test/resources/empty.geojson @@ -0,0 +1 @@ +{}