Declutter point features (cities/peaks) at low zoom

CMakeLists.txt

@@ -136,6 +136,7 @@ file(GLOB tilemaker_src_files
 	src/config_validator.cpp
 	src/coordinates.cpp
 	src/coordinates_geom.cpp
+	src/declutter.cpp
 	src/external/streamvbyte_decode.c
 	src/external/streamvbyte_encode.c
 	src/external/streamvbyte_zigzag.c

Makefile

@@ -100,6 +100,7 @@ tilemaker: \
 	src/config_validator.o \
 	src/coordinates_geom.o \
 	src/coordinates.o \
+	src/declutter.o \
 	src/external/streamvbyte_decode.o \
 	src/external/streamvbyte_encode.o \
 	src/external/streamvbyte_zigzag.o \

docs/CONFIGURATION.md

@@ -86,6 +86,9 @@ You can add optional parameters to layers:
 * `combine_lines_below` - whether to merge all linestrings in the tile with the same attributes. If not defined, the global setting `combine_below` will be used.
 * `combine_points` - merge points with the same attributes (defaults to `true`: specify `false` to disable)
 * `z_order_ascending` - sort features in ascending order by a numeric value set in the Lua processing script (defaults to `true`: specify `false` for descending order)
+* `declutter_below` - thin out point features below this zoom level, so only the most important appear when zoomed out (see 'Decluttering point features' below)
+* `declutter_distance` - how far apart (in pixels) decluttered features should be kept - defaults to 40
+* `declutter_threshold` - the score a feature needs to appear at the layer's `minzoom`; this halves at each subsequent zoom level - defaults to 0 (no threshold)
 
 `write_to` enables you to combine different layer specs within one outputted layer. For example:
 
@@ -98,6 +101,30 @@ This would combine the `roads` (z12-14) and `low_roads` (z9-11) layers into a si
 
 (See also 'Shapefiles and GeoJSON' below.)
 
+### Decluttering point features
+
+For more effective maps at lower zoom levels, tilemaker allows you to prioritise the most important point features - such as the largest cities or the highest peaks.
+
+In your Lua profile, give each point a `Score()`, typically calculated from OSM tags such as `population` or `ele`. Then set three `declutter` properties on the layer (see above) to determine which zoom levels this applies at, and the distance between each selected feature.
+
+    "place": {
+      "minzoom": 4, "maxzoom": 14,
+      "declutter_below": 11, "declutter_distance": 40, "declutter_threshold": 50000
+    }
+
+```lua
+    local place = Find("place")
+    if place=="city" or place=="town" or place=="village" then
+      Layer("place", false)
+      Attribute("name", Find("name"))
+      Score(tonumber(Find("population")) or 0)
+    end
+```
+
+Starting at the layer's `minzoom`, tilemaker takes the features whose score is at least `declutter_threshold` (highest score first), and places each one that isn't within `declutter_distance` of a feature already placed. It then moves up a zoom level, halving the threshold so that less important features are included. Anything still unplaced by the time we reach `declutter_below` is written from that zoom level upwards.
+
+A feature whose `MinZoom()` is above the zoom being considered isn't eligible for it. Only point features are decluttered, not lines or polygons.
+
 ### Additional metadata
 
 Tilemaker writes a `json` metadata field containing a `vector_layers` key, whose value is an array of JSON objects describing each layer and its attributes. This is part of the MBTiles 1.3 spec and required by certain clients.
@@ -171,6 +198,7 @@ To do that, you use these methods:
 * `IsMultiPolygon()`: returns true if the current object is a multipolygon.
 * `ZOrder(number)`: Set a numeric value (default 0) used to sort features within a layer. Use this feature to ensure a proper rendering order if the rendering engine itself does not support sorting. Sorting is not supported across layers merged with `write_to`. Features with different z-order are not merged if `combine_below`, `combine_lines_below` or `combine_polygons_below` is used. Use this in conjunction with `feature_limit` to only write the most important (highest z-order) features within a tile. (Values can be -50,000,000 to 50,000,000 and are lossy, particularly beyond -1000 to 1000.)
 * `MinZoom(zoom)`: set the minimum zoom level (0-15) at which this object will be written. Note that the JSON layer configuration minimum still applies (so `:MinZoom(5)` will have no effect if your layer only starts at z6).
+* `Score(number)`: set how important this point feature is (a whole number; anything else is rounded down), so that tilemaker can decide which zoom level to show it from. Only meaningful for points (i.e. features written with `Layer` for a node, or `LayerAsCentroid`) in a layer with `declutter_below` set - see 'Decluttering point features' below.
 * `Length()` and `Area()`: return the length (metres)/area (square metres) of the current object. Requires Boost 1.67+.
 * `Centroid()`: return the lat/lon of the centre of the current object as a two-element Lua table (element 1 is lat, 2 is lon).
 
@@ -233,6 +261,8 @@ Limited Lua transformations are available for these files. You can supply an `at
 
 To set the minimum zoom level at which an individual feature is rendered, use `attribute_function` to set a `_minzoom` value in your return table.
 
+Similarly, to declutter shapefile/GeoJSON points, set a `_score` value in your return table and `declutter_below` on the layer.
+
 Shapefiles/GeoJSON **must** be in WGS84 projection, i.e. pure latitude/longitude. (Use ogr2ogr to reproject them if your source material is in a different projection.) They will be clipped to the bounds of the first .pbf that you import, unless you specify otherwise with a `bounding_box` setting in your JSON file.
 
 ### Lua spatial queries

include/declutter.h

@@ -0,0 +1,56 @@
+/*! \file */ 
+#ifndef _DECLUTTER_H
+#define _DECLUTTER_H
+
+#include <mutex>
+#include <vector>
+#include "coordinates.h"
+#include "output_object.h"
+
+class TileDataSource;
+struct LayerDef;
+
+/// A point feature held back from the tile index until its minimum zoom has been decided
+struct DeclutterEntry {
+	OutputObject oo;
+	LatpLon point;
+	uint64_t id;
+	int32_t score;
+	bool fromShapefile;
+};
+
+/**
+	\brief Thins out point features so only the most important ones appear at low zooms.
+
+	Layers with declutter_below set don't index their point features as they're read: the
+	features are parked here, together with the score their Lua profile gave them with
+	Score(). Once everything has been read, apply() works up through the zoom levels,
+	placing the highest-scoring features first and holding back any that fall too close to
+	one already placed, then hands them all to the tile index.
+
+	Ranking globally rather than per z6 tile means features either side of a z6 boundary
+	still compete with each other.
+*/
+class Declutter {
+
+public:
+	/// Note which layers are decluttered - call before any add()
+	void configure(const std::vector<LayerDef>& layers);
+
+	bool inUse() const { return anyLayers; }
+	bool isDecluttered(uint_least8_t layer) const { return anyLayers && layerEnabled[layer]; }
+
+	/// Park a point feature (thread-safe)
+	void add(const OutputObject& oo, LatpLon point, uint64_t id, int32_t score, bool fromShapefile);
+
+	/// Assign minimum zooms, then write everything to the tile index
+	void apply(const std::vector<LayerDef>& layers, TileDataSource& osmSource, TileDataSource& shpSource);
+
+private:
+	bool anyLayers = false;
+	std::vector<bool> layerEnabled;
+	std::vector<std::vector<DeclutterEntry>> entries;	// one per layer
+	std::vector<std::mutex> mutexes;					//  |
+};
+
+#endif //_DECLUTTER_H

include/geojson_processor.h

@@ -43,7 +43,7 @@ class GeoJSONProcessor {
 	template <bool Flag, typename T>
 	std::vector<Point> pointsFromGeoJSONArray(const rapidjson::GenericArray<Flag, T> &arr);
 
-	AttributeIndex readProperties(const rapidjson::Value &pr, bool &hasName, std::string &name, LayerDef &layer, unsigned &minzoom);
+	AttributeIndex readProperties(const rapidjson::Value &pr, bool &hasName, std::string &name, LayerDef &layer, unsigned &minzoom, int32_t &score);
 };
 
 #endif //_GEOJSON_PROCESSOR_H

include/osm_lua_processing.h

@@ -11,6 +11,7 @@
 #include "osm_store.h"
 #include "shared_data.h"
 #include "output_object.h"
+#include "declutter.h"
 #include "shp_mem_tiles.h"
 #include "osm_mem_tiles.h"
 #include "helpers.h"
@@ -57,6 +58,7 @@ class OsmLuaProcessing {
 		const class ShpMemTiles &shpMemTiles, 
 		class OsmMemTiles &osmMemTiles,
 		AttributeStore &attributeStore,
+		class Declutter &declutter,
 		bool materializeGeometries,
 		bool isFirst
 	);
@@ -204,6 +206,7 @@ class OsmLuaProcessing {
 	void AttributeInteger(const std::string &key, const int val, const char minzoom);
 	void MinZoom(const double z);
 	void ZOrder(const double z);
+	void Score(const double score);
 
 	// Relation scan support
 
@@ -265,6 +268,7 @@ class OsmLuaProcessing {
 		relationAccepted = false;
 		relationList.clear();
 		relationSubscript = -1;
+		declutterOutputs.clear();
 		lastStoredGeometryId = 0;
 		isWay = false;
 		isRelation = false;
@@ -273,6 +277,8 @@ class OsmLuaProcessing {
 
 	void removeAttributeIfNeeded(const std::string& key);
 
+	void noteIfDecluttered(LatpLon point);
+
 	const inline Point getPoint() {
 		return Point(lon/10000000.0,latp/10000000.0);
 	}
@@ -291,6 +297,7 @@ class OsmLuaProcessing {
 	const class ShpMemTiles &shpMemTiles;
 	class OsmMemTiles &osmMemTiles;
 	AttributeStore &attributeStore;			// key/value store
+	class Declutter &declutter;				// point features held back for decluttering
 
 	int64_t originalOsmID;					///< Original OSM object ID
 	bool isWay, isRelation, isClosed;		///< Way, node, relation?
@@ -321,6 +328,11 @@ class OsmLuaProcessing {
 	class LayerDefinition &layers;
 
 	std::vector<std::pair<OutputObject, AttributeSet>> outputs;		// All output objects that have been created
+
+	// Point features in decluttered layers: the exact location we'd have indexed them at,
+	// plus the score from Score(). Sparse - most objects add nothing here.
+	struct PendingDeclutter { uint32_t outputIndex; LatpLon point; int32_t score; };
+	std::vector<PendingDeclutter> declutterOutputs;
 	std::vector<std::string> outputKeys;
 	const PbfReader::Relation* currentRelation;
 	const boost::container::flat_map<std::string, std::string>* currentPostScanTags; // for postScan only

include/shared_data.h

@@ -39,6 +39,10 @@ struct LayerDef {
 	std::string indexName;
 	std::map<std::string, uint> attributeMap; // string 0, number 1, bool 2
 	bool writeTo;
+	// Decluttering: set after addLayer(), so these live at the end of the struct
+	uint declutterBelow = 0;			// zoom below which point features are thinned out by Score()
+	double declutterDistance = 40;		// how far apart (in 256px screen pixels) they should be kept
+	double declutterThreshold = 0;		// score needed at the layer's minzoom (halves at each zoom)
 
 	const bool useColumn(std::string &col) {
 		return allSourceColumns || (std::find(sourceColumns.begin(), sourceColumns.end(), col) != sourceColumns.end());

include/shp_mem_tiles.h

@@ -3,13 +3,14 @@
 #define _SHP_MEM_TILES
 
 #include "tile_data.h"
+#include "declutter.h"
 
 extern bool verbose;
 
 class ShpMemTiles : public TileDataSource
 {
 public:
-	ShpMemTiles(size_t threadNum, uint indexZoom);
+	ShpMemTiles(size_t threadNum, uint indexZoom, class Declutter& declutter);
 
 	std::string name() const override { return "shp"; }
 
@@ -25,6 +26,7 @@ class ShpMemTiles : public TileDataSource
 		bool hasName,
 		const std::string& name,
 		uint minzoom,
+		int32_t score,
 		AttributeIndex attrIdx
 	);
 
@@ -62,6 +64,7 @@ class ShpMemTiles : public TileDataSource
 	}
 
 private:
+	class Declutter& declutter;
 	std::vector<OutputObject> indexedGeometries;				// prepared boost::geometry objects (from shapefiles)
 	std::map<uint, std::string> indexedGeometryNames;			//  | optional names for each one
 	std::map<std::string, RTree> indices;			// Spatial indices, boost::geometry::index objects for shapefile indices

include/shp_processor.h

@@ -41,11 +41,11 @@ class ShpProcessor {
 	AttributeIndex readShapefileAttributes(DBFHandle dbf, int recordNum,
 	                                       std::unordered_map<int,std::string> &columnMap,
 	                                       std::unordered_map<int,int> &columnTypeMap,
-	                                       LayerDef &layer, uint &minzoom);
+	                                       LayerDef &layer, uint &minzoom, int32_t &score);
 
 	// Process an individual shapefile record
 	void processShapeGeometry(SHPObject* shape, AttributeIndex attrIdx, 
-	                          const LayerDef &layer, uint layerNum, bool hasName, const std::string &name);
+	                          const LayerDef &layer, uint layerNum, bool hasName, const std::string &name, int32_t score);
 };
 
 #endif //_SHP_PROCESSOR_H

resources/config-schema.json

@@ -59,6 +59,9 @@
 					"simplify_ratio": { "type": "number" },
 					"filter_below": { "type": "integer", "minimum": 0 },
 					"filter_area": { "type": "number" },
+					"declutter_below": { "type": "integer", "minimum": 0 },
+					"declutter_distance": { "type": "number" },
+					"declutter_threshold": { "type": "number" },
 					"feature_limit": { "type": "integer", "minimum": 0 },
 					"feature_limit_below": { "type": "integer", "minimum": 0 },
 					"combine_points": { "type": "boolean" },

src/declutter.cpp

@@ -0,0 +1,110 @@
+#include "declutter.h"
+
+#include <algorithm>
+#include <iostream>
+#include <boost/geometry/index/rtree.hpp>
+
+#include "shared_data.h"
+#include "tile_data.h"
+
+namespace bgi = boost::geometry::index;
+
+// Features are ranked in world units (0-1 across the map in each direction), so a
+// separation of n screen pixels at zoom z is simply n/(256<<z).
+typedef boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> WorldPoint;
+
+void Declutter::configure(const std::vector<LayerDef>& layers) {
+	layerEnabled.assign(layers.size(), false);
+	for (size_t i = 0; i < layers.size(); i++) {
+		if (layers[i].declutterBelow == 0) continue;
+		layerEnabled[i] = true;
+		anyLayers = true;
+	}
+	if (!anyLayers) return;
+	entries.resize(layers.size());
+	mutexes = std::vector<std::mutex>(layers.size());
+}
+
+void Declutter::add(const OutputObject& oo, LatpLon point, uint64_t id, int32_t score, bool fromShapefile) {
+	std::lock_guard<std::mutex> lock(mutexes[oo.layer]);
+	entries[oo.layer].push_back({ oo, point, id, score, fromShapefile });
+}
+
+// Work up through the zoom levels, giving each feature the lowest zoom at which it both
+// clears the score threshold and isn't crowded out by a feature already placed. The
+// threshold halves at each zoom, and the separation is constant in screen terms, so each
+// zoom admits both lower-scoring and more tightly packed features than the one before.
+static void assignMinZooms(const LayerDef& layer, std::vector<DeclutterEntry>& list) {
+	// Highest score first; ties broken on id and position so the ranking doesn't depend on
+	// the order features happened to be read in
+	std::sort(list.begin(), list.end(), [](const DeclutterEntry& a, const DeclutterEntry& b) {
+		if (a.score != b.score) return a.score > b.score;
+		if (a.id != b.id) return a.id < b.id;
+		if (a.point.latp != b.point.latp) return a.point.latp < b.point.latp;
+		return a.point.lon < b.point.lon;
+	});
+
+	bgi::rtree<WorldPoint, bgi::quadratic<16>> placed;
+	std::vector<bool> done(list.size(), false);
+	size_t remaining = list.size();
+
+	double threshold = layer.declutterThreshold;
+	for (uint z = layer.minzoom; z < layer.declutterBelow && remaining > 0; z++, threshold /= 2) {
+		const double separation = layer.declutterDistance / (256.0 * (1u << z));
+
+		for (size_t i = 0; i < list.size(); i++) {
+			DeclutterEntry& e = list[i];
+			if (e.score < threshold) break;		// sorted by descending score, so nothing further qualifies
+			if (done[i] || e.oo.minZoom > z) continue;
+
+			// Tweak for the Cheltenham case (two significant cities next to each other) - don't let it be repeatedly pushed out by Gloucester
+			double allowed = separation;
+			if (threshold > 0 && e.score > threshold * 5) allowed = separation / (e.score / threshold / 3.0);
+
+			const double x = lon2tilexf(e.point.lon / 10000000.0, 0);
+			const double y = latp2tileyf(e.point.latp / 10000000.0, 0);
+			const WorldPoint p(x, y);
+
+			// Anything within `allowed` of p is inside this box, so we only have to measure
+			// the handful of features the box picks up
+			const boost::geometry::model::box<WorldPoint> around(
+				WorldPoint(x - allowed, y - allowed), WorldPoint(x + allowed, y + allowed));
+			bool crowded = false;
+			for (auto it = placed.qbegin(bgi::intersects(around)); it != placed.qend() && !crowded; ++it)
+				crowded = boost::geometry::distance(p, *it) < allowed;
+			if (crowded) continue;
+
+			e.oo.setMinZoom(z);
+			done[i] = true;
+			remaining--;
+			placed.insert(p);
+		}
+	}
+
+	// Anything that never won a place appears from declutter_below upwards
+	for (size_t i = 0; i < list.size(); i++)
+		if (!done[i] && list[i].oo.minZoom < layer.declutterBelow)
+			list[i].oo.setMinZoom(layer.declutterBelow);
+}
+
+void Declutter::apply(const std::vector<LayerDef>& layers, TileDataSource& osmSource, TileDataSource& shpSource) {
+	if (!anyLayers) return;
+
+	for (size_t i = 0; i < entries.size(); i++) {
+		if (entries[i].empty()) continue;
+		std::cout << "Decluttering " << entries[i].size() << " features in layer " << layers[i].name << ":" << std::flush;
+		assignMinZooms(layers[i], entries[i]);
+
+		std::vector<size_t> perZoom(16, 0);	// minZoom is a 4-bit field
+		for (const auto& e : entries[i]) perZoom[e.oo.minZoom]++;
+		for (size_t z = 0; z < perZoom.size(); z++)
+			if (perZoom[z] > 0) std::cout << " z" << z << ":" << perZoom[z];
+		std::cout << std::endl;
+
+		for (const auto& e : entries[i]) {
+			TileDataSource& source = e.fromShapefile ? shpSource : osmSource;
+			source.addObjectToSmallIndex(latpLon2index(e.point, source.getIndexZoom()), e.oo, e.id);
+		}
+		std::vector<DeclutterEntry>().swap(entries[i]);
+	}
+}