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

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
@@ -120,7 +120,7 @@ testargs = \
 %.json.check:
 	./tippecanoe -q -a@ -f -o [email protected] $(call testargs,$(patsubst %.json.check,%,$(word 4,$(subst /, ,$@)))) $(foreach suffix,$(suffixes),$(sort $(wildcard $(subst $(SPACE),/,$(wordlist 1,2,$(subst /, ,$@)))/*.$(suffix)))) < /dev/null
 	./tippecanoe-decode -x generator [email protected] > [email protected]
-	cmp [email protected] $(patsubst %.check,%,$@)
+	diff [email protected] $(patsubst %.check,%,$@)
 	rm [email protected] [email protected]
 
 # Don't test overflow with geobuf, because it fails (https://github.com/mapbox/geobuf/issues/87)
@@ -137,7 +137,7 @@ fewer-tests: tippecanoe tippecanoe-decode geobuf-test raw-tiles-test parallel-te
 	for i in $(wildcard $(subst $(SPACE),/,$(wordlist 1,2,$(subst /, ,$@)))/*.json.gz); do gzip -dc $$i | ./tippecanoe-json-tool -w | ./node_modules/geobuf/bin/json2geobuf > $$i.geobuf; done
 	./tippecanoe -q -a@ -f -o [email protected] $(call testargs,$(patsubst %.json.checkbuf,%,$(word 4,$(subst /, ,$@)))) $(foreach suffix,$(suffixes),$(addsuffix .geobuf,$(sort $(wildcard $(subst $(SPACE),/,$(wordlist 1,2,$(subst /, ,$@)))/*.$(suffix))))) < /dev/null
 	./tippecanoe-decode -x generator [email protected] | sed 's/checkbuf/check/g' | sed 's/\.geobuf//g' > [email protected]
-	cmp [email protected] $(patsubst %.checkbuf,%,$@)
+	diff [email protected] $(patsubst %.checkbuf,%,$@)
 	rm [email protected] [email protected]
 
 parallel-test: $(eval SHELL:=$(ADVSHELL))
@@ -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 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

@@ -341,6 +341,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.
 
@@ -765,6 +767,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\fCextend\-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