Merge pull request #170 from mapbox/sane-polygon

Split complex polygons into multiple features
2026-10-02 08:25:40 +02:00 · 2016-02-11 14:29:33 -08:00
parent a1da515dc9 84a2a1bd9d
commit eecdf7747a
24 changed files with 371 additions and 71 deletions
@@ -1,3 +1,8 @@
+## 1.9.0
+
+* Claim vector tile version 2 in mbtiles
+* Split too-complex polygons into multiple features
+
 ## 1.8.1

 * Bug fixes to maxzoom, and more tests
@@ -114,6 +114,7 @@ Options
 * -pk: Don't limit tiles to 500K bytes
 * -pd: Dynamically drop some fraction of features from large tiles to keep them under the 500K size limit. It will probably look ugly at the tile boundaries.
 * -pi: Preserve the original input order of features as the drawing order instead of ordering geographically. (This is implemented as a restoration of the original order at the end, so that dot-dropping is still geographic, which means it also undoes -ao).
+ * -pp: Don't split complex polygons (over 700 vertices after simplification) into multiple features.
 * -q: Work quietly instead of reporting progress

 Example
@@ -202,6 +203,9 @@ have their probability diffused, so that some of them will be drawn as a square
 this minimum size and others will not be drawn at all, preserving the total area that
 all of them should have had together.

+Any polygons that have over 700 vertices after line simplification will be split into
+multiple features so they can be rendered efficiently, unless you use -pp to prevent this.
+
 Features in the same tile that share the same type and attributes are coalesced
 together into a single geometry. You are strongly encouraged to use -x to exclude
 any unnecessary properties to reduce wasted file size.
@@ -341,21 +341,19 @@ drawvec close_poly(drawvec &geom) {
 	return out;
 }

-static bool inside(draw d, int edge, long long area, long long buffer) {
-	long long clip_buffer = buffer * area / 256;
-
+static bool inside(draw d, int edge, long long minx, long long miny, long long maxx, long long maxy) {
 	switch (edge) {
 	case 0:  // top
-		return d.y > -clip_buffer;
+		return d.y > miny;

 	case 1:  // right
-		return d.x < area + clip_buffer;
+		return d.x < maxx;

 	case 2:  // bottom
-		return d.y < area + clip_buffer;
+		return d.y < maxy;

 	case 3:  // left
-		return d.x > -clip_buffer;
+		return d.x > minx;
 	}

 	fprintf(stderr, "internal error inside\n");
@@ -376,24 +374,22 @@ static draw get_line_intersection(draw p0, draw p1, draw p2, draw p3) {
 	return draw(VT_LINETO, p0.x + (t * s1_x), p0.y + (t * s1_y));
 }

-static draw intersect(draw a, draw b, int edge, long long area, long long buffer) {
-	long long clip_buffer = buffer * area / 256;
-
+static draw intersect(draw a, draw b, int edge, long long minx, long long miny, long long maxx, long long maxy) {
 	switch (edge) {
 	case 0:  // top
-		return get_line_intersection(a, b, draw(VT_MOVETO, -clip_buffer, -clip_buffer), draw(VT_MOVETO, area + clip_buffer, -clip_buffer));
+		return get_line_intersection(a, b, draw(VT_MOVETO, minx, miny), draw(VT_MOVETO, maxx, miny));
 		break;

 	case 1:  // right
-		return get_line_intersection(a, b, draw(VT_MOVETO, area + clip_buffer, -clip_buffer), draw(VT_MOVETO, area + clip_buffer, area + clip_buffer));
+		return get_line_intersection(a, b, draw(VT_MOVETO, maxx, miny), draw(VT_MOVETO, maxx, maxy));
 		break;

 	case 2:  // bottom
-		return get_line_intersection(a, b, draw(VT_MOVETO, area + clip_buffer, area + clip_buffer), draw(VT_MOVETO, -clip_buffer, area + clip_buffer));
+		return get_line_intersection(a, b, draw(VT_MOVETO, maxx, maxy), draw(VT_MOVETO, minx, maxy));
 		break;

 	case 3:  // left
-		return get_line_intersection(a, b, draw(VT_MOVETO, -clip_buffer, area + clip_buffer), draw(VT_MOVETO, -clip_buffer, -clip_buffer));
+		return get_line_intersection(a, b, draw(VT_MOVETO, minx, maxy), draw(VT_MOVETO, minx, miny));
 		break;
 	}

@@ -402,14 +398,9 @@ static draw intersect(draw a, draw b, int edge, long long area, long long buffer
 }

 // http://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm
-static drawvec clip_poly1(drawvec &geom, int z, int detail, int buffer) {
+static drawvec clip_poly1(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) {
 	drawvec out = geom;

-	long long area = 0xFFFFFFFF;
-	if (z != 0) {
-		area = 1LL << (32 - z);
-	}
-
 	for (int edge = 0; edge < 4; edge++) {
 		if (out.size() > 0) {
 			drawvec in = out;
@@ -420,13 +411,13 @@ static drawvec clip_poly1(drawvec &geom, int z, int detail, int buffer) {
 			for (unsigned e = 0; e < in.size(); e++) {
 				draw E = in[e];

-				if (inside(E, edge, area, buffer)) {
-					if (!inside(S, edge, area, buffer)) {
-						out.push_back(intersect(S, E, edge, area, buffer));
+				if (inside(E, edge, minx, miny, maxx, maxy)) {
+					if (!inside(S, edge, minx, miny, maxx, maxy)) {
+						out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy));
 					}
 					out.push_back(E);
-				} else if (inside(S, edge, area, buffer)) {
-					out.push_back(intersect(S, E, edge, area, buffer));
+				} else if (inside(S, edge, minx, miny, maxx, maxy)) {
+					out.push_back(intersect(S, E, edge, minx, miny, maxx, maxy));
 				}

 				S = E;
@@ -445,7 +436,9 @@ static drawvec clip_poly1(drawvec &geom, int z, int detail, int buffer) {
 		}

 		if (out.size() < 3) {
-			fprintf(stderr, "Polygon degenerated to a line segment\n");
+			// fprintf(stderr, "Polygon degenerated to a line segment\n");
+			out.clear();
+			return out;
 		}

 		out[0].op = VT_MOVETO;
@@ -457,11 +450,7 @@ static drawvec clip_poly1(drawvec &geom, int z, int detail, int buffer) {
 	return out;
 }

-drawvec simple_clip_poly(drawvec &geom, int z, int detail, int buffer) {
-	if (z == 0) {
-		return geom;
-	}
-
+drawvec simple_clip_poly(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) {
 	drawvec out;

 	for (unsigned i = 0; i < geom.size(); i++) {
@@ -477,7 +466,7 @@ drawvec simple_clip_poly(drawvec &geom, int z, int detail, int buffer) {
 			for (unsigned k = i; k < j; k++) {
 				tmp.push_back(geom[k]);
 			}
-			tmp = clip_poly1(tmp, z, detail, buffer);
+			tmp = clip_poly1(tmp, minx, miny, maxx, maxy);
 			if (tmp.size() > 0) {
 				if (tmp[0].x != tmp[tmp.size() - 1].x || tmp[0].y != tmp[tmp.size() - 1].y) {
 					fprintf(stderr, "Internal error: Polygon ring not closed\n");
@@ -498,6 +487,17 @@ drawvec simple_clip_poly(drawvec &geom, int z, int detail, int buffer) {
 	return out;
 }

+drawvec simple_clip_poly(drawvec &geom, int z, int detail, int buffer) {
+	long long area = 0xFFFFFFFF;
+	if (z != 0) {
+		area = 1LL << (32 - z);
+	}
+
+	long long clip_buffer = buffer * area / 256;
+
+	return simple_clip_poly(geom, -clip_buffer, -clip_buffer, area + clip_buffer, area + clip_buffer);
+}
+
 drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double *accum_area) {
 	drawvec out;
 	long long pixel = (1 << (32 - detail - z)) * 2;
@@ -943,3 +943,63 @@ drawvec fix_polygon(drawvec &geom) {

 	return out;
 }
+
+std::vector<drawvec> chop_polygon(std::vector<drawvec> &geoms) {
+	while (1) {
+		bool again = false;
+		std::vector<drawvec> out;
+
+		for (unsigned i = 0; i < geoms.size(); i++) {
+			if (geoms[i].size() > 700) {
+				long long midx = 0, midy = 0, count = 0;
+				long long maxx = LONG_LONG_MIN, maxy = LONG_LONG_MIN, minx = LONG_LONG_MAX, miny = LONG_LONG_MAX;
+
+				for (unsigned j = 0; j < geoms[i].size(); j++) {
+					if (geoms[i][j].op == VT_MOVETO || geoms[i][j].op == VT_LINETO) {
+						midx += geoms[i][j].x;
+						midy += geoms[i][j].y;
+						count++;
+
+						if (geoms[i][j].x > maxx) {
+							maxx = geoms[i][j].x;
+						}
+						if (geoms[i][j].y > maxy) {
+							maxy = geoms[i][j].y;
+						}
+						if (geoms[i][j].x < minx) {
+							minx = geoms[i][j].x;
+						}
+						if (geoms[i][j].y < miny) {
+							miny = geoms[i][j].y;
+						}
+					}
+				}
+
+				midx /= count;
+				midy /= count;
+
+				if (maxy - miny > maxx - minx) {
+					// printf("clipping y to %lld %lld %lld %lld\n", minx, miny, maxx, midy);
+					out.push_back(simple_clip_poly(geoms[i], minx, miny, maxx, midy));
+					// printf("          and %lld %lld %lld %lld\n", minx, midy, maxx, maxy);
+					out.push_back(simple_clip_poly(geoms[i], minx, midy, maxx, maxy));
+				} else {
+					// printf("clipping x to %lld %lld %lld %lld\n", minx, miny, midx, maxy);
+					out.push_back(simple_clip_poly(geoms[i], minx, miny, midx, maxy));
+					// printf("          and %lld %lld %lld %lld\n", midx, midy, maxx, maxy);
+					out.push_back(simple_clip_poly(geoms[i], midx, miny, maxx, maxy));
+				}
+
+				again = true;
+			} else {
+				out.push_back(geoms[i]);
+			}
+		}
+
+		if (!again) {
+			return out;
+		}
+
+		geoms = out;
+	}
+}
@@ -29,3 +29,4 @@ int quick_check(long long *bbox, int z, int detail, long long buffer);
 drawvec simplify_lines(drawvec &geom, int z, int detail);
 drawvec reorder_lines(drawvec &geom);
 drawvec fix_polygon(drawvec &geom);
+std::vector<drawvec> chop_polygon(std::vector<drawvec> &geoms);
@@ -8,6 +8,8 @@ features. This is a tool for making maps from huge datasets
 .PP
 [Build Status](https://travis\-ci.org/mapbox/tippecanoe.svg)
 \[la]https://travis-ci.org/mapbox/tippecanoe\[ra]
+[Coverage Status](https://coveralls.io/repos/mapbox/tippecanoe/badge.svg?branch=master&service=github)
+\[la]https://coveralls.io/github/mapbox/tippecanoe?branch=master\[ra]
 .SH Intent
 .PP
 The goal of Tippecanoe is to enable making a scale\-independent view of your data,
@@ -143,6 +145,8 @@ If you use \-rg, it will guess a drop rate that will keep at most 50,000 feature
 .IP \(bu 2
 \-pi: Preserve the original input order of features as the drawing order instead of ordering geographically. (This is implemented as a restoration of the original order at the end, so that dot\-dropping is still geographic, which means it also undoes \-ao).
 .IP \(bu 2
+\-pp: Don't split complex polygons (over 700 vertices after simplification) into multiple features.
+.IP \(bu 2
 \-q: Work quietly instead of reporting progress
 .RE
 .SH Example
@@ -232,6 +236,9 @@ have their probability diffused, so that some of them will be drawn as a square
 this minimum size and others will not be drawn at all, preserving the total area that
 all of them should have had together.
 .PP
+Any polygons that have over 700 vertices after line simplification will be split into
+multiple features so they can be rendered efficiently, unless you use \-pp to prevent this.
+.PP
 Features in the same tile that share the same type and attributes are coalesced
 together into a single geometry. You are strongly encouraged to use \-x to exclude
 any unnecessary properties to reduce wasted file size.
@@ -155,7 +155,7 @@ void mbtiles_write_metadata(sqlite3 *outdb, const char *fname, char **layername,
 	}
 	sqlite3_free(sql);

-	sql = sqlite3_mprintf("INSERT INTO metadata (name, value) VALUES ('version', %d);", 1);
+	sql = sqlite3_mprintf("INSERT INTO metadata (name, value) VALUES ('version', %d);", 2);
 	if (sqlite3_exec(outdb, sql, NULL, NULL, &err) != SQLITE_OK) {
 		fprintf(stderr, "set version : %s\n", err);
 		if (!forcetable) {
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/highzoom/out/-z30.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "in" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/multilayer/out/-ltogether_-z3.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "together" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "separate",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "lines" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/ne_110m_admin_0_countries/out/-z4_-yname.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "in" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/ne_110m_admin_1_states_provinces_lines/out/-X_-z4.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "in" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "1",
 "name": "tests/ne_110m_admin_1_states_provinces_lines/out/-lcountries_-P_-Z1_-z7_-b4_-xfeaturecla_-xscalerank_-acrol_-ps.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 1, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "countries" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/ne_110m_admin_1_states_provinces_lines/out/-z5_-ymapcolor13_-ymapcolor9_-pSi_-d8_-D16.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "in" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/ne_110m_populated_places/out/-yNAME.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "in" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/ne_110m_populated_places/out/-yNAME_-z5.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "in" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/ne_110m_populated_places/out/-yNAME_-z5_-B3.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "in" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/ne_110m_populated_places/out/-yNAME_-z5_-r1.5.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "in" }, "features": [
@@ -8,7 +8,7 @@
 "minzoom": "0",
 "name": "tests/nullisland/out/-b0_-z4.json.check.mbtiles",
 "type": "overlay",
-"version": "1"
+"version": "2"
 }, "features": [
 { "type": "FeatureCollection", "properties": { "zoom": 0, "x": 0, "y": 0 }, "features": [
 { "type": "FeatureCollection", "properties": { "layer": "in" }, "features": [
@@ -468,7 +468,7 @@ void rewrite(drawvec &geom, int z, int nextzoom, int maxzoom, long long *bbox, u
 }

 struct partial {
-	drawvec geom;
+	std::vector<drawvec> geoms;
 	long long layer;
 	char *meta;
 	signed char t;
@@ -495,8 +495,8 @@ void *partial_feature_worker(void *v) {
 	std::vector<struct partial> *partials = a->partials;

 	for (unsigned i = a->task; i < (*partials).size(); i += a->tasks) {
-		drawvec geom = (*partials)[i].geom;
-		(*partials)[i].geom.clear();  // avoid keeping two copies in memory
+		drawvec geom = (*partials)[i].geoms[0];  // XXX assumption of a single geometry at the beginning
+		(*partials)[i].geoms.clear();		 // avoid keeping two copies in memory
 		signed char t = (*partials)[i].t;
 		int z = (*partials)[i].z;
 		int line_detail = (*partials)[i].line_detail;
@@ -526,17 +526,26 @@ void *partial_feature_worker(void *v) {

 		to_tile_scale(geom, z, line_detail);

+		std::vector<drawvec> geoms;
+		geoms.push_back(geom);
+
+		if (t == VT_POLYGON && !prevent[P_POLYGON_SPLIT]) {
+			geoms = chop_polygon(geoms);
+		}
+
 		if (t == VT_POLYGON) {
 			// Scaling may have made the polygon degenerate.
 			// Give Clipper a chance to try to fix it.
-			geom = clean_or_clip_poly(geom, 0, 0, 0, false);
-			geom = close_poly(geom);
+			for (unsigned i = 0; i < geoms.size(); i++) {
+				geoms[i] = clean_or_clip_poly(geoms[i], 0, 0, 0, false);
+				geoms[i] = close_poly(geoms[i]);
+			}
 		}

 		// Worth skipping this if not coalescing anyway?
-		if (geom.size() > 0) {
-			(*partials)[i].index = encode(geom[0].x, geom[0].y);
-			(*partials)[i].index2 = encode(geom[geom.size() - 1].x, geom[geom.size() - 1].y);
+		if (geoms.size() > 0 && geoms[0].size() > 0) {
+			(*partials)[i].index = encode(geoms[0][0].x, geoms[0][0].y);
+			(*partials)[i].index2 = encode(geoms[0][geoms[0].size() - 1].x, geoms[0][geoms[0].size() - 1].y);

 			// Anything numbered below the start of the line
 			// can't possibly be the next feature.
@@ -549,7 +558,7 @@ void *partial_feature_worker(void *v) {
 			(*partials)[i].index2 = 0;
 		}

-		(*partials)[i].geom = geom;
+		(*partials)[i].geoms = geoms;
 	}

 	return NULL;
@@ -801,7 +810,7 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi

 			if (geom.size() > 0) {
 				partial p;
-				p.geom = geom;
+				p.geoms.push_back(geom);
 				p.layer = layer;
 				p.meta = meta;
 				p.t = t;
@@ -851,27 +860,31 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi

 		// This is serial because decode_meta() unifies duplicates
 		for (unsigned i = 0; i < partials.size(); i++) {
-			drawvec geom = partials[i].geom;
-			partials[i].geom.clear();  // avoid keeping two copies in memory
+			std::vector<drawvec> geoms = partials[i].geoms;
+			partials[i].geoms.clear();  // avoid keeping two copies in memory
 			long long layer = partials[i].layer;
-			char *meta = partials[i].meta;
 			signed char t = partials[i].t;
 			int segment = partials[i].segment;
 			long long original_seq = partials[i].original_seq;

-			if (t == VT_POINT || to_feature(geom, NULL)) {
-				struct coalesce c;
+			// A complex polygon may have been split up into multiple geometries.
+			// Break them out into multiple features if necessary.
+			for (unsigned j = 0; j < geoms.size(); j++) {
+				if (t == VT_POINT || to_feature(geoms[j], NULL)) {
+					struct coalesce c;
+					char *meta = partials[i].meta;

-				c.type = t;
-				c.index = partials[i].index;
-				c.index2 = partials[i].index2;
-				c.geom = geom;
-				c.metasrc = meta;
-				c.coalesced = false;
-				c.original_seq = original_seq;
+					c.type = t;
+					c.index = partials[i].index;
+					c.index2 = partials[i].index2;
+					c.geom = geoms[j];
+					c.metasrc = meta;
+					c.coalesced = false;
+					c.original_seq = original_seq;

-				decode_meta(&meta, stringpool + pool_off[segment], keys[layer], values[layer], file_keys[layer], &c.meta);
-				features[layer].push_back(c);
+					decode_meta(&meta, stringpool + pool_off[segment], keys[layer], values[layer], file_keys[layer], &c.meta);
+					features[layer].push_back(c);
+				}
 			}
 		}

@@ -899,7 +912,7 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi
 				}
 #endif

-				if (additional[A_COALESCE] && out.size() > 0 && out[y].geom.size() + features[j][x].geom.size() < 20000 && coalcmp(&features[j][x], &out[y]) == 0 && features[j][x].type != VT_POINT) {
+				if (additional[A_COALESCE] && out.size() > 0 && out[y].geom.size() + features[j][x].geom.size() < 700 && coalcmp(&features[j][x], &out[y]) == 0 && features[j][x].type != VT_POINT) {
 					unsigned z;
 					for (z = 0; z < features[j][x].geom.size(); z++) {
 						out[y].geom.push_back(features[j][x].geom[z]);
@@ -909,6 +922,7 @@ long long write_tile(char **geoms, char *metabase, char *stringpool, int z, unsi
 					out.push_back(features[j][x]);
 				}
 			}
+
 			features[j] = out;

 			out.clear();
@@ -60,4 +60,6 @@ static int prevent_options[] = {
 	P_DYNAMIC_DROP,
 #define P_INPUT_ORDER ((int) 'i')
 	P_INPUT_ORDER,
+#define P_POLYGON_SPLIT ((int) 'p')
+	P_POLYGON_SPLIT,
 };
@@ -1 +1 @@
-#define VERSION "tippecanoe v1.8.1\n"
+#define VERSION "tippecanoe v1.9.0\n"