Fix "Hole lies outside shell" for reversed islands

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
 

src/osmcoastline.cpp

@@ -38,8 +38,10 @@
 #include <ogr_core.h>
 #include <ogr_geometry.h>
 
+#include <algorithm>
 #include <cassert>
 #include <cerrno>
+#include <cstddef>
 #include <cstdlib>
 #include <cstring>
 #include <exception>
@@ -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<bool> find_misplaced_holes(const OGRPolygon& polygon) {
+    const auto num_rings = static_cast<std::size_t>(polygon.getNumInteriorRings());
+
+    std::vector<OGREnvelope> envelopes(num_rings);
+    std::vector<OGRPoint> first_points(num_rings);
+    for (std::size_t i = 0; i < num_rings; ++i) {
+        const OGRLinearRing* ring = polygon.getInteriorRing(static_cast<int>(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<bool> 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<bool> 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<int>(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<OGRPolygon> fix_misplaced_holes(std::unique_ptr<OGRPolygon> polygon,
+                                                              polygon_vector_type* polygons,
+                                                              OutputDatabase& output,
+                                                              unsigned int* turned_around) {
+    const std::vector<bool> misplaced = find_misplaced_holes(*polygon);
+
+    if (std::find(misplaced.cbegin(), misplaced.cend(), true) == misplaced.cend()) {
+        return polygon;
+    }
+
+    auto fixed_polygon = std::make_unique<OGRPolygon>();
+    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<std::size_t>(i)]) {
+            fixed_polygon->addRingDirectly(interior_ring->clone());
+            continue;
+        }
+
+        auto ring = std::unique_ptr<OGRLinearRing>(interior_ring->clone());
+        ring->reversePoints();
+
+        auto ls = std::unique_ptr<OGRLineString>(OGRGeometryFactory::forceToLineString(ring->clone())->toLineString());
+        output.add_error_line(std::move(ls), "direction");
+
+        auto island = std::make_unique<OGRPolygon>();
+        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<OGRPolygon> 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<OGRLineString>(OGRGeometryFactory::forceToLineString(ring)->toLineString());
+    output.add_error_line(std::move(ls), "invalid");
+
+    std::unique_ptr<OGRGeometry> buf0{polygon->Buffer(0)};
+    if (buf0 && buf0->getGeometryType() == wkbPolygon && buf0->IsValid()) {
+        buf0->assignSpatialReference(srs.wgs84());
+        polygons->push_back(static_cast_unique_ptr<OGRPolygon>(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<OGRGeometry> 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<OGRPolygon> p{static_cast<OGRPolygon*>(geom)};
-        if (p->IsValid()) {
-            polygons->push_back(std::move(p));
-        } else {
-            auto* ring = p->getExteriorRing()->clone();
-            auto ls = std::unique_ptr<OGRLineString>(OGRGeometryFactory::forceToLineString(ring)->toLineString());
-            output.add_error_line(std::move(ls), "invalid");
-            std::unique_ptr<OGRGeometry> buf0{p->Buffer(0)};
-            if (buf0 && buf0->getGeometryType() == wkbPolygon && buf0->IsValid()) {
-                buf0->assignSpatialReference(srs.wgs84());
-                polygons->push_back(static_cast_unique_ptr<OGRPolygon>(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<OGRGeometry*> 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<OGRPolygon>(std::move(mega_geometry)));
-        } else {
-            std::cerr << "Ignoring invalid polygon geometry.\n";
-            (*errors)++;
-        }
+        add_polygon_to(&polygons, static_cast_unique_ptr<OGRPolygon>(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,15 +492,16 @@ 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};
 
             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;
 

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'
+
+#-----------------------------------------------------------------------------

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'
+
+#-----------------------------------------------------------------------------