Another try at fixing longitude wraparound for bins (#261)

* Another try at fixing longitude wraparound for bins

* Gonna get it right this time

* Forgot to update the comment

* The filter case was not supposed to reinterpret geometry

* Copy antimeridian-crossing geometries to the other side too

* Getting closer to getting antimeridian-crossing polygons right

* Update changelog and version

CHANGELOG.md

@@ -1,3 +1,7 @@
+# 2.62.1
+
+* More work in progress on binning point features in overzoom
+
 # 2.62.0
 
 * Fix another bad interaction, this time between dropping-as-needed and --limit-tile-feature-count

Makefile

@@ -371,6 +371,14 @@ overzoom-test: tippecanoe-overzoom
 	./tippecanoe-decode tests/pbf/0-0-0-pop-2-0-1.pbf.out 2 0 1 > tests/pbf/0-0-0-pop-2-0-1.pbf.out.json.check
 	cmp tests/pbf/0-0-0-pop-2-0-1.pbf.out.json.check tests/pbf/0-0-0-pop-2-0-1.pbf.out.json
 	rm tests/pbf/0-0-0-pop-2-0-1.pbf.out tests/pbf/0-0-0-pop-2-0-1.pbf.out.json.check
+	./tippecanoe-overzoom -o tests/pbf/0-0-0-pop-1-1-0.pbf.out --assign-to-bins tests/pbf/h3-1-1-0.geojson tests/pbf/0-0-0.pbf 1/1/0 1/1/0
+	./tippecanoe-decode tests/pbf/0-0-0-pop-1-1-0.pbf.out 1 1 0 > tests/pbf/0-0-0-pop-1-1-0.pbf.out.json.check
+	cmp tests/pbf/0-0-0-pop-1-1-0.pbf.out.json.check tests/pbf/0-0-0-pop-1-1-0.pbf.out.json
+	rm tests/pbf/0-0-0-pop-1-1-0.pbf.out tests/pbf/0-0-0-pop-1-1-0.pbf.out.json.check
+	./tippecanoe-overzoom -o tests/pbf/0-0-0-pop-0-0-0.pbf.out --assign-to-bins tests/pbf/h3-0-0-0.geojson tests/pbf/0-0-0.pbf 0/0/0 0/0/0
+	./tippecanoe-decode tests/pbf/0-0-0-pop-0-0-0.pbf.out 0 0 0 > tests/pbf/0-0-0-pop-0-0-0.pbf.out.json.check
+	cmp tests/pbf/0-0-0-pop-0-0-0.pbf.out.json.check tests/pbf/0-0-0-pop-0-0-0.pbf.out.json
+	rm tests/pbf/0-0-0-pop-0-0-0.pbf.out tests/pbf/0-0-0-pop-0-0-0.pbf.out.json.check
 
 join-test: tippecanoe tippecanoe-decode tile-join
 	./tippecanoe -q -f -z12 -o tests/join-population/tabblock_06001420.mbtiles -YALAND10:'Land area' -L'{"file": "tests/join-population/tabblock_06001420.json", "description": "population"}'

overzoom.cpp

@@ -197,7 +197,7 @@ int main(int argc, char **argv) {
 			exit(EXIT_OPEN);
 		}
 
-		bins = parse_layers(f, nz, nx, ny, 1LL << detail);
+		bins = parse_layers(f, nz, nx, ny, 1LL << detail, true);
 		fclose(f);
 	}
 

plugin.cpp

@@ -82,7 +82,7 @@ std::vector<mvt_layer> parse_layers(int fd, int z, unsigned x, unsigned y, std::
 		exit(EXIT_OPEN);
 	}
 
-	std::vector<mvt_layer> out = parse_layers(f, z, x, y, extent);
+	std::vector<mvt_layer> out = parse_layers(f, z, x, y, extent, false);
 
 	if (fclose(f) != 0) {
 		perror("fclose postfilter output");

read_json.cpp

@@ -178,7 +178,7 @@ static std::vector<mvt_geometry> to_feature(drawvec &geom) {
 	return out;
 }
 
-std::vector<mvt_layer> parse_layers(FILE *fp, int z, unsigned x, unsigned y, int extent) {
+std::vector<mvt_layer> parse_layers(FILE *fp, int z, unsigned x, unsigned y, int extent, bool fix_longitudes) {
 	std::map<std::string, mvt_layer> ret;
 	std::shared_ptr<std::string> tile_stringpool = std::make_shared<std::string>();
 
@@ -278,6 +278,57 @@ std::vector<mvt_layer> parse_layers(FILE *fp, int z, unsigned x, unsigned y, int
 
 		drawvec dv;
 		parse_geometry(t, coordinates, dv, VT_MOVETO, "Filter output", jp->line, j);
+
+		// handle longitude wraparound
+		//
+		// this is supposed to be data for a single tile,
+		// so any jump from the left hand side edge of the world
+		// to the right edge, or vice versa, is unexpected,
+		// so move it to the other side.
+
+		if (fix_longitudes && mb_geometry[t] == VT_POLYGON) {
+			const long long quarter_world = 1LL << 30;
+			const long long world = 1LL << 32;
+
+			bool copy_to_left = false;
+			bool copy_to_right = false;
+
+			for (size_t i = 0; i < dv.size(); i++) {
+				// is this vertex on a different side of the world
+				// than the first vertex? then shift this one to match
+				if (i > 0) {
+					if ((dv[0].x < quarter_world) && (dv[i].x > 3 * quarter_world)) {
+						dv[i].x -= world;
+					}
+					if ((dv[0].x > 3 * quarter_world) && (dv[i].x < quarter_world)) {
+						dv[i].x += world;
+					}
+				}
+
+				// does it stick off the edge of the world?
+				// then we need another copy on the other side of the world
+				if (dv[i].x < 0) {
+					copy_to_right = true;
+				}
+				if (dv[i].x > world) {
+					copy_to_left = true;
+				}
+			}
+
+			if (copy_to_left) {
+				size_t n = dv.size();
+				for (size_t i = 0; i < n; i++) {
+					dv.emplace_back(dv[i].op, dv[i].x - world, (long long) dv[i].y);
+				}
+			}
+			if (copy_to_right) {
+				size_t n = dv.size();
+				for (size_t i = 0; i < n; i++) {
+					dv.emplace_back(dv[i].op, dv[i].x + world, (long long) dv[i].y);
+				}
+			}
+		}
+
 		if (mb_geometry[t] == VT_POLYGON) {
 			dv = fix_polygon(dv, false, false);
 		}
@@ -286,26 +337,19 @@ std::vector<mvt_layer> parse_layers(FILE *fp, int z, unsigned x, unsigned y, int
 		for (size_t i = 0; i < dv.size(); i++) {
 			long long scale = 1LL << (32 - z);
 
-			// offset
+			// offset to tile
 			dv[i].x -= scale * x;
 			dv[i].y -= scale * y;
 
-			// handle longitude wraparound
+			// scale to tile
-			if (dv[i].x > 2 * scale && dv[i].x - (1LL << 32) > -3 * scale) {
-				dv[i].x -= 1LL << 32;
-			}
-			if (dv[i].x < -3 * scale && dv[i].x + (1LL << 32) < 2 * scale) {
-				dv[i].x += 1LL << 32;
-			}
-
-			// scale
 			dv[i].x = std::round(dv[i].x * extent / (double) scale);
 			dv[i].y = std::round(dv[i].y * extent / (double) scale);
 		}
 
 		if (mb_geometry[t] == VT_POLYGON) {
-			// we can try scaling up because these are tile coordinates
+			// don't try scaling up because we may have coordinates
-			dv = clean_or_clip_poly(dv, 0, 0, false, true);
+			// on the other side of the world
+			dv = clean_or_clip_poly(dv, z, 256, true, false);
 			if (dv.size() < 3) {
 				dv.clear();
 			}

read_json.hpp

@@ -12,6 +12,6 @@ extern int mb_geometry[GEOM_TYPES];
 
 void json_context(json_object *j);
 void parse_geometry(int t, json_object *j, drawvec &out, int op, const char *fname, int line, json_object *feature);
-std::vector<mvt_layer> parse_layers(FILE *fp, int z, unsigned x, unsigned y, int extent);
+std::vector<mvt_layer> parse_layers(FILE *fp, int z, unsigned x, unsigned y, int extent, bool fix_longitudes);
 
 serial_val stringify_value(json_object *value, const char *reading, int line, json_object *feature);

version.hpp

@@ -1,6 +1,6 @@
 #ifndef VERSION_HPP
 #define VERSION_HPP
 
-#define VERSION "v2.62.0"
+#define VERSION "v2.62.1"
 
 #endif