Add an option for clustering features

CHANGELOG.md

@@ -1,3 +1,8 @@
+## 1.27.5
+
+* Add --cluster-densest-as-needed to cluster features
+* Add --maximum-tile-features to set the maximum number of features in a tile
+
 ## 1.27.4
 
 * Support CSV point input

README.md

@@ -212,6 +212,7 @@ tippecanoe -z5 -o filtered.mbtiles -j '{ "ne_10m_admin_0_countries": [ "all", [
  * `-an` or `--drop-smallest-as-needed`: Dynamically drop the smallest features (physically smallest: the shortest lines or the smallest polygons) from each zoom level to keep large tiles under the 500K size limit. This option will not work for point features.
  * `-aN` or `--coalesce-smallest-as-needed`: Dynamically combine the smallest features (physically smallest: the shortest lines or the smallest polygons) from each zoom level into other nearby features to keep large tiles under the 500K size limit. This option will not work for point features, and will probably not help very much with LineStrings. It is mostly intended for polygons, to maintain the full original area covered by polygons while still reducing the feature count somehow. The attributes of the small polygons are *not* preserved into the combined features, only their geometry.
  * `-pd` or `--force-feature-limit`: 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. (This is like `-ad` but applies to each tile individually, not to the entire zoom level.) You probably don't want to use this.
+ * `-aC` or `--cluster-densest-as-needed`: If a tile is too large, try to reduce its size by increasing the minimum spacing between features, and leaving one placeholder feature from each group.  The remaining feature will be given a `"cluster": true` attribute to indicate that it represents a cluster and a `"point_count"` attribute to indicate the number of features that were clustered into it.
 
 ### Dropping tightly overlapping features
 
@@ -255,6 +256,7 @@ tippecanoe -z5 -o filtered.mbtiles -j '{ "ne_10m_admin_0_countries": [ "all", [
 ### Setting or disabling tile size limits
 
  * `-M` _bytes_ or `--maximum-tile-bytes=`_bytes_: Use the specified number of _bytes_ as the maximum compressed tile size instead of 500K.
+ * `-O` _features_ or `--maximum-tile-featuress=`_features_: Use the specified number of _features_ as the maximum in a tile instead of 200,000.
  * `-pf` or `--no-feature-limit`: Don't limit tiles to 200,000 features
  * `-pk` or `--no-tile-size-limit`: Don't limit tiles to 500K bytes
  * `-pC` or `--no-tile-compression`: Don't compress the PBF vector tile data.

main.cpp

@@ -68,6 +68,7 @@ int quiet_progress = 0;
 int geometry_scale = 0;
 double simplification = 1;
 size_t max_tile_size = 500000;
+size_t max_tile_features = 200000;
 
 int prevent[256];
 int additional[256];
@@ -2277,6 +2278,7 @@ int main(int argc, char **argv) {
 		{"drop-smallest-as-needed", no_argument, &additional[A_DROP_SMALLEST_AS_NEEDED], 1},
 		{"coalesce-smallest-as-needed", no_argument, &additional[A_COALESCE_SMALLEST_AS_NEEDED], 1},
 		{"force-feature-limit", no_argument, &prevent[P_DYNAMIC_DROP], 1},
+		{"cluster-densest-as-needed", no_argument, &additional[A_CLUSTER_DENSEST_AS_NEEDED], 1},
 
 		{"Dropping tightly overlapping features", 0, 0, 0},
 		{"gamma", required_argument, 0, 'g'},
@@ -2315,6 +2317,7 @@ int main(int argc, char **argv) {
 
 		{"Setting or disabling tile size limits", 0, 0, 0},
 		{"maximum-tile-bytes", required_argument, 0, 'M'},
+		{"maximum-tile-features", required_argument, 0, 'O'},
 		{"no-feature-limit", no_argument, &prevent[P_FEATURE_LIMIT], 1},
 		{"no-tile-size-limit", no_argument, &prevent[P_KILOBYTE_LIMIT], 1},
 		{"no-tile-compression", no_argument, &prevent[P_TILE_COMPRESSION], 1},
@@ -2601,6 +2604,10 @@ int main(int argc, char **argv) {
 			max_tile_size = atoll(optarg);
 			break;
 
+		case 'O':
+			max_tile_features = atoll(optarg);
+			break;
+
 		case 'c':
 			postfilter = optarg;
 			break;

main.hpp

@@ -27,6 +27,7 @@ extern size_t CPUS;
 extern size_t TEMP_FILES;
 
 extern size_t max_tile_size;
+extern size_t max_tile_features;
 
 int mkstemp_cloexec(char *name);
 FILE *fopen_oflag(const char *name, const char *mode, int oflag);

man/tippecanoe.1

@@ -246,6 +246,8 @@ compensate for the larger marker, or \fB\fC\-Bf\fR\fInumber\fP to allow at most
 \fB\fC\-aN\fR or \fB\fC\-\-coalesce\-smallest\-as\-needed\fR: Dynamically combine the smallest features (physically smallest: the shortest lines or the smallest polygons) from each zoom level into other nearby features to keep large tiles under the 500K size limit. This option will not work for point features, and will probably not help very much with LineStrings. It is mostly intended for polygons, to maintain the full original area covered by polygons while still reducing the feature count somehow. The attributes of the small polygons are \fInot\fP preserved into the combined features, only their geometry.
 .IP \(bu 2
 \fB\fC\-pd\fR or \fB\fC\-\-force\-feature\-limit\fR: 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. (This is like \fB\fC\-ad\fR but applies to each tile individually, not to the entire zoom level.) You probably don't want to use this.
+.IP \(bu 2
+\fB\fC\-aC\fR or \fB\fC\-\-cluster\-densest\-as\-needed\fR: If a tile is too large, try to reduce its size by increasing the minimum spacing between features, and leaving one placeholder feature from each group.  The remaining feature will be given a \fB\fC"cluster": true\fR attribute to indicate that it represents a cluster and a \fB\fC"point_count"\fR attribute to indicate the number of features that were clustered into it.
 .RE
 .SS Dropping tightly overlapping features
 .RS
@@ -287,9 +289,9 @@ the line or polygon within one tile unit of its proper location. You can probabl
 .IP \(bu 2
 \fB\fC\-pi\fR or \fB\fC\-\-preserve\-input\-order\fR: 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 \fB\fC\-ao\fR).
 .IP \(bu 2
-\fB\fC\-ao\fR or \fB\fC\-\-reorder\fR: Reorder features to put ones with the same properties in sequence, to try to get them to coalesce. You probably don't want to use this.
+\fB\fC\-ao\fR or \fB\fC\-\-reorder\fR: Reorder features to put ones with the same properties in sequence, to try to get them to coalesce. You probably want to use this if you use \fB\fC\-\-coalesce\fR\&.
 .IP \(bu 2
-\fB\fC\-ac\fR or \fB\fC\-\-coalesce\fR: Coalesce adjacent line and polygon features that have the same properties. You probably don't want to use this.
+\fB\fC\-ac\fR or \fB\fC\-\-coalesce\fR: Coalesce adjacent line and polygon features that have the same properties. This can be useful if you have lots of small polygons with identical attributes and you would like to merge them together.
 .IP \(bu 2
 \fB\fC\-ar\fR or \fB\fC\-\-reverse\fR: Try reversing the directions of lines to make them coalesce and compress better. You probably don't want to use this.
 .RE
@@ -308,6 +310,8 @@ the line or polygon within one tile unit of its proper location. You can probabl
 .IP \(bu 2
 \fB\fC\-M\fR \fIbytes\fP or \fB\fC\-\-maximum\-tile\-bytes=\fR\fIbytes\fP: Use the specified number of \fIbytes\fP as the maximum compressed tile size instead of 500K.
 .IP \(bu 2
+\fB\fC\-O\fR \fIfeatures\fP or \fB\fC\-\-maximum\-tile\-featuress=\fR\fIfeatures\fP: Use the specified number of \fIfeatures\fP as the maximum in a tile instead of 200,000.
+.IP \(bu 2
 \fB\fC\-pf\fR or \fB\fC\-\-no\-feature\-limit\fR: Don't limit tiles to 200,000 features
 .IP \(bu 2
 \fB\fC\-pk\fR or \fB\fC\-\-no\-tile\-size\-limit\fR: Don't limit tiles to 500K bytes

mbtiles.cpp

@@ -561,7 +561,7 @@ std::map<std::string, layermap_entry> merge_layermaps(std::vector<std::map<std::
 
 	for (size_t i = 0; i < maps.size(); i++) {
 		for (auto map = maps[i].begin(); map != maps[i].end(); ++map) {
-			if (map->second.points + map->second.lines + map->second.polygons == 0) {
+			if (map->second.points + map->second.lines + map->second.polygons + map->second.retain == 0) {
 				continue;
 			}
 

mbtiles.hpp

@@ -29,6 +29,7 @@ struct layermap_entry {
 	size_t points = 0;
 	size_t lines = 0;
 	size_t polygons = 0;
+	size_t retain = 0;  // keep for tilestats, even if no features directly here
 
 	layermap_entry(size_t _id) {
 		id = _id;

options.hpp

@@ -19,6 +19,7 @@
 #define A_GRID_LOW_ZOOMS ((int) 'L')
 #define A_DETECT_WRAPAROUND ((int) 'w')
 #define A_EXTEND_ZOOMS ((int) 'e')
+#define A_CLUSTER_DENSEST_AS_NEEDED ((int) 'C')
 
 #define P_SIMPLIFY ((int) 's')
 #define P_SIMPLIFY_LOW ((int) 'S')

serial.cpp

@@ -487,7 +487,7 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf) {
 	long long midy = (sf.bbox[1] / 2 + sf.bbox[3] / 2) & ((1LL << 32) - 1);
 	bbox_index = encode(midx, midy);
 
-	if (additional[A_DROP_DENSEST_AS_NEEDED] || additional[A_CALCULATE_FEATURE_DENSITY] || additional[A_INCREASE_GAMMA_AS_NEEDED] || sst->uses_gamma) {
+	if (additional[A_DROP_DENSEST_AS_NEEDED] || additional[A_CLUSTER_DENSEST_AS_NEEDED] || additional[A_CALCULATE_FEATURE_DENSITY] || additional[A_INCREASE_GAMMA_AS_NEEDED] || sst->uses_gamma) {
 		sf.index = bbox_index;
 	} else {
 		sf.index = 0;

tile.cpp

@@ -1385,6 +1385,34 @@ void *run_prefilter(void *v) {
 	return NULL;
 }
 
+void add_tilestats(std::string const &layername, int z, std::vector<std::map<std::string, layermap_entry>> *layermaps, size_t tiling_seg, std::vector<std::vector<std::string>> *layer_unmaps, std::string const &key, serial_val const &val) {
+	std::map<std::string, layermap_entry> &layermap = (*layermaps)[tiling_seg];
+	if (layermap.count(layername) == 0) {
+		layermap_entry lme = layermap_entry(layermap.size());
+		lme.minzoom = z;
+		lme.maxzoom = z;
+		lme.retain = 1;
+
+		layermap.insert(std::pair<std::string, layermap_entry>(layername, lme));
+
+		if (lme.id >= (*layer_unmaps)[tiling_seg].size()) {
+			(*layer_unmaps)[tiling_seg].resize(lme.id + 1);
+			(*layer_unmaps)[tiling_seg][lme.id] = layername;
+		}
+	}
+	auto fk = layermap.find(layername);
+	if (fk == layermap.end()) {
+		fprintf(stderr, "Internal error: layer %s not found\n", layername.c_str());
+		exit(EXIT_FAILURE);
+	}
+
+	type_and_string attrib;
+	attrib.type = val.type;
+	attrib.string = val.s;
+
+	add_to_file_keys(fk->second.file_keys, key, attrib);
+}
+
 long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *stringpool, int z, unsigned tx, unsigned ty, int detail, int min_detail, sqlite3 *outdb, const char *outdir, int buffer, const char *fname, FILE **geomfile, int minzoom, int maxzoom, double todo, volatile long long *along, long long alongminus, double gamma, int child_shards, long long *meta_off, long long *pool_off, unsigned *initial_x, unsigned *initial_y, volatile int *running, double simplification, std::vector<std::map<std::string, layermap_entry>> *layermaps, std::vector<std::vector<std::string>> *layer_unmaps, size_t tiling_seg, size_t pass, size_t passes, unsigned long long mingap, long long minextent, double fraction, const char *prefilter, const char *postfilter, write_tile_args *arg) {
 	int line_detail;
 	double merge_fraction = 1;
@@ -1429,6 +1457,7 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 		double scale = (double) (1LL << (64 - 2 * (z + 8)));
 		double gap = 0, density_gap = 0;
 		double spacing = 0;
+		size_t clustered = 0;
 
 		long long original_features = 0;
 		long long unclipped_features = 0;
@@ -1540,7 +1569,34 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 				}
 			}
 
-			if (additional[A_DROP_DENSEST_AS_NEEDED]) {
+			if (additional[A_CLUSTER_DENSEST_AS_NEEDED]) {
+				indices.push_back(sf.index);
+				if (sf.index - merge_previndex < mingap) {
+					clustered++;
+					continue;
+				} else {
+					if (clustered > 0) {
+						std::string layername = (*layer_unmaps)[sf.segment][sf.layer];
+						serial_val sv, sv2;
+
+						sf.full_keys.push_back("clustered");
+						sv.type = mvt_bool;
+						sv.s = "true";
+						sf.full_values.push_back(sv);
+
+						add_tilestats(layername, z, layermaps, tiling_seg, layer_unmaps, "clustered", sv);
+
+						sf.full_keys.push_back("point_count");
+						sv2.type = mvt_double;
+						sv2.s = std::to_string(clustered + 1);
+						sf.full_values.push_back(sv2);
+
+						add_tilestats(layername, z, layermaps, tiling_seg, layer_unmaps, "point_count", sv2);
+					}
+
+					clustered = 0;
+				}
+			} else if (additional[A_DROP_DENSEST_AS_NEEDED]) {
 				indices.push_back(sf.index);
 				if (sf.index - merge_previndex < mingap) {
 					continue;
@@ -1884,7 +1940,7 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 			fprintf(stderr, "Is your data in the wrong projection? It should be in WGS84/EPSG:4326.\n");
 		}
 
-		long long totalsize = 0;
+		size_t totalsize = 0;
 		for (auto layer_iterator = layers.begin(); layer_iterator != layers.end(); ++layer_iterator) {
 			std::vector<coalesce> &layer_features = layer_iterator->second;
 			totalsize += layer_features.size();
@@ -1899,11 +1955,11 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 		}
 
 		if (totalsize > 0 && tile.layers.size() > 0) {
-			if (totalsize > 200000 && !prevent[P_FEATURE_LIMIT]) {
-				fprintf(stderr, "tile %d/%u/%u has %lld features, >200000    \n", z, tx, ty, totalsize);
+			if (totalsize > max_tile_features && !prevent[P_FEATURE_LIMIT]) {
+				fprintf(stderr, "tile %d/%u/%u has %zu features, >%zu    \n", z, tx, ty, totalsize, max_tile_features);
 
 				if (has_polygons && additional[A_MERGE_POLYGONS_AS_NEEDED] && merge_fraction > .05 && merge_successful) {
-					merge_fraction = merge_fraction * 200000 / tile.layers.size() * 0.95;
+					merge_fraction = merge_fraction * max_tile_features / tile.layers.size() * 0.95;
 					if (!quiet) {
 						fprintf(stderr, "Going to try merging %0.2f%% of the polygons to make it fit\n", 100 - merge_fraction * 100);
 					}
@@ -1926,8 +1982,8 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 					}
 					line_detail++;  // to keep it the same when the loop decrements it
 					continue;
-				} else if (additional[A_DROP_DENSEST_AS_NEEDED]) {
-					mingap_fraction = mingap_fraction * 200000.0 / totalsize * 0.90;
+				} else if (additional[A_DROP_DENSEST_AS_NEEDED] || additional[A_CLUSTER_DENSEST_AS_NEEDED]) {
+					mingap_fraction = mingap_fraction * max_tile_features / totalsize * 0.90;
 					unsigned long long mg = choose_mingap(indices, mingap_fraction);
 					if (mg <= mingap) {
 						mg = mingap * 1.5;
@@ -1943,7 +1999,7 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 					line_detail++;
 					continue;
 				} else if (additional[A_DROP_SMALLEST_AS_NEEDED] || additional[A_COALESCE_SMALLEST_AS_NEEDED]) {
-					minextent_fraction = minextent_fraction * 200000.0 / totalsize * 0.90;
+					minextent_fraction = minextent_fraction * max_tile_features / totalsize * 0.90;
 					long long m = choose_minextent(extents, minextent_fraction);
 					if (m != minextent) {
 						minextent = m;
@@ -1958,7 +2014,7 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 						continue;
 					}
 				} else if (prevent[P_DYNAMIC_DROP] || additional[A_DROP_FRACTION_AS_NEEDED]) {
-					fraction = fraction * 200000 / totalsize * 0.95;
+					fraction = fraction * max_tile_features / totalsize * 0.95;
 					if (!quiet) {
 						fprintf(stderr, "Going to try keeping %0.2f%% of the features to make it fit\n", fraction * 100);
 					}
@@ -2010,7 +2066,7 @@ long long write_tile(FILE *geoms, long long *geompos_in, char *metabase, char *s
 						fprintf(stderr, "Going to try gamma of %0.3f to make it fit\n", gamma);
 					}
 					line_detail++;  // to keep it the same when the loop decrements it
-				} else if (additional[A_DROP_DENSEST_AS_NEEDED]) {
+				} else if (additional[A_DROP_DENSEST_AS_NEEDED] || additional[A_CLUSTER_DENSEST_AS_NEEDED]) {
 					mingap_fraction = mingap_fraction * max_tile_size / compressed.size() * 0.90;
 					unsigned long long mg = choose_mingap(indices, mingap_fraction);
 					if (mg <= mingap) {
@@ -2321,7 +2377,7 @@ int traverse_zooms(int *geomfd, off_t *geom_size, char *metabase, char *stringpo
 		int err = INT_MAX;
 
 		size_t start = 1;
-		if (additional[A_INCREASE_GAMMA_AS_NEEDED] || additional[A_DROP_DENSEST_AS_NEEDED] || additional[A_DROP_FRACTION_AS_NEEDED] || additional[A_DROP_SMALLEST_AS_NEEDED] || additional[A_COALESCE_SMALLEST_AS_NEEDED]) {
+		if (additional[A_INCREASE_GAMMA_AS_NEEDED] || additional[A_DROP_DENSEST_AS_NEEDED] || additional[A_CLUSTER_DENSEST_AS_NEEDED] || additional[A_DROP_FRACTION_AS_NEEDED] || additional[A_DROP_SMALLEST_AS_NEEDED] || additional[A_COALESCE_SMALLEST_AS_NEEDED]) {
 			start = 0;
 		}
 

version.hpp

@@ -1,6 +1,6 @@
 #ifndef VERSION_HPP
 #define VERSION_HPP
 
-#define VERSION "tippecanoe v1.27.4\n"
+#define VERSION "tippecanoe v1.27.5\n"
 
 #endif