Distinguish duplicate feature locations during the index sort

CHANGELOG.md

@@ -1,3 +1,7 @@
+# 2.81.0
+
+* Add --distinguish-duplicates option
+
 # 2.80.0
 
 * Remove undocumented command-line options

README.md

@@ -488,6 +488,7 @@ the same layer, enclose them in an `all` expression so they will all be evaluate
  * `-aS` or `--coalesce-fraction-as-needed`: Dynamically combine a fraction of features from each zoom level into other nearby features to keep large tiles under the 500K size limit. (Again, mostly useful for polygons.)
  * `-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 `"clustered": true` attribute to indicate that it represents a cluster, a `"point_count"` attribute to indicate the number of features that were clustered into it, and a `"sqrt_point_count"` attribute to indicate the relative width of a feature to represent the cluster. If the features being clustered are points, the representative feature will be located at the average of the original points' locations; otherwise, one of the original features will be left as the representative.
+ * `--distinguish-duplicates`: Treat features that share a location with some other feature as belonging to a series of sub-layers, so that they will be dropped or coalesced within their own sub-layer's feature sequence instead of being treated as infinitely dense and dropped or coalesced first. Up to 50 duplicates of each location are distinguished; any beyond that are all deferred to one final sub-layer.
 
 ### Dropping tightly overlapping features
 

main.cpp

@@ -348,7 +348,123 @@ int calc_feature_minzoom(struct index *ix, struct drop_state *ds, int maxzoom, d
 	return feature_minzoom;
 }
 
-static void merge(struct mergelist *merges, size_t nmerges, unsigned char *map, FILE *indexfile, int bytes, char *geom_map, FILE *geom_out, std::atomic<long long> *geompos, long long *progress, long long *progress_max, long long *progress_reported, int maxzoom, double gamma, struct drop_state *ds) {
+// The features that share a location with some other feature are deferred to
+// a later pass through the feature sequence, so that they are not treated as
+// infinitely dense and therefore dropped or coalesced before any of the
+// features whose locations are distinct.
+//
+// The number of duplicates that are distinguished at each location is limited,
+// because with a large number of features that all share the same location,
+// some of them will still have to be dropped or coalesced in some sequence,
+// and because each pass costs a temporary file. The duplicates beyond the
+// limit are all deferred to one final pass.
+#define MAX_DUPLICATES_DISTINGUISHED 50
+
+struct duplicate_deferral {
+	const char *tmpdir = NULL;
+
+	// The greatest number of extra passes that features can be deferred to,
+	// or 0 if duplicate locations are not being distinguished at all
+	size_t max_passes = 0;
+
+	// One temporary file per pass beyond the first: pass n is files[n - 1]
+	std::vector<FILE *> files;
+
+	// State for detecting runs of features that share the same index
+	unsigned long long previndex = 0;
+	bool have_previndex = false;
+	size_t pass = 0;
+};
+
+// Return the temporary file that features deferred to the specified pass
+// are accumulated in, creating it if this is the first feature to be
+// deferred that far.
+static FILE *deferral_file(struct duplicate_deferral *dd, size_t pass) {
+	while (dd->files.size() < pass) {
+		std::string name = std::string(dd->tmpdir) + "/dup.XXXXXXXX";
+		std::vector<char> tmpl(name.begin(), name.end());
+		tmpl.push_back('\0');
+
+		int fd = mkstemp_cloexec(tmpl.data());
+		if (fd < 0) {
+			perror("temporary file for duplicate features");
+			exit(EXIT_OPEN);
+		}
+		unlink(tmpl.data());
+
+		FILE *fp = fdopen(fd, "wb+");
+		if (fp == NULL) {
+			perror("reopen temporary file for duplicate features");
+			exit(EXIT_OPEN);
+		}
+
+		dd->files.push_back(fp);
+	}
+
+	return dd->files[pass - 1];
+}
+
+// Write one feature, and the index record that locates it, to the sorted
+// output. The features arrive here in index order, so a feature whose index
+// matches the previous feature's is a duplicate location, and is set aside
+// to be appended to the output later by replay_deferrals().
+static void write_feature(struct index ix, char *geom_map, FILE *geom_out, std::atomic<long long> *geompos, FILE *indexfile, int maxzoom, double gamma, struct drop_state *ds, struct duplicate_deferral *dd) {
+	if (dd->max_passes > 0) {
+		if (dd->have_previndex && ix.ix == dd->previndex) {
+			if (dd->pass < dd->max_passes) {
+				dd->pass++;
+			}
+		} else {
+			dd->pass = 0;
+		}
+
+		dd->previndex = ix.ix;
+		dd->have_previndex = true;
+	}
+
+	// MAGIC: This knows that the feature minzoom is the last byte of the serialized feature
+	// and is writing one byte less and then adding the byte for the minzoom.
+
+	long long len = ix.end - ix.start - 1;
+	int feature_minzoom = calc_feature_minzoom(&ix, ds, maxzoom, gamma);
+
+	if (dd->pass > 0) {
+		FILE *fp = deferral_file(dd, dd->pass);
+
+		// The index record has to be written now, to keep the index in the
+		// same order that the features were read in, but its location will
+		// not be known until the geometry is appended to the output, so
+		// remember where to come back and correct it.
+
+		long long indexoff = ftello(indexfile);
+		if (indexoff < 0) {
+			perror("ftello index");
+			exit(EXIT_SEEK);
+		}
+
+		unsigned char minzoom_byte = feature_minzoom;
+		if (fwrite(&indexoff, sizeof(indexoff), 1, fp) != 1 ||
+		    fwrite(&ix, sizeof(ix), 1, fp) != 1 ||
+		    fwrite(geom_map + ix.start, 1, len, fp) != (size_t) len ||
+		    fwrite(&minzoom_byte, 1, 1, fp) != 1) {
+			perror("write duplicate feature");
+			exit(EXIT_WRITE);
+		}
+	} else {
+		long long pos = *geompos;
+
+		fwrite_check(geom_map + ix.start, 1, len, geom_out, geompos, "merge geometry");
+		serialize_byte(geom_out, feature_minzoom, geompos, "merge geometry");
+
+		ix.start = pos;
+		ix.end = *geompos;
+	}
+
+	std::atomic<long long> indexpos(0);
+	fwrite_check(&ix, sizeof(struct index), 1, indexfile, &indexpos, "merge temporary");
+}
+
+static void merge(struct mergelist *merges, size_t nmerges, unsigned char *map, FILE *indexfile, int bytes, char *geom_map, FILE *geom_out, std::atomic<long long> *geompos, long long *progress, long long *progress_max, long long *progress_reported, int maxzoom, double gamma, struct drop_state *ds, struct duplicate_deferral *dd) {
 	struct mergelist *head = NULL;
 
 	for (size_t i = 0; i < nmerges; i++) {
@@ -361,14 +477,8 @@ static void merge(struct mergelist *merges, size_t nmerges, unsigned char *map,
 
 	while (head != NULL) {
 		struct index ix = *((struct index *) (map + head->start));
-		long long pos = *geompos;
-
-		// MAGIC: This knows that the feature minzoom is the last byte of the serialized feature
-		// and is writing one byte less and then adding the byte for the minzoom.
 
-		fwrite_check(geom_map + ix.start, 1, ix.end - ix.start - 1, geom_out, geompos, "merge geometry");
-		int feature_minzoom = calc_feature_minzoom(&ix, ds, maxzoom, gamma);
-		serialize_byte(geom_out, feature_minzoom, geompos, "merge geometry");
+		write_feature(ix, geom_map, geom_out, geompos, indexfile, maxzoom, gamma, ds, dd);
 
 		// Count this as an 75%-accomplishment, since we already 25%-counted it
 		*progress += (ix.end - ix.start) * 3 / 4;
@@ -378,10 +488,6 @@ static void merge(struct mergelist *merges, size_t nmerges, unsigned char *map,
 			*progress_reported = 100 * *progress / *progress_max;
 		}
 
-		ix.start = pos;
-		ix.end = *geompos;
-		std::atomic<long long> indexpos;
-		fwrite_check(&ix, bytes, 1, indexfile, &indexpos, "merge temporary");
 		head->start += bytes;
 
 		struct mergelist *m = head;
@@ -394,6 +500,71 @@ static void merge(struct mergelist *merges, size_t nmerges, unsigned char *map,
 	}
 }
 
+// Append the features that were deferred because they shared a location with
+// some other feature, one pass at a time, and correct the index records that
+// were written for them with the locations they have ended up at.
+static void replay_deferrals(struct duplicate_deferral *dd, FILE *geom_out, std::atomic<long long> *geompos, FILE *indexfile) {
+	if (dd->files.size() == 0) {
+		return;
+	}
+
+	// so that the corrections below aren't overwritten by buffered index writes
+	if (fflush(indexfile) != 0) {
+		perror("flush index");
+		exit(EXIT_WRITE);
+	}
+	int indexfd = fileno(indexfile);
+
+	for (size_t i = 0; i < dd->files.size(); i++) {
+		FILE *fp = dd->files[i];
+
+		if (fseeko(fp, 0, SEEK_SET) < 0) {
+			perror("rewind duplicate features");
+			exit(EXIT_SEEK);
+		}
+
+		while (true) {
+			long long indexoff;
+			struct index ix;
+
+			if (fread(&indexoff, sizeof(indexoff), 1, fp) != 1) {
+				if (ferror(fp)) {
+					perror("read duplicate features");
+					exit(EXIT_READ);
+				}
+				break;  // end of this pass
+			}
+			if (fread(&ix, sizeof(ix), 1, fp) != 1) {
+				fprintf(stderr, "Short read of duplicate feature index\n");
+				exit(EXIT_READ);
+			}
+
+			std::string s;
+			s.resize(ix.end - ix.start);
+			if (s.size() != 0 && fread(&s[0], 1, s.size(), fp) != s.size()) {
+				fprintf(stderr, "Short read of duplicate feature geometry\n");
+				exit(EXIT_READ);
+			}
+
+			long long pos = *geompos;
+			fwrite_check(s.data(), 1, s.size(), geom_out, geompos, "deferred geometry");
+			ix.start = pos;
+			ix.end = *geompos;
+			if (pwrite(indexfd, &ix, sizeof(ix), indexoff) != (ssize_t) sizeof(ix)) {
+				perror("correct duplicate feature index");
+				exit(EXIT_WRITE);
+			}
+		}
+
+		if (fclose(fp) != 0) {
+			perror("close duplicate features");
+			exit(EXIT_CLOSE);
+		}
+	}
+
+	dd->files.clear();
+}
+
 struct sort_arg {
 	int task;
 	int cpus;
@@ -741,7 +912,7 @@ void start_parsing(int fd, STREAM *fp, long long offset, long long len, std::ato
 	parser_created = true;
 }
 
-void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int splits, long long mem, const char *tmpdir, long long *availfiles, FILE *geomfile, FILE *indexfile, std::atomic<long long> *geompos_out, long long *progress, long long *progress_max, long long *progress_reported, int maxzoom, int basezoom, double droprate, double gamma, struct drop_state *ds) {
+void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int splits, long long mem, const char *tmpdir, long long *availfiles, FILE *geomfile, FILE *indexfile, std::atomic<long long> *geompos_out, long long *progress, long long *progress_max, long long *progress_reported, int maxzoom, int basezoom, double droprate, double gamma, struct drop_state *ds, struct duplicate_deferral *dd) {
 	// Arranged as bits to facilitate subdividing again if a subdivided file is still huge
 	int splitbits = log(splits) / log(2);
 	splits = 1 << splitbits;
@@ -960,7 +1131,7 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 				madvise(geommap, geomst.st_size, MADV_RANDOM);
 				madvise(geommap, geomst.st_size, MADV_WILLNEED);
 
-				merge(merges, nmerges, (unsigned char *) indexmap, indexfile, bytes, geommap, geomfile, geompos_out, progress, progress_max, progress_reported, maxzoom, gamma, ds);
+				merge(merges, nmerges, (unsigned char *) indexmap, indexfile, bytes, geommap, geomfile, geompos_out, progress, progress_max, progress_reported, maxzoom, gamma, ds, dd);
 
 				madvise(indexmap, indexst.st_size, MADV_DONTNEED);
 				if (munmap(indexmap, indexst.st_size) < 0) {
@@ -991,11 +1162,8 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 
 				for (size_t a = 0; a < indexst.st_size / sizeof(struct index); a++) {
 					struct index ix = indexmap[a];
-					long long pos = *geompos_out;
 
-					fwrite_check(geommap + ix.start, ix.end - ix.start, 1, geomfile, geompos_out, "geom");
-					int feature_minzoom = calc_feature_minzoom(&ix, ds, maxzoom, gamma);
-					serialize_byte(geomfile, feature_minzoom, geompos_out, "merge geometry");
+					write_feature(ix, geommap, geomfile, geompos_out, indexfile, maxzoom, gamma, ds, dd);
 
 					// Count this as an 75%-accomplishment, since we already 25%-counted it
 					*progress += (ix.end - ix.start) * 3 / 4;
@@ -1004,11 +1172,6 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 						fflush(stderr);
 						*progress_reported = 100 * *progress / *progress_max;
 					}
-
-					ix.start = pos;
-					ix.end = *geompos_out;
-					std::atomic<long long> indexpos;
-					fwrite_check(&ix, sizeof(struct index), 1, indexfile, &indexpos, "index");
 				}
 
 				madvise(indexmap, indexst.st_size, MADV_DONTNEED);
@@ -1029,7 +1192,7 @@ void radix1(int *geomfds_in, int *indexfds_in, int inputs, int prefix, int split
 				// counter backward but will be an honest estimate of the work remaining.
 				*progress_max += geomst.st_size / 4;
 
-				radix1(&geomfds[i], &indexfds[i], 1, prefix + splitbits, *availfiles / 4, mem, tmpdir, availfiles, geomfile, indexfile, geompos_out, progress, progress_max, progress_reported, maxzoom, basezoom, droprate, gamma, ds);
+				radix1(&geomfds[i], &indexfds[i], 1, prefix + splitbits, *availfiles / 4, mem, tmpdir, availfiles, geomfile, indexfile, geompos_out, progress, progress_max, progress_reported, maxzoom, basezoom, droprate, gamma, ds, dd);
 				already_closed = 1;
 			}
 		}
@@ -1087,6 +1250,24 @@ void radix(std::vector<struct reader> &readers, int nreaders, FILE *geomfile, FI
 			       - 4			 // top-level geom and index output, both FILE and fd
 			       - 3;			 // stdin, stdout, stderr
 
+	struct duplicate_deferral dd;
+	dd.tmpdir = tmpdir;
+
+	if (additional[A_DISTINGUISH_DUPLICATES]) {
+		// one file for each pass that duplicate features can be deferred to,
+		// including the final pass for the ones beyond the limit,
+		// but not so many that there are none left to sort with
+		dd.max_passes = MAX_DUPLICATES_DISTINGUISHED + 1;
+		if ((long long) dd.max_passes > availfiles - 8) {
+			dd.max_passes = std::max(0LL, availfiles - 8);
+
+			if (dd.max_passes == 0) {
+				fprintf(stderr, "Warning: not enough available files to distinguish duplicate locations\n");
+			}
+		}
+		availfiles -= dd.max_passes;
+	}
+
 	// 4 because for each we have output and input FILE and fd for geom and index
 	int splits = availfiles / 4;
 
@@ -1114,12 +1295,14 @@ void radix(std::vector<struct reader> &readers, int nreaders, FILE *geomfile, FI
 
 	long long progress = 0, progress_max = geom_total, progress_reported = -1;
 	long long availfiles_before = availfiles;
-	radix1(geomfds, indexfds, nreaders, 0, splits, mem, tmpdir, &availfiles, geomfile, indexfile, geompos, &progress, &progress_max, &progress_reported, maxzoom, basezoom, droprate, gamma, ds);
+	radix1(geomfds, indexfds, nreaders, 0, splits, mem, tmpdir, &availfiles, geomfile, indexfile, geompos, &progress, &progress_max, &progress_reported, maxzoom, basezoom, droprate, gamma, ds, &dd);
 
 	if (availfiles - 2 * nreaders != availfiles_before) {
 		fprintf(stderr, "Internal error: miscounted available file descriptors: %lld vs %lld\n", availfiles - 2 * nreaders, availfiles);
 		exit(EXIT_IMPOSSIBLE);
 	}
+
+	replay_deferrals(&dd, geomfile, geompos, indexfile);
 }
 
 void choose_first_zoom(long long *file_bbox, long long *file_bbox1, long long *file_bbox2, std::vector<struct reader> &readers, unsigned *iz, unsigned *ix, unsigned *iy, int minzoom, int buffer) {
@@ -2716,7 +2899,9 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
 			}
 		} else {
 			for (long long ip = 0; ip < indices; ip++) {
-				if (ip > 0 && map[ip].start != map[ip - 1].end) {
+				// The features are not consecutive in the geometry if some of them
+				// were deferred to a later pass for sharing a location with another
+				if (ip > 0 && map[ip].start != map[ip - 1].end && !additional[A_DISTINGUISH_DUPLICATES]) {
 					fprintf(stderr, "Mismatched index at %lld: %lld vs %lld\n", ip, map[ip].start, map[ip].end);
 				}
 				int feature_minzoom = calc_feature_minzoom(&map[ip], ds, maxzoom, gamma);
@@ -3100,6 +3285,7 @@ int main(int argc, char **argv) {
 		{"force-feature-limit", no_argument, &prevent[P_DYNAMIC_DROP], 1},
 		{"cluster-densest-as-needed", no_argument, &additional[A_CLUSTER_DENSEST_AS_NEEDED], 1},
 		{"keep-point-cluster-position", no_argument, &additional[A_KEEP_POINT_CLUSTER_POSITION], 1},
+		{"distinguish-duplicates", no_argument, &additional[A_DISTINGUISH_DUPLICATES], 1},
 
 		{"Dropping tightly overlapping features", 0, 0, 0},
 		{"gamma", required_argument, 0, 'g'},

man/tippecanoe.1

@@ -621,6 +621,8 @@ preference to retaining points in sparse areas and dropping points in dense area
 \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"clustered": true\fR attribute to indicate that it represents a cluster, a \fB\fC"point_count"\fR attribute to indicate the number of features that were clustered into it, and a \fB\fC"sqrt_point_count"\fR attribute to indicate the relative width of a feature to represent the cluster. If the features being clustered are points, the representative feature will be located at the average of the original points' locations; otherwise, one of the original features will be left as the representative.
+.IP \(bu 2
+\fB\fC\-\-distinguish\-duplicates\fR: Treat features that share a location with some other feature as belonging to a series of sub\-layers, so that they will be dropped or coalesced within their own sub\-layer's feature sequence instead of being treated as infinitely dense and dropped or coalesced first. Up to 50 duplicates of each location are distinguished; any beyond that are all deferred to one final sub\-layer.
 .RE
 .SS Dropping tightly overlapping features
 .RS

options.hpp

@@ -29,6 +29,7 @@
 #define A_VISVALINGAM ((int) 'v')
 #define A_DETECT_WRAPAROUND ((int) 'w')
 #define A_KEEP_POINT_CLUSTER_POSITION ((int) 'a')
+#define A_DISTINGUISH_DUPLICATES ((int) 'u')
 #define A_DROP_BY_ATTRIBUTE_AS_NEEDED ((int) 'A')
 
 #define P_TILE_COMPRESSION ((int) 'C')

serial.cpp

@@ -743,6 +743,7 @@ int serialize_feature(struct serialization_state *sst, serial_feature &sf, std::
 	    prevent[P_DYNAMIC_DROP] ||
 	    additional[A_INCREASE_GAMMA_AS_NEEDED] ||
 	    additional[A_GENERATE_POLYGON_LABEL_POINTS] ||
+	    additional[A_DISTINGUISH_DUPLICATES] ||
 	    sst->uses_gamma ||
 	    retain_points_multiplier > 1 ||
 	    preserve_multiplier_density_threshold > 0 ||