Add per-layer invalid_polygon_repair_scope option to repair clipped-and-quantised polygons

docs/CONFIGURATION.md

@@ -78,6 +78,7 @@ You can add optional parameters to layers:
 * `simplify_length` - how much to simplify features (in kilometers) on the zoom level `simplify_below-1`, preceding `simplify_level`
 * `simplify_ratio` - (optional: the default value is 2.0) the actual simplify level will be `simplify_level * pow(simplify_ratio, (simplify_below-1) - <current zoom>)`
 * `simplify_algorithm` - which simplification algorithm to use (defaults to Douglas-Peucker; you can also specify `"visvalingam"`, which can be better for landuse and similar polygons, or `"buildings"`, which preserves rectilinear shapes)
+* `invalid_polygon_repair_scope` - which invalid output polygons are repaired before being written. `"simplified_only"` (default) repairs just those that simplification may have broken; `"all"` is a superset and additionally repairs polygons that were merely clipped and quantised onto the integer tile grid. Use `"all"` for layers built from large clipped polygons such as coastline/ocean shapefiles, whose tiles are otherwise written out invalid and can be mis-rendered by strict renderers. For those, a repair is only accepted if it preserves the polygon's area, so it can never fill a hole and make an island disappear.
 * `filter_below` - filter areas by minimum size below this zoom level
 * `filter_area` - minimum size (in square degrees of longitude) for the zoom level `filter_below-1`
 * `feature_limit` - restrict the number of features written to each tile

include/geom.h

@@ -82,7 +82,13 @@ void make_valid(MultiPolygon &mp);
 // make_valid first (preserves area), then a zero-width buffer as a last
 // resort. Returns true if mp is valid afterwards; on failure mp is left as the
 // best-effort input so callers never regress.
-bool repair_multi_polygon(MultiPolygon &mp);
+//
+// `strictArea` additionally rejects any repair that GROWS a polygon, i.e. that
+// filled a hole. Use it for geometry that was not simplified: such geometry is
+// only invalid because of the integer quantisation onto the tile grid, so a
+// correct repair must not change its area. Simplified geometry legitimately
+// changes area and must keep the lenient bound.
+bool repair_multi_polygon(MultiPolygon &mp, bool strictArea = false);
 
 void union_many(std::vector<MultiPolygon> &mps);
 

include/shared_data.h

@@ -43,6 +43,8 @@ struct LayerDef {
 	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)
+	// Geometry repair: also set after addLayer()
+	uint repairScope = 0;				// REPAIR_SIMPLIFIED_ONLY (default) or REPAIR_ALL
 
 	const bool useColumn(std::string &col) {
 		return allSourceColumns || (std::find(sourceColumns.begin(), sourceColumns.end(), col) != sourceColumns.end());
@@ -51,6 +53,16 @@ struct LayerDef {
 	static const uint DOUGLAS_PEUCKER = 0;
 	static const uint VISVALINGAM = 1;
 	static const uint BUILDINGS = 2;
+
+	// Which invalid output polygons are repaired before being written to a tile.
+	// REPAIR_SIMPLIFIED_ONLY (default) repairs only polygons that simplification
+	// may have broken - the historic behaviour. REPAIR_ALL is a superset: it also
+	// repairs polygons that were merely clipped and quantised onto the integer
+	// tile grid (e.g. ocean shapefiles). For those, a repair must not GROW the
+	// polygon, because growing means a hole was filled and an island would
+	// silently disappear.
+	static const uint REPAIR_SIMPLIFIED_ONLY = 0;
+	static const uint REPAIR_ALL = 1;
 };
 
 ///\brief Defines layers used in map rendering

resources/config-schema.json

@@ -69,6 +69,7 @@
 					"combine_polygons_below": { "type": "integer", "minimum": 0 },
 					"z_order_ascending": { "type": "boolean" },
 					"simplify_algorithm": { "type": "string" },
+					"invalid_polygon_repair_scope": { "type": "string", "enum": [ "simplified_only", "all" ] },
 					"source": { "type": "string" },
 					"source_columns": {
 						"oneOf": [

src/geom.cpp

@@ -7,6 +7,7 @@
 #include <boost/geometry/strategies/buffer.hpp>
 
 #include "geometry/correct.hpp"
+#include <limits>
 
 #include <boost/range/adaptor/transformed.hpp>
 #include <boost/range/irange.hpp>
@@ -148,14 +149,27 @@ void make_valid(MultiPolygon &mp)
 
 // Repair a single (possibly invalid) polygon in an area-preserving way and
 // append the resulting valid polygon(s) to `out`. Returns true on success.
-// `minArea` is the lower bound on the repaired area we are willing to accept.
-static bool repair_one_polygon(const Polygon &p, double minArea, MultiPolygon &out)
+// A repair is only accepted if the resulting area stays within
+// [minArea, maxArea]:
+//  - the LOWER bound catches a collapse (dissolve/buffer dropping most of the
+//    covered area on huge inputs);
+//  - the UPPER bound catches a FILLED HOLE, which is what used to make an ocean
+//    polygon swallow a whole island. Measured on real data, the separation is
+//    large: genuine repairs of quantised ocean polygons change the area by
+//    <=0.01 %, whereas filling one island hole grew it by 32 %.
+static bool repair_one_polygon(const Polygon &p, double minArea, double maxArea, MultiPolygon &out)
 {
+	auto acceptable = [&](const MultiPolygon &candidate) {
+		if (!geom::is_valid(candidate)) return false;
+		const double a = std::abs(geom::area(candidate));
+		return a >= minArea && a <= maxArea;
+	};
+
 	// 1) Dissolve (resolves self-intersections of this single polygon).
 	try {
 		MultiPolygon fixed;
 		geometry::correct(p, fixed, 1E-12);
-		if (geom::is_valid(fixed) && std::abs(geom::area(fixed)) >= minArea) {
+		if (acceptable(fixed)) {
 			for (auto &fp : fixed) out.push_back(std::move(fp));
 			return true;
 		}
@@ -174,7 +188,7 @@ static bool repair_one_polygon(const Polygon &p, double minArea, MultiPolygon &o
 
 		geom::buffer(p, buffered, distanceStrategy, sideStrategy, joinStrategy, endStrategy, pointStrategy);
 		geom::correct(buffered);
-		if (geom::is_valid(buffered) && std::abs(geom::area(buffered)) >= minArea) {
+		if (acceptable(buffered)) {
 			for (auto &bp : buffered) out.push_back(std::move(bp));
 			return true;
 		}
@@ -185,7 +199,7 @@ static bool repair_one_polygon(const Polygon &p, double minArea, MultiPolygon &o
 	return false;
 }
 
-bool repair_multi_polygon(MultiPolygon &mp)
+bool repair_multi_polygon(MultiPolygon &mp, bool strictArea)
 {
 	if (geom::is_valid(mp)) return true;
 
@@ -207,13 +221,19 @@ bool repair_multi_polygon(MultiPolygon &mp)
 			out.push_back(p);
 			continue;
 		}
-		// Lenient threshold: resolving a self-intersection legitimately changes a
-		// single polygon's (shoelace) area, so anything down to half the original
-		// is accepted. Per-polygon repair cannot trigger the cross-polygon union
-		// that previously caused the catastrophic ~99% collapse, so this only
-		// rejects a genuine local collapse.
-		const double minArea = 0.5 * std::abs(geom::area(p));
-		if (!repair_one_polygon(p, minArea, out)) {
+		// Lenient lower bound: resolving a self-intersection legitimately changes
+		// a single polygon's (shoelace) area, so anything down to half the
+		// original is accepted. Per-polygon repair cannot trigger the
+		// cross-polygon union that previously caused the catastrophic ~99%
+		// collapse, so this only rejects a genuine local collapse.
+		// The upper bound is only applied in strict mode, i.e. for geometry that
+		// was not simplified. Applying it to simplified geometry would be wrong:
+		// it rejects legitimate repairs there and measurably made things worse
+		// (5 -> 20 invalid features in one tile).
+		const double origArea = std::abs(geom::area(p));
+		const double maxArea = strictArea ? 1.01 * origArea
+		                                  : std::numeric_limits<double>::infinity();
+		if (!repair_one_polygon(p, 0.5 * origArea, maxArea, out)) {
 			out.push_back(p);
 			allValid = false;
 		}

src/shared_data.cpp

@@ -349,6 +349,20 @@ void Config::readConfig(rapidjson::Document &jsonConfig, bool &hasClippingBox, B
 				source, sourceColumns, allSourceColumns, indexed, indexName,
 				writeTo);
 
+		// Which invalid polygons to repair: "simplified_only" (default) repairs
+		// what simplification may have broken; "all" additionally repairs
+		// clipped-and-quantised geometry such as ocean shapefiles, which is
+		// otherwise written out invalid.
+		if (it->value.HasMember("invalid_polygon_repair_scope")) {
+			const string scope = it->value["invalid_polygon_repair_scope"].GetString();
+			if (scope == "all") {
+				layers.layers[layerNum].repairScope = LayerDef::REPAIR_ALL;
+			} else if (scope != "simplified_only") {
+				cerr << "Unknown invalid_polygon_repair_scope \"" << scope << "\" in layer " << layerName
+				     << "; expected \"simplified_only\" or \"all\"" << endl;
+			}
+		}
+
 		// Decluttering (thinning out point features by the score their profile gives them)
 		if (it->value.HasMember("declutter_below")) {
 			LayerDef &layerDef = layers.layers[layerNum];

src/tile_worker.cpp

@@ -221,6 +221,7 @@ void writeMultiPolygon(
 	unsigned zoom,
 	double simplifyLevel,
 	unsigned simplifyAlgo,
+	unsigned repairScope,
 	const MultiPolygon& mp
 ) {
 	bbox.scaleGeometry(scaledMultiPolygon, mp);
@@ -243,8 +244,19 @@ void writeMultiPolygon(
 
 	geom::correct(current);
 
+	// Simplification can turn a valid input into a self-intersecting/spiky one;
+	// such polygons are silently dropped by many renderers (missing features), so
+	// those are always repaired. Geometry that was NOT simplified is only invalid
+	// because of the integer quantisation onto the tile grid; repairing that as
+	// well is opt-in per layer ("invalid_polygon_repair_scope": "all"), and uses
+	// the strict area guard so a repair can never fill a hole and hide an island.
+	const bool mayRepair = (simplifyLevel > 0 || repairScope == LayerDef::REPAIR_ALL);
+
+	// is_valid() runs full self-intersection detection, so only evaluate it when
+	// the answer can change what we write: for a layer that is neither simplified
+	// nor opted into repair the block below is a no-op.
 	geom::validity_failure_type failure;
-	if (!geom::is_valid(current, failure)) {
+	if ((mayRepair || verbose) && !geom::is_valid(current, failure)) {
 		if (verbose) {
 			cout << "output multipolygon has " << boost_validity_error(failure) << endl;
 
@@ -253,12 +265,9 @@ void writeMultiPolygon(
 			else
 				cout << "input multipolygon valid" << endl;
 		}
-
-		if (simplifyLevel > 0) {
-			// Simplification can turn a valid input into a self-intersecting/spiky
-			// one; such polygons are silently dropped by many renderers (missing
-			// features). Repair (dissolve, then zero-width buffer) before writing.
-			bool repaired = repair_multi_polygon(current);
+
+		if (mayRepair) {
+			bool repaired = repair_multi_polygon(current, /*strictArea=*/ simplifyLevel == 0);
 
 			if (geom::is_empty(current))
 				return;
@@ -308,6 +317,7 @@ void ProcessObjects(
 	class SharedData& sharedData,
 	double simplifyLevel,
 	unsigned simplifyAlgo,
+	unsigned repairScope,
 	double filterArea,
 	bool combinePoints,
 	bool combineLines,
@@ -412,7 +422,7 @@ void ProcessObjects(
 			if (oo.oo.geomType == LINESTRING_ || oo.oo.geomType == MULTILINESTRING_)
 				writeMultiLinestring(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get<MultiLinestring>(g));
 			else if (oo.oo.geomType == POLYGON_)
-				writeMultiPolygon(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, boost::get<MultiPolygon>(g));
+				writeMultiPolygon(attributeStore, sharedData, vtLayer, bbox, oo, zoom, simplifyLevel, simplifyAlgo, repairScope, boost::get<MultiPolygon>(g));
 		}
 	}
 }
@@ -495,7 +505,7 @@ void ProcessLayer(
 			if (ld.featureLimit>0 && end-ooListSameLayer.first>ld.featureLimit && zoom<ld.featureLimitBelow) end = ooListSameLayer.first+ld.featureLimit;
 			ProcessObjects(sources[i], attributeStore, 
 				ooListSameLayer.first, end, sharedData, 
-				simplifyLevel, ld.simplifyAlgo,
+				simplifyLevel, ld.simplifyAlgo, ld.repairScope,
 				filterArea, ld.combinePoints, zoom < ld.combineLinesBelow, zoom < ld.combinePolygonsBelow, zoom, bbox, vtLayer);
 		}
 	}