From 79a14a72e9051ee2db93b31865e08e2064b115fa Mon Sep 17 00:00:00 2001 From: Lee Kindness Date: Tue, 18 Aug 2026 20:32:24 +0100 Subject: [PATCH] Fix "Hole lies outside shell" for reversed islands organizePolygons() decides which rings are holes based on their direction only and uses just a bounding box test when assigning a hole to a polygon. A ring whose coastline was mapped the wrong way round can therefore end up as a hole of a polygon it is not inside of, or as a hole nested inside another hole. GEOS then reports "Hole lies outside shell" or "Holes are nested" and the whole polygon, possibly a whole continent, is invalid and ends up in the error_lines table. The existing direction fix runs after the polygons have been created and only looks at exterior rings, so it never caught these cases. Those rings are not holes at all but land mapped the wrong way round. They are now taken out of the polygon, turned around into land polygons of their own and reported in the error_lines table with the error "direction". While at it, invalid polygons that are not part of a multipolygon are now repaired in the same way as those that are, instead of being dropped. Fixes #41 Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 10 + src/osmcoastline.cpp | 193 +++++++++++++++--- ...nvalid-direction-island-in-bbox-of-many.sh | 65 ++++++ test/t/invalid-direction-island-in-bbox.sh | 63 ++++++ test/t/valid-inland-sea-with-island.sh | 30 ++- 5 files changed, 326 insertions(+), 35 deletions(-) create mode 100755 test/t/invalid-direction-island-in-bbox-of-many.sh create mode 100755 test/t/invalid-direction-island-in-bbox.sh diff --git a/CHANGELOG.md b/CHANGELOG.md index 55d604d..913fe4f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,16 @@ This project adheres to [Semantic Versioning](https://semver.org/). ### Fixed +- Coastlines mapped the wrong way round can end up as holes of a land polygon + they are not inside of ("Hole lies outside shell") or as holes inside + another hole ("Holes are nested"). This made the (possibly very large) land + polygon invalid. Those rings are now turned around into land polygons of + their own and reported in the `error_lines` table with the error + `direction` (#41). +- Land polygons that are a single (invalid) polygon are now repaired in the + same way as polygons that are part of a multipolygon instead of being + dropped. + ## [2.5.0] - 2026-01-18 diff --git a/src/osmcoastline.cpp b/src/osmcoastline.cpp index 9568f2a..545dd0c 100644 --- a/src/osmcoastline.cpp +++ b/src/osmcoastline.cpp @@ -38,8 +38,10 @@ #include #include +#include #include #include +#include #include #include #include @@ -71,10 +73,166 @@ const unsigned int max_warnings = 500; /* ================================================== */ +/** + * Is the point inside the ring? The envelope of the ring has to be given, + * because calculating it is expensive and we need it several times. + * + * This uses the simple point-in-ring test instead of a GEOS operation, + * because it has to work on invalid geometries, too. + */ +[[nodiscard]] bool point_in_ring(const OGRLinearRing* ring, const OGREnvelope& envelope, const OGRPoint& point) { + if (point.getX() < envelope.MinX || point.getX() > envelope.MaxX || + point.getY() < envelope.MinY || point.getY() > envelope.MaxY) { + return false; + } + + return ring->isPointInRing(&point, false); +} + +/** + * Find the interior rings of the polygon that can not be holes of that + * polygon, because they are outside the exterior ring or inside another + * hole. Only the first point of each interior ring is looked at, so rings + * that are only partly outside are not detected. Those cases show up as + * intersections and are reported elsewhere. + */ +[[nodiscard]] std::vector find_misplaced_holes(const OGRPolygon& polygon) { + const auto num_rings = static_cast(polygon.getNumInteriorRings()); + + std::vector envelopes(num_rings); + std::vector first_points(num_rings); + for (std::size_t i = 0; i < num_rings; ++i) { + const OGRLinearRing* ring = polygon.getInteriorRing(static_cast(i)); + assert(ring); + ring->getEnvelope(&envelopes[i]); + ring->getPoint(0, &first_points[i]); + } + + const OGRLinearRing* exterior_ring = polygon.getExteriorRing(); + assert(exterior_ring); + OGREnvelope exterior_envelope; + exterior_ring->getEnvelope(&exterior_envelope); + + // rings outside the exterior ring ("Hole lies outside shell") + std::vector outside(num_rings, false); + for (std::size_t i = 0; i < num_rings; ++i) { + outside[i] = !point_in_ring(exterior_ring, exterior_envelope, first_points[i]); + } + + // rings inside another hole ("Holes are nested"). Rings that are outside + // the exterior ring are no holes at all, so they are not considered as + // the containing ring here. + std::vector misplaced{outside}; + for (std::size_t i = 0; i < num_rings; ++i) { + if (misplaced[i]) { + continue; + } + for (std::size_t j = 0; j < num_rings; ++j) { + if (i != j && !outside[j] && envelopes[j].Contains(envelopes[i]) && + point_in_ring(polygon.getInteriorRing(static_cast(j)), envelopes[j], first_points[i])) { + misplaced[i] = true; + break; + } + } + } + + return misplaced; +} + +/** + * organizePolygons() decides which rings are holes based on their direction + * and only checks the bounding box when assigning a hole to a polygon. So a + * ring whose coastline was mapped the wrong way round can end up as a hole + * of a polygon it isn't inside of or as a hole inside another hole. GEOS + * then reports "Hole lies outside shell" or "Holes are nested" and the whole + * polygon (possibly a whole continent) is invalid. + * + * Those rings are not holes but land mapped the wrong way round. Take them + * out of the polygon, report them, and turn them around into land polygons + * of their own. + */ +[[nodiscard]] std::unique_ptr fix_misplaced_holes(std::unique_ptr polygon, + polygon_vector_type* polygons, + OutputDatabase& output, + unsigned int* turned_around) { + const std::vector misplaced = find_misplaced_holes(*polygon); + + if (std::find(misplaced.cbegin(), misplaced.cend(), true) == misplaced.cend()) { + return polygon; + } + + auto fixed_polygon = std::make_unique(); + fixed_polygon->addRingDirectly(polygon->getExteriorRing()->clone()); + + for (int i = 0; i < polygon->getNumInteriorRings(); ++i) { + const OGRLinearRing* interior_ring = polygon->getInteriorRing(i); + assert(interior_ring); + + if (!misplaced[static_cast(i)]) { + fixed_polygon->addRingDirectly(interior_ring->clone()); + continue; + } + + auto ring = std::unique_ptr(interior_ring->clone()); + ring->reversePoints(); + + auto ls = std::unique_ptr(OGRGeometryFactory::forceToLineString(ring->clone())->toLineString()); + output.add_error_line(std::move(ls), "direction"); + + auto island = std::make_unique(); + island->addRingDirectly(ring.release()); + island->assignSpatialReference(srs.wgs84()); + polygons->push_back(std::move(island)); + + ++(*turned_around); + } + + fixed_polygon->assignSpatialReference(srs.wgs84()); + return fixed_polygon; +} + +/** + * Add the polygon to the list of polygons, trying to fix it if it isn't + * valid. + */ +void add_polygon_to(polygon_vector_type* polygons, + std::unique_ptr polygon, + OutputDatabase& output, + unsigned int* warnings, unsigned int* errors, + unsigned int* turned_around) { + if (polygon->IsValid()) { + polygons->push_back(std::move(polygon)); + return; + } + + if (polygon->getNumInteriorRings() > 0) { + polygon = fix_misplaced_holes(std::move(polygon), polygons, output, turned_around); + if (polygon->IsValid()) { + polygons->push_back(std::move(polygon)); + return; + } + } + + auto* ring = polygon->getExteriorRing()->clone(); + auto ls = std::unique_ptr(OGRGeometryFactory::forceToLineString(ring)->toLineString()); + output.add_error_line(std::move(ls), "invalid"); + + std::unique_ptr buf0{polygon->Buffer(0)}; + if (buf0 && buf0->getGeometryType() == wkbPolygon && buf0->IsValid()) { + buf0->assignSpatialReference(srs.wgs84()); + polygons->push_back(static_cast_unique_ptr(std::move(buf0))); + (*warnings)++; + } else { + std::cerr << "Ignoring invalid polygon geometry.\n"; + (*errors)++; + } +} + void add_polygons_in_multi_to(polygon_vector_type *polygons, std::unique_ptr mega_geometry, OutputDatabase& output, - unsigned int* warnings, unsigned int* errors) { + unsigned int* warnings, unsigned int* errors, + unsigned int* turned_around) { // This isn't an owning pointer on purpose. We are going to "steal" parts // of the geometry a few lines below but only mark them as unowned farther // below when we are calling removeGeometry() on it. If this was an @@ -87,22 +245,7 @@ void add_polygons_in_multi_to(polygon_vector_type *polygons, assert(geom); assert(geom->getGeometryType() == wkbPolygon); std::unique_ptr p{static_cast(geom)}; - if (p->IsValid()) { - polygons->push_back(std::move(p)); - } else { - auto* ring = p->getExteriorRing()->clone(); - auto ls = std::unique_ptr(OGRGeometryFactory::forceToLineString(ring)->toLineString()); - output.add_error_line(std::move(ls), "invalid"); - std::unique_ptr buf0{p->Buffer(0)}; - if (buf0 && buf0->getGeometryType() == wkbPolygon && buf0->IsValid()) { - buf0->assignSpatialReference(srs.wgs84()); - polygons->push_back(static_cast_unique_ptr(std::move(buf0))); - (*warnings)++; - } else { - std::cerr << "Ignoring invalid polygon geometry.\n"; - (*errors)++; - } - } + add_polygon_to(polygons, std::move(p), output, warnings, errors, turned_around); } mega_multipolygon->removeGeometry(-1, FALSE); @@ -112,7 +255,7 @@ void add_polygons_in_multi_to(polygon_vector_type *polygons, /** * This function assembles all the coastline rings into one huge multipolygon. */ -polygon_vector_type create_polygons(CoastlineRingCollection& coastline_rings, OutputDatabase& output, unsigned int* warnings, unsigned int* errors) { +polygon_vector_type create_polygons(CoastlineRingCollection& coastline_rings, OutputDatabase& output, unsigned int* warnings, unsigned int* errors, unsigned int* turned_around) { std::vector all_polygons = coastline_rings.add_polygons_to_vector(); if (all_polygons.empty()) { @@ -136,16 +279,11 @@ polygon_vector_type create_polygons(CoastlineRingCollection& coastline_rings, Ou polygon_vector_type polygons; if (mega_geometry->getGeometryType() == wkbPolygon) { - if (mega_geometry->IsValid()) { - polygons.push_back(static_cast_unique_ptr(std::move(mega_geometry))); - } else { - std::cerr << "Ignoring invalid polygon geometry.\n"; - (*errors)++; - } + add_polygon_to(&polygons, static_cast_unique_ptr(std::move(mega_geometry)), output, warnings, errors, turned_around); } else if (mega_geometry->getGeometryType() != wkbMultiPolygon) { throw std::runtime_error{"mega geometry isn't a (multi)polygon. Something is very wrong!"}; } else { - add_polygons_in_multi_to(&polygons, std::move(mega_geometry), output, warnings, errors); + add_polygons_in_multi_to(&polygons, std::move(mega_geometry), output, warnings, errors, turned_around); } return polygons; @@ -354,7 +492,8 @@ int main(int argc, char *argv[]) { if (options.output_polygons != output_polygon_type::none || options.output_lines) { try { vout << "Create polygons...\n"; - CoastlinePolygons coastline_polygons{create_polygons(coastline_rings, *output_database, &warnings, &errors), \ + unsigned int turned_around = 0; + CoastlinePolygons coastline_polygons{create_polygons(coastline_rings, *output_database, &warnings, &errors, &turned_around), \ *output_database, \ options.bbox_overlap, \ options.max_points_in_polygon}; @@ -362,7 +501,7 @@ int main(int argc, char *argv[]) { stats.land_polygons_before_split = coastline_polygons.num_polygons(); vout << "Fixing coastlines going the wrong way...\n"; - stats.rings_turned_around = coastline_polygons.fix_direction(); + stats.rings_turned_around = turned_around + coastline_polygons.fix_direction(); vout << " Turned " << stats.rings_turned_around << " polygons around.\n"; warnings += stats.rings_turned_around; diff --git a/test/t/invalid-direction-island-in-bbox-of-many.sh b/test/t/invalid-direction-island-in-bbox-of-many.sh new file mode 100755 index 0000000..330df6f --- /dev/null +++ b/test/t/invalid-direction-island-in-bbox-of-many.sh @@ -0,0 +1,65 @@ +#!/bin/sh +#----------------------------------------------------------------------------- +# +# Same as invalid-direction-island-in-bbox, but with a second, correctly +# mapped island. This means organizePolygons() returns a multipolygon which +# is handled by a different code path. +# See https://github.com/osmcode/osmcoastline/issues/41 +# +#----------------------------------------------------------------------------- + +# shellcheck source=test/init.sh +. "$1/test/init.sh" + +set -x + +#----------------------------------------------------------------------------- + +"$BIN_DIR/src/nodegrid2opl" << 'NODES' >"$INPUT" + + 0---------------1 + | | + | | a---b + | 3--------2 | | + | | d---c + | | 6---7 + | | | | + | | 9---8 + | | + 5------4 + +NODES + +cat <<'OSM' >>"$INPUT" +w200 v1 Tnatural=coastline Nn100,n105,n104,n103,n102,n101,n100 +w201 v1 Tnatural=coastline Nn106,n107,n108,n109,n106 +w202 v1 Tnatural=coastline Nn110,n113,n112,n111,n110 +OSM + +#----------------------------------------------------------------------------- + +"$OSMC" --verbose --overwrite --srs="$SRID" --output-database="$DB" "$INPUT" >"$LOG" 2>&1 +RC=$? +set -e + +test $RC -eq 1 + +grep 'Turned 1 polygons around.$' "$LOG" + +grep '^There were 1 warnings.$' "$LOG" +grep '^There were 0 errors.$' "$LOG" + +check_count land_polygons 3; +check_count error_points 0; +check_count error_lines 1; + +echo "SELECT InsertEpsgSrid(4326);" | $SQL + +# the island that was turned around +echo "SELECT AsText(Transform(geometry, 4326)) FROM land_polygons;" | $SQL \ + | grep -F 'POLYGON((1.15 1.94, 1.19 1.94, 1.19 1.92, 1.15 1.92, 1.15 1.94))' + +echo "SELECT AsText(Transform(geometry, 4326)), osm_id, error FROM error_lines;" | $SQL \ + | grep -F '|0|direction' + +#----------------------------------------------------------------------------- diff --git a/test/t/invalid-direction-island-in-bbox.sh b/test/t/invalid-direction-island-in-bbox.sh new file mode 100755 index 0000000..fd05867 --- /dev/null +++ b/test/t/invalid-direction-island-in-bbox.sh @@ -0,0 +1,63 @@ +#!/bin/sh +#----------------------------------------------------------------------------- +# +# Invalid island with coastline going the wrong direction inside the +# bounding box of (but outside) a larger land polygon. Used to result in a +# "Hole lies outside shell" error from GEOS which made the larger polygon +# invalid. See https://github.com/osmcode/osmcoastline/issues/41 +# +#----------------------------------------------------------------------------- + +# shellcheck source=test/init.sh +. "$1/test/init.sh" + +set -x + +#----------------------------------------------------------------------------- + +"$BIN_DIR/src/nodegrid2opl" << 'NODES' >"$INPUT" + + 0---------------1 + | | + | | + | 3--------2 + | | + | | 6---7 + | | | | + | | 9---8 + | | + 5------4 + +NODES + +cat <<'OSM' >>"$INPUT" +w200 v1 Tnatural=coastline Nn100,n105,n104,n103,n102,n101,n100 +w201 v1 Tnatural=coastline Nn106,n107,n108,n109,n106 +OSM + +#----------------------------------------------------------------------------- + +"$OSMC" --verbose --overwrite --srs="$SRID" --output-database="$DB" "$INPUT" >"$LOG" 2>&1 +RC=$? +set -e + +test $RC -eq 1 + +grep 'Turned 1 polygons around.$' "$LOG" + +grep '^There were 1 warnings.$' "$LOG" +grep '^There were 0 errors.$' "$LOG" + +check_count land_polygons 2; +check_count error_points 0; +check_count error_lines 1; + +echo "SELECT InsertEpsgSrid(4326);" | $SQL + +echo "SELECT AsText(Transform(geometry, 4326)) FROM land_polygons;" | $SQL \ + | grep -F 'POLYGON((1.15 1.94, 1.19 1.94, 1.19 1.92, 1.15 1.92, 1.15 1.94))' + +echo "SELECT AsText(Transform(geometry, 4326)), osm_id, error FROM error_lines;" | $SQL \ + | grep -F '|0|direction' + +#----------------------------------------------------------------------------- diff --git a/test/t/valid-inland-sea-with-island.sh b/test/t/valid-inland-sea-with-island.sh index c728468..70bf005 100755 --- a/test/t/valid-inland-sea-with-island.sh +++ b/test/t/valid-inland-sea-with-island.sh @@ -37,21 +37,35 @@ OSM RC=$? set -e -test $RC -eq 2 +test $RC -eq 1 + +# The island in the inland sea is mapped the wrong way round, so it ends up +# as a hole nested inside the hole formed by the inland sea. It is turned +# around into a land polygon of its own and reported as "direction". The +# inland sea is still reported as "questionable". +grep 'Turned 1 polygons around.$' "$LOG" if [ "$SRID" = "4326" ]; then - grep 'Found 3 rings in input data.$' "$LOG" - grep '^There were 3 warnings.$' "$LOG" - check_count error_lines 3; + grep 'Found 1 rings in input data.$' "$LOG" + grep '^There were 2 warnings.$' "$LOG" + check_count error_lines 2; else # "questionables" are not checked in 3857 - grep '^There were 0 warnings.$' "$LOG" - check_count error_lines 0; + grep '^There were 1 warnings.$' "$LOG" + check_count error_lines 1; fi -grep '^There were 1 errors.$' "$LOG" +grep '^There were 0 errors.$' "$LOG" -check_count land_polygons 0; +check_count land_polygons 2; check_count error_points 0; +echo "SELECT InsertEpsgSrid(4326);" | $SQL + +echo "SELECT AsText(Transform(geometry, 4326)) FROM land_polygons;" | $SQL \ + | grep -F 'POLYGON((1.13 1.96, 1.14 1.96, 1.14 1.95, 1.13 1.96))' + +echo "SELECT AsText(Transform(geometry, 4326)), osm_id, error FROM error_lines;" | $SQL \ + | grep -F '|0|direction' + #-----------------------------------------------------------------------------