Add --clip-bounding-box option to clip input geometry

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 * Fix null pointer crash when reading filter output that does not
   tag features with their extent
+* Add `--clip-bounding-box` option to clip input geometry
 
 ## 1.32.0
 

README.md

@@ -447,6 +447,7 @@ the same layer, enclose them in an `all` expression so they will all be evaluate
  * `-aw` or `--detect-longitude-wraparound`: Detect when adjacent points within a feature jump to the other side of the world, and try to fix the geometry.
  * `-pw` or `--use-source-polygon-winding`: Instead of respecting GeoJSON polygon ring order, use the original polygon winding in the source data to distinguish inner (clockwise) and outer (counterclockwise) polygon rings.
  * `-pW` or `--reverse-source-polygon-winding`: Instead of respecting GeoJSON polygon ring order, use the opposite of the original polygon winding in the source data to distinguish inner (counterclockwise) and outer (clockwise) polygon rings.
+ * `--clip-bounding-box=`*minlon*`,`*minlat*`,`*maxlon*`,`*maxlat*: Clip all features to the specified bounding box.
 
 ### Setting or disabling tile size limits
 

geometry.cpp

@@ -595,16 +595,20 @@ drawvec reduce_tiny_poly(drawvec &geom, int z, int detail, bool *reduced, double
 }
 
 drawvec clip_point(drawvec &geom, int z, long long buffer) {
-	drawvec out;
-
 	long long min = 0;
 	long long area = 1LL << (32 - z);
 
 	min -= buffer * area / 256;
 	area += buffer * area / 256;
 
+	return clip_point(geom, min, min, area, area);
+}
+
+drawvec clip_point(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) {
+	drawvec out;
+
 	for (size_t i = 0; i < geom.size(); i++) {
-		if (geom[i].x >= min && geom[i].y >= min && geom[i].x <= area && geom[i].y <= area) {
+		if (geom[i].x >= minx && geom[i].y >= miny && geom[i].x <= maxx && geom[i].y <= maxy) {
 			out.push_back(geom[i]);
 		}
 	}
@@ -646,13 +650,17 @@ bool point_within_tile(long long x, long long y, int z) {
 }
 
 drawvec clip_lines(drawvec &geom, int z, long long buffer) {
-	drawvec out;
-
 	long long min = 0;
 	long long area = 1LL << (32 - z);
 	min -= buffer * area / 256;
 	area += buffer * area / 256;
 
+	return clip_lines(geom, min, min, area, area);
+}
+
+drawvec clip_lines(drawvec &geom, long long minx, long long miny, long long maxx, long long maxy) {
+	drawvec out;
+
 	for (size_t i = 0; i < geom.size(); i++) {
 		if (i > 0 && (geom[i - 1].op == VT_MOVETO || geom[i - 1].op == VT_LINETO) && geom[i].op == VT_LINETO) {
 			double x1 = geom[i - 1].x;
@@ -661,7 +669,7 @@ drawvec clip_lines(drawvec &geom, int z, long long buffer) {
 			double x2 = geom[i - 0].x;
 			double y2 = geom[i - 0].y;
 
-			int c = clip(&x1, &y1, &x2, &y2, min, min, area, area);
+			int c = clip(&x1, &y1, &x2, &y2, minx, miny, maxx, maxy);
 
 			if (c > 1) {  // clipped
 				out.push_back(draw(VT_MOVETO, x1, y1));

geometry.hpp

@@ -76,4 +76,8 @@ void check_polygon(drawvec &geom);
 double get_area(drawvec &geom, size_t i, size_t j);
 double get_mp_area(drawvec &geom);
 
+drawvec simple_clip_poly(drawvec &geom, long long x1, long long y1, long long x2, long long y2);
+drawvec clip_lines(drawvec &geom, long long x1, long long y1, long long x2, long long y2);
+drawvec clip_point(drawvec &geom, long long x1, long long y1, long long x2, long long y2);
+
 #endif

main.cpp

@@ -93,6 +93,8 @@ long long MAX_FILES;
 static long long diskfree;
 char **av;
 
+std::vector<clipbbox> clipbboxes;
+
 void checkdisk(std::vector<struct reader> *r) {
 	long long used = 0;
 	for (size_t i = 0; i < r->size(); i++) {
@@ -2594,6 +2596,7 @@ int main(int argc, char **argv) {
 		{"detect-longitude-wraparound", no_argument, &additional[A_DETECT_WRAPAROUND], 1},
 		{"use-source-polygon-winding", no_argument, &prevent[P_USE_SOURCE_POLYGON_WINDING], 1},
 		{"reverse-source-polygon-winding", no_argument, &prevent[P_REVERSE_SOURCE_POLYGON_WINDING], 1},
+		{"clip-bounding-box", required_argument, 0, '~'},
 
 		{"Filtering tile contents", 0, 0, 0},
 		{"prefilter", required_argument, 0, 'C'},
@@ -2682,6 +2685,17 @@ int main(int argc, char **argv) {
 				max_tilestats_sample_values = atoi(optarg);
 			} else if (strcmp(opt, "tile-stats-values-limit") == 0) {
 				max_tilestats_values = atoi(optarg);
+			} else if (strcmp(opt, "clip-bounding-box") == 0) {
+				clipbbox clip;
+				if (sscanf(optarg, "%lf,%lf,%lf,%lf", &clip.lon1, &clip.lat1, &clip.lon2, &clip.lat2) == 4) {
+					clipbboxes.push_back(clip);
+				} else {
+					fprintf(stderr, "%s: Can't parse bounding box --%s=%s\n", argv[0], opt, optarg);
+					exit(EXIT_FAILURE);
+				}
+			} else {
+				fprintf(stderr, "%s: Unrecognized option --%s\n", argv[0], opt);
+				exit(EXIT_FAILURE);
 			}
 			break;
 		}
@@ -3003,6 +3017,14 @@ int main(int argc, char **argv) {
 		}
 	}
 
+	// Wait until here to project the bounding box, so that the behavior is
+	// the same no matter what order the projection and bounding box are
+	// specified in
+	for (auto &c : clipbboxes) {
+		projection->project(c.lon1, c.lat1, 32, &c.minx, &c.maxy);
+		projection->project(c.lon2, c.lat2, 32, &c.maxx, &c.miny);
+	}
+
 	if (max_tilestats_sample_values < max_tilestats_values) {
 		max_tilestats_sample_values = max_tilestats_values;
 	}

main.hpp

@@ -18,6 +18,20 @@ struct index {
 	}
 };
 
+struct clipbbox {
+	double lon1;
+	double lat1;
+	double lon2;
+	double lat2;
+
+	long long minx;
+	long long miny;
+	long long maxx;
+	long long maxy;
+};
+
+extern std::vector<clipbbox> clipbboxes;
+
 void checkdisk(std::vector<struct reader> *r);
 
 extern int geometry_scale;

man/tippecanoe.1

@@ -560,6 +560,8 @@ the line or polygon within one tile unit of its proper location. You can probabl
 \fB\fC\-pw\fR or \fB\fC\-\-use\-source\-polygon\-winding\fR: Instead of respecting GeoJSON polygon ring order, use the original polygon winding in the source data to distinguish inner (clockwise) and outer (counterclockwise) polygon rings.
 .IP \(bu 2
 \fB\fC\-pW\fR or \fB\fC\-\-reverse\-source\-polygon\-winding\fR: Instead of respecting GeoJSON polygon ring order, use the opposite of the original polygon winding in the source data to distinguish inner (counterclockwise) and outer (clockwise) polygon rings.
+.IP \(bu 2
+\fB\fC\-\-clip\-bounding\-box=\fR\fIminlon\fP\fB\fC,\fR\fIminlat\fP\fB\fC,\fR\fImaxlon\fP\fB\fC,\fR\fImaxlat\fP: Clip all features to the specified bounding box.
 .RE
 .SS Setting or disabling tile size limits
 .RS

serial.cpp

@@ -387,6 +387,16 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf) {
 		sf.geometry = fix_polygon(sf.geometry);
 	}
 
+	for (auto &c : clipbboxes) {
+		if (sf.t == VT_POLYGON) {
+			sf.geometry = simple_clip_poly(sf.geometry, c.minx >> geometry_scale, c.miny >> geometry_scale, c.maxx >> geometry_scale, c.maxy >> geometry_scale);
+		} else if (sf.t == VT_LINE) {
+			sf.geometry = clip_lines(sf.geometry, c.minx >> geometry_scale, c.miny >> geometry_scale, c.maxx >> geometry_scale, c.maxy >> geometry_scale);
+		} else if (sf.t == VT_POINT) {
+			sf.geometry = clip_point(sf.geometry, c.minx >> geometry_scale, c.miny >> geometry_scale, c.maxx >> geometry_scale, c.maxy >> geometry_scale);
+		}
+	}
+
 	if (!sf.has_id) {
 		if (additional[A_GENERATE_IDS]) {
 			sf.has_id = true;