diff --git a/CHANGELOG.md b/CHANGELOG.md index 55d604d..14f1867 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,15 @@ This project adheres to [Semantic Versioning](https://semver.org/). ### Fixed +- Islands mapped with the coastline going the wrong way can end up as holes + of a land polygon they are not inside of ("Hole lies outside shell"). 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..43a07dc 100644 --- a/src/osmcoastline.cpp +++ b/src/osmcoastline.cpp @@ -40,6 +40,7 @@ #include #include +#include #include #include #include @@ -71,10 +72,124 @@ const unsigned int max_warnings = 500; /* ================================================== */ +/** + * Is the specified ring outside the specified (exterior) ring? Only the + * first point of the inner ring is checked, so this will not detect rings + * that are only partly outside. Those cases are reported as intersections + * elsewhere. This uses the simple point-in-ring test instead of a GEOS + * operation, because it has to work on invalid geometries, too. + */ +bool ring_is_outside(const OGRLinearRing* exterior_ring, const OGRLinearRing* ring) { + OGRPoint point; + ring->getPoint(0, &point); + return !exterior_ring->isPointInRing(&point, TRUE); +} + +/** + * organizePolygons() decides which rings are holes based on their direction + * and only checks the bounding box when assigning a hole to a polygon. So an + * island whose coastline was mapped the wrong way round can end up as a hole + * of a polygon it isn't inside of. GEOS then reports "Hole lies outside + * shell" and the whole polygon (possibly a whole continent) is invalid. + * + * Take those rings out of the polygon, report them, and turn them around + * into land polygons of their own. + */ +std::unique_ptr fix_holes_outside_shell(std::unique_ptr polygon, + polygon_vector_type* polygons, + OutputDatabase& output, + unsigned int* turned_around) { + const OGRLinearRing* exterior_ring = polygon->getExteriorRing(); + assert(exterior_ring); + + std::vector is_outside; + is_outside.reserve(polygon->getNumInteriorRings()); + unsigned int num_outside = 0; + for (int i = 0; i < polygon->getNumInteriorRings(); ++i) { + const OGRLinearRing* interior_ring = polygon->getInteriorRing(i); + assert(interior_ring); + is_outside.push_back(ring_is_outside(exterior_ring, interior_ring)); + if (is_outside.back()) { + ++num_outside; + } + } + + if (num_outside == 0) { + return polygon; + } + + auto fixed_polygon = std::make_unique(); + fixed_polygon->addRingDirectly(exterior_ring->clone()); + + for (int i = 0; i < polygon->getNumInteriorRings(); ++i) { + const OGRLinearRing* interior_ring = polygon->getInteriorRing(i); + assert(interior_ring); + + if (!is_outside[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_holes_outside_shell(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 +202,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 +212,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 +236,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 +449,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 +458,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..9abe136 100755 --- a/test/t/valid-inland-sea-with-island.sh +++ b/test/t/valid-inland-sea-with-island.sh @@ -37,21 +37,24 @@ OSM RC=$? set -e -test $RC -eq 2 +test $RC -eq 1 if [ "$SRID" = "4326" ]; then - grep 'Found 3 rings in input data.$' "$LOG" + grep 'Found 2 rings in input data.$' "$LOG" grep '^There were 3 warnings.$' "$LOG" check_count error_lines 3; 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; +# 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. The invalid +# polygon is repaired, which drops the island. +check_count land_polygons 1; check_count error_points 0; #-----------------------------------------------------------------------------