Add an option to extend zooms if still dropping, but with a limit (#131)

* Add an option to extend zooms if still dropping, but with a limit

* At least when to overzoom, even if not actually doing it yet

* Refactor to give tile-join access to overzoom()

* Didn't work, but *might* have worked

* OK, it did something now

* Ah, there's the bug!

* Hook up pmtiles and dirtiles as overzooming sources

* Add command line option to enable or disable overzooming

* Add (currently broken) test of overzooming in tile-join

* Slightly more abstraction for the tile-join readers

* Factor out duplicated code

* Move construction into a constructor

* More changing accessors to methods

* Reduce magic

* Start tracking a list of the tiles at maxzoom

* I think it worked?

* Add missing #include

* Fix sequence of overzoomed tiles (Y sorts backwards for TMS)

* Don't spend memory on overzooming when we aren't going to use it

* Diff rather than cmp, in the hope of figuring out this broken test

* Keep full coordinate precision if we might extend zooms

* Try a slightly different byte limit

* Make drop-densest more consistent across tile boundaries

* Also affects this test

* Does it behave any differently if it can extend forever?

* I think the discrepancy is a thread-safety problem here

* Revert "Does it behave any differently if it can extend forever?"

This reverts commit 0dff0a0.

* Lost this change to the test

* This time for sure!

* Revert "Also affects this test"

This reverts commit cd1f7c2.

* Revert "Make drop-densest more consistent across tile boundaries"

This reverts commit 563f7d2.

* Revert "Try a slightly different byte limit"

This reverts commit 2e27121.

* Add some more explanatory comments

* Amend the join-test to detect my current bug

* Allow overzooming to complete the zoom if it ever starts

* Forgot to correct the test

* Update changelog and version

* Cleanups from code review

* Remove version number from fixture to fix test

CHANGELOG.md

@@ -1,3 +1,8 @@
+# 2.30.0
+
+* Add --extend-zooms-if-still-dropping-maximum option
+* Add --overzoom option to tile-join
+
 # 2.29.0
 
 * Add tippecanoe-overzoom tool

Makefile

@@ -67,7 +67,7 @@ tippecanoe-enumerate: enumerate.o
 tippecanoe-decode: decode.o projection.o mvt.o write_json.o text.o jsonpull/jsonpull.o dirtiles.o pmtiles_file.o
 	$(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3
 
-tile-join: tile-join.o projection.o mbtiles.o mvt.o memfile.o dirtiles.o jsonpull/jsonpull.o text.o evaluator.o csv.o write_json.o pmtiles_file.o
+tile-join: tile-join.o projection.o mbtiles.o mvt.o memfile.o dirtiles.o jsonpull/jsonpull.o text.o evaluator.o csv.o write_json.o pmtiles_file.o clip.o
 	$(CXX) $(PG) $(LIBS) $(FINAL_FLAGS) $(CXXFLAGS) -o $@ $^ $(LDFLAGS) -lm -lz -lsqlite3 -lpthread
 
 tippecanoe-json-tool: jsontool.o jsonpull/jsonpull.o csv.o text.o geojson-loop.o
@@ -404,8 +404,16 @@ join-test: tile-join
 	./tippecanoe-decode -x generator -x generator_options -x name -x description tests/join-population/empty.dirtiles > tests/join-population/empty.out.json.check
 	cmp tests/join-population/empty.out.json.check tests/join-population/empty.out.json
 	rm -rf tests/join-population/empty.dirtiles tests/join-population/empty.out.dirtiles tests/join-population/empty.out.json.check
-
-
+	#
+	# Test overzooming of tilesets with different maxzooms
+	#
+	mkdir -p tests/ne_110m_ocean/join
+	./tippecanoe -q -z2 -f -o tests/ne_110m_ocean/join/ocean.mbtiles tests/ne_110m_ocean/in.json
+	./tippecanoe -q -z4 -d8 -y name -f -o tests/ne_110m_ocean/join/countries.mbtiles tests/ne_110m_admin_0_countries/in.json.gz
+	./tile-join --overzoom -f -o tests/ne_110m_ocean/join/joined.mbtiles tests/ne_110m_ocean/join/ocean.mbtiles tests/ne_110m_ocean/join/countries.mbtiles
+	./tippecanoe-decode -x generator tests/ne_110m_ocean/join/joined.mbtiles > tests/ne_110m_ocean/join/joined.mbtiles.json.check
+	cmp tests/ne_110m_ocean/join/joined.mbtiles.json.check tests/ne_110m_ocean/join/joined.mbtiles.json
+	rm -f tests/ne_110m_ocean/join/ocean.mbtiles tests/ne_110m_ocean/join/countries.mbtiles tests/ne_110m_ocean/join/joined.mbtiles tests/ne_110m_ocean/join/joined.mbtiles.json.check
 
 join-filter-test:
 	# Comes out different from the direct tippecanoe run because null attributes are lost

README.md

@@ -345,6 +345,8 @@ Parallel processing will also be automatic if the input file is in FlatGeobuf fo
  * `-ae` or `--extend-zooms-if-still-dropping`: Increase the maxzoom if features are still being dropped at that zoom level.
    The detail and simplification options that ordinarily apply only to the maximum zoom level will apply both to the originally
    specified maximum zoom and to any levels added beyond that.
+ * `--extend-zooms-if-still-dropping-maximum=`_count_: Increase the maxzoom if features are still being dropped at that zoom level
+   by up to _count_ zoom levels.
  * `-R` _zoom_`/`_x_`/`_y_ or `--one-tile=`_zoom_`/`_x_`/`_y_: Set the minzoom and maxzoom to _zoom_ and produce only
    the single specified tile at that zoom level.
 
@@ -769,6 +771,11 @@ The options are:
  * `-f` or `--force`: Remove *out.mbtiles* if it already exists.
  * `-r` or `--read-from`: list of input mbtiles to read from.
 
+### Overzooming
+
+ * `--overzoom`: If one of the source tilesets has a larger maxzoom than the others, scale up tiles from the tilesets with the lower maxzooms so they will all have the same maxzoom in the output tileset.
+ * `--buffer=`_pixels_ or `-b` _pixels_: Set the size of the tile buffer in the overzoomed tiles.
+
 ### Tileset description and attribution
 
  * `-A` *attribution* or `--attribution=`*attribution*: Set the attribution string.

clip.cpp

@@ -5,6 +5,8 @@
 #include <mapbox/geometry/snap_rounding.hpp>
 #include "geometry.hpp"
 #include "errors.hpp"
+#include "compression.hpp"
+#include "mvt.hpp"
 
 drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) {
 	drawvec out;
@@ -565,3 +567,127 @@ drawvec close_poly(drawvec &geom) {
 
 	return out;
 }
+
+std::string overzoom(std::string s, int oz, int ox, int oy, int nz, int nx, int ny,
+		     int detail, int buffer, std::set<std::string> const &keep) {
+	mvt_tile tile, outtile;
+	bool was_compressed;
+
+	try {
+		if (!tile.decode(s, was_compressed)) {
+			fprintf(stderr, "Couldn't parse tile %d/%u/%u\n", oz, ox, oy);
+			exit(EXIT_MVT);
+		}
+	} catch (std::exception const &e) {
+		fprintf(stderr, "PBF decoding error in tile %d/%u/%u\n", oz, ox, oy);
+		exit(EXIT_PROTOBUF);
+	}
+
+	for (auto const &layer : tile.layers) {
+		mvt_layer outlayer = mvt_layer();
+
+		int det = detail;
+		if (det <= 0) {
+			det = std::round(log(layer.extent) / log(2));
+		}
+
+		outlayer.name = layer.name;
+		outlayer.version = layer.version;
+		outlayer.extent = 1LL << det;
+
+		for (auto const &feature : layer.features) {
+			mvt_feature outfeature;
+			drawvec geom;
+			int t = feature.type;
+
+			// Convert feature geometry to world coordinates
+
+			long long tilesize = 1LL << (32 - oz);	// source tile size in world coordinates
+			draw ring_closure(0, 0, 0);
+
+			for (auto const &g : feature.geometry) {
+				if (g.op == mvt_closepath) {
+					geom.push_back(ring_closure);
+				} else {
+					geom.emplace_back(g.op,
+							  g.x * tilesize / layer.extent + ox * tilesize,
+							  g.y * tilesize / layer.extent + oy * tilesize);
+
+					if (g.op == mvt_moveto) {
+						ring_closure = geom.back();
+						ring_closure.op = mvt_lineto;
+					}
+				}
+			}
+
+			// Now offset from world coordinates to output tile coordinates,
+			// but retain world scale, because that is what tippecanoe clipping expects
+
+			long long outtilesize = 1LL << (32 - nz);  // destination tile size in world coordinates
+			for (auto &g : geom) {
+				g.x -= nx * outtilesize;
+				g.y -= ny * outtilesize;
+			}
+
+			// Clip to output tile
+
+			if (t == VT_LINE) {
+				geom = clip_lines(geom, nz, buffer);
+			} else if (t == VT_POLYGON) {
+				geom = simple_clip_poly(geom, nz, buffer);
+			} else if (t == VT_POINT) {
+				geom = clip_point(geom, nz, buffer);
+			}
+
+			// Scale to output tile extent
+
+			to_tile_scale(geom, nz, det);
+
+			// Clean geometries
+
+			geom = remove_noop(geom, t, 0);
+			if (t == VT_POLYGON) {
+				geom = clean_or_clip_poly(geom, 0, 0, false);
+				geom = close_poly(geom);
+			}
+
+			// Add geometry to output feature
+
+			outfeature.type = t;
+			for (auto const &g : geom) {
+				outfeature.geometry.emplace_back(g.op, g.x, g.y);
+			}
+
+			// ID and attributes, if it didn't get clipped away
+
+			if (outfeature.geometry.size() > 0) {
+				if (feature.has_id) {
+					outfeature.has_id = true;
+					outfeature.id = feature.id;
+				}
+
+				for (size_t i = 0; i + 1 < feature.tags.size(); i += 2) {
+					if (keep.size() == 0 || keep.find(layer.keys[feature.tags[i]]) != keep.end()) {
+						outlayer.tag(outfeature, layer.keys[feature.tags[i]], layer.values[feature.tags[i + 1]]);
+					}
+				}
+
+				outlayer.features.push_back(outfeature);
+			}
+		}
+
+		if (outlayer.features.size() > 0) {
+			outtile.layers.push_back(outlayer);
+		}
+	}
+
+	if (outtile.layers.size() > 0) {
+		std::string pbf = outtile.encode();
+		std::string compressed;
+		compress(pbf, compressed, true);
+
+		return compressed;
+	} else {
+		return "";
+	}
+}

geometry.hpp

@@ -3,6 +3,8 @@
 
 #include <vector>
 #include <atomic>
+#include <set>
+#include <string>
 #include <sqlite3.h>
 #include <stdio.h>
 
@@ -88,4 +90,7 @@ drawvec clip_point(drawvec &geom, long long x1, long long y1, long long x2, long
 void visvalingam(drawvec &ls, size_t start, size_t end, double threshold, size_t retain);
 int pnpoly(const drawvec &vert, size_t start, size_t nvert, long long testx, long long testy);
 
+std::string overzoom(std::string s, int oz, int ox, int oy, int nz, int nx, int ny,
+		     int detail, int buffer, std::set<std::string> const &keep);
+
 #endif

main.cpp

@@ -91,6 +91,7 @@ size_t limit_tile_feature_count_at_maxzoom = 0;
 unsigned int drop_denser = 0;
 std::map<std::string, serial_val> set_attributes;
 unsigned long long preserve_point_density_threshold = 0;
+long long extend_zooms_max = 0;
 
 std::vector<order_field> order_by;
 bool order_reverse;
@@ -2892,6 +2893,7 @@ int main(int argc, char **argv) {
 		{"minimum-zoom", required_argument, 0, 'Z'},
 		{"smallest-maximum-zoom-guess", required_argument, 0, '~'},
 		{"extend-zooms-if-still-dropping", no_argument, &additional[A_EXTEND_ZOOMS], 1},
+		{"extend-zooms-if-still-dropping-maximum", required_argument, 0, '~'},
 		{"one-tile", required_argument, 0, 'R'},
 
 		{"Tile resolution", 0, 0, 0},
@@ -3138,6 +3140,8 @@ int main(int argc, char **argv) {
 				}
 			} else if (strcmp(opt, "preserve-point-density-threshold") == 0) {
 				preserve_point_density_threshold = atoll_require(optarg, "Preserve point density threshold");
+			} else if (strcmp(opt, "extend-zooms-if-still-dropping-maximum") == 0) {
+				extend_zooms_max = atoll_require(optarg, "Maximum number by which to extend zooms");
 			} else {
 				fprintf(stderr, "%s: Unrecognized option --%s\n", argv[0], opt);
 				exit(EXIT_ARGS);
@@ -3585,7 +3589,7 @@ int main(int argc, char **argv) {
 		}
 	}
 
-	if (extra_detail >= 0) {
+	if (extra_detail >= 0 || additional[A_EXTEND_ZOOMS] || extend_zooms_max > 0) {
 		geometry_scale = 0;
 	} else {
 		geometry_scale = 32 - (full_detail + maxzoom);

main.hpp

@@ -61,6 +61,7 @@ extern int tiny_polygon_size;
 extern size_t limit_tile_feature_count;
 extern size_t limit_tile_feature_count_at_maxzoom;
 extern std::map<std::string, serial_val> set_attributes;
+extern long long extend_zooms_max;
 
 struct order_field {
 	std::string name;

man/tippecanoe.1

@@ -403,6 +403,9 @@ Parallel processing will also be automatic if the input file is in FlatGeobuf fo
 The detail and simplification options that ordinarily apply only to the maximum zoom level will apply both to the originally
 specified maximum zoom and to any levels added beyond that.
 .IP \(bu 2
+\fB\fC\-\-extend\-zooms\-if\-still\-dropping\-maximum=\fR\fIcount\fP: Increase the maxzoom if features are still being dropped at that zoom level
+by up to \fIcount\fP zoom levels.
+.IP \(bu 2
 \fB\fC\-R\fR \fIzoom\fP\fB\fC/\fR\fIx\fP\fB\fC/\fR\fIy\fP or \fB\fC\-\-one\-tile=\fR\fIzoom\fP\fB\fC/\fR\fIx\fP\fB\fC/\fR\fIy\fP: Set the minzoom and maxzoom to \fIzoom\fP and produce only
 the single specified tile at that zoom level.
 .RE
@@ -965,6 +968,13 @@ The options are:
 .IP \(bu 2
 \fB\fC\-r\fR or \fB\fC\-\-read\-from\fR: list of input mbtiles to read from.
 .RE
+.SS Overzooming
+.RS
+.IP \(bu 2
+\fB\fC\-\-overzoom\fR: If one of the source tilesets has a larger maxzoom than the others, scale up tiles from the tilesets with the lower maxzooms so they will all have the same maxzoom in the output tileset.
+.IP \(bu 2
+\fB\fC\-\-buffer=\fR\fIpixels\fP or \fB\fC\-b\fR \fIpixels\fP: Set the size of the tile buffer in the overzoomed tiles.
+.RE
 .SS Tileset description and attribution
 .RS
 .IP \(bu 2