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