Skip to content

Commit 73c8ffc

Browse files
jotojoto
authored
Merge pull request #53 from wangi/hole-lies-outside-shell
Fix "Hole lies outside shell" for reversed islands
2 parents 2c6dd17 + 79a14a7 commit 73c8ffc

5 files changed

Lines changed: 326 additions & 35 deletions

‎CHANGELOG.md‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,16 @@ This project adheres to [Semantic Versioning](https://semver.org/).
1212

1313
### Fixed
1414

15+
- Coastlines mapped the wrong way round can end up as holes of a land polygon
16+
they are not inside of ("Hole lies outside shell") or as holes inside
17+
another hole ("Holes are nested"). This made the (possibly very large) land
18+
polygon invalid. Those rings are now turned around into land polygons of
19+
their own and reported in the `error_lines` table with the error
20+
`direction` (#41).
21+
- Land polygons that are a single (invalid) polygon are now repaired in the
22+
same way as polygons that are part of a multipolygon instead of being
23+
dropped.
24+
1525

1626
## [2.5.0] - 2026-01-18
1727

‎src/osmcoastline.cpp‎

Lines changed: 166 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -38,8 +38,10 @@
3838
#include <ogr_core.h>
3939
#include <ogr_geometry.h>
4040

41+
#include <algorithm>
4142
#include <cassert>
4243
#include <cerrno>
44+
#include <cstddef>
4345
#include <cstdlib>
4446
#include <cstring>
4547
#include <exception>
@@ -71,10 +73,166 @@ const unsigned int max_warnings = 500;
7173

7274
/* ================================================== */
7375

76+
/**
77+
* Is the point inside the ring? The envelope of the ring has to be given,
78+
* because calculating it is expensive and we need it several times.
79+
*
80+
* This uses the simple point-in-ring test instead of a GEOS operation,
81+
* because it has to work on invalid geometries, too.
82+
*/
83+
[[nodiscard]] bool point_in_ring(const OGRLinearRing* ring, const OGREnvelope& envelope, const OGRPoint& point) {
84+
if (point.getX() < envelope.MinX || point.getX() > envelope.MaxX ||
85+
point.getY() < envelope.MinY || point.getY() > envelope.MaxY) {
86+
return false;
87+
}
88+
89+
return ring->isPointInRing(&point, false);
90+
}
91+
92+
/**
93+
* Find the interior rings of the polygon that can not be holes of that
94+
* polygon, because they are outside the exterior ring or inside another
95+
* hole. Only the first point of each interior ring is looked at, so rings
96+
* that are only partly outside are not detected. Those cases show up as
97+
* intersections and are reported elsewhere.
98+
*/
99+
[[nodiscard]] std::vector<bool> find_misplaced_holes(const OGRPolygon& polygon) {
100+
const auto num_rings = static_cast<std::size_t>(polygon.getNumInteriorRings());
101+
102+
std::vector<OGREnvelope> envelopes(num_rings);
103+
std::vector<OGRPoint> first_points(num_rings);
104+
for (std::size_t i = 0; i < num_rings; ++i) {
105+
const OGRLinearRing* ring = polygon.getInteriorRing(static_cast<int>(i));
106+
assert(ring);
107+
ring->getEnvelope(&envelopes[i]);
108+
ring->getPoint(0, &first_points[i]);
109+
}
110+
111+
const OGRLinearRing* exterior_ring = polygon.getExteriorRing();
112+
assert(exterior_ring);
113+
OGREnvelope exterior_envelope;
114+
exterior_ring->getEnvelope(&exterior_envelope);
115+
116+
// rings outside the exterior ring ("Hole lies outside shell")
117+
std::vector<bool> outside(num_rings, false);
118+
for (std::size_t i = 0; i < num_rings; ++i) {
119+
outside[i] = !point_in_ring(exterior_ring, exterior_envelope, first_points[i]);
120+
}
121+
122+
// rings inside another hole ("Holes are nested"). Rings that are outside
123+
// the exterior ring are no holes at all, so they are not considered as
124+
// the containing ring here.
125+
std::vector<bool> misplaced{outside};
126+
for (std::size_t i = 0; i < num_rings; ++i) {
127+
if (misplaced[i]) {
128+
continue;
129+
}
130+
for (std::size_t j = 0; j < num_rings; ++j) {
131+
if (i != j && !outside[j] && envelopes[j].Contains(envelopes[i]) &&
132+
point_in_ring(polygon.getInteriorRing(static_cast<int>(j)), envelopes[j], first_points[i])) {
133+
misplaced[i] = true;
134+
break;
135+
}
136+
}
137+
}
138+
139+
return misplaced;
140+
}
141+
142+
/**
143+
* organizePolygons() decides which rings are holes based on their direction
144+
* and only checks the bounding box when assigning a hole to a polygon. So a
145+
* ring whose coastline was mapped the wrong way round can end up as a hole
146+
* of a polygon it isn't inside of or as a hole inside another hole. GEOS
147+
* then reports "Hole lies outside shell" or "Holes are nested" and the whole
148+
* polygon (possibly a whole continent) is invalid.
149+
*
150+
* Those rings are not holes but land mapped the wrong way round. Take them
151+
* out of the polygon, report them, and turn them around into land polygons
152+
* of their own.
153+
*/
154+
[[nodiscard]] std::unique_ptr<OGRPolygon> fix_misplaced_holes(std::unique_ptr<OGRPolygon> polygon,
155+
polygon_vector_type* polygons,
156+
OutputDatabase& output,
157+
unsigned int* turned_around) {
158+
const std::vector<bool> misplaced = find_misplaced_holes(*polygon);
159+
160+
if (std::find(misplaced.cbegin(), misplaced.cend(), true) == misplaced.cend()) {
161+
return polygon;
162+
}
163+
164+
auto fixed_polygon = std::make_unique<OGRPolygon>();
165+
fixed_polygon->addRingDirectly(polygon->getExteriorRing()->clone());
166+
167+
for (int i = 0; i < polygon->getNumInteriorRings(); ++i) {
168+
const OGRLinearRing* interior_ring = polygon->getInteriorRing(i);
169+
assert(interior_ring);
170+
171+
if (!misplaced[static_cast<std::size_t>(i)]) {
172+
fixed_polygon->addRingDirectly(interior_ring->clone());
173+
continue;
174+
}
175+
176+
auto ring = std::unique_ptr<OGRLinearRing>(interior_ring->clone());
177+
ring->reversePoints();
178+
179+
auto ls = std::unique_ptr<OGRLineString>(OGRGeometryFactory::forceToLineString(ring->clone())->toLineString());
180+
output.add_error_line(std::move(ls), "direction");
181+
182+
auto island = std::make_unique<OGRPolygon>();
183+
island->addRingDirectly(ring.release());
184+
island->assignSpatialReference(srs.wgs84());
185+
polygons->push_back(std::move(island));
186+
187+
++(*turned_around);
188+
}
189+
190+
fixed_polygon->assignSpatialReference(srs.wgs84());
191+
return fixed_polygon;
192+
}
193+
194+
/**
195+
* Add the polygon to the list of polygons, trying to fix it if it isn't
196+
* valid.
197+
*/
198+
void add_polygon_to(polygon_vector_type* polygons,
199+
std::unique_ptr<OGRPolygon> polygon,
200+
OutputDatabase& output,
201+
unsigned int* warnings, unsigned int* errors,
202+
unsigned int* turned_around) {
203+
if (polygon->IsValid()) {
204+
polygons->push_back(std::move(polygon));
205+
return;
206+
}
207+
208+
if (polygon->getNumInteriorRings() > 0) {
209+
polygon = fix_misplaced_holes(std::move(polygon), polygons, output, turned_around);
210+
if (polygon->IsValid()) {
211+
polygons->push_back(std::move(polygon));
212+
return;
213+
}
214+
}
215+
216+
auto* ring = polygon->getExteriorRing()->clone();
217+
auto ls = std::unique_ptr<OGRLineString>(OGRGeometryFactory::forceToLineString(ring)->toLineString());
218+
output.add_error_line(std::move(ls), "invalid");
219+
220+
std::unique_ptr<OGRGeometry> buf0{polygon->Buffer(0)};
221+
if (buf0 && buf0->getGeometryType() == wkbPolygon && buf0->IsValid()) {
222+
buf0->assignSpatialReference(srs.wgs84());
223+
polygons->push_back(static_cast_unique_ptr<OGRPolygon>(std::move(buf0)));
224+
(*warnings)++;
225+
} else {
226+
std::cerr << "Ignoring invalid polygon geometry.\n";
227+
(*errors)++;
228+
}
229+
}
230+
74231
void add_polygons_in_multi_to(polygon_vector_type *polygons,
75232
std::unique_ptr<OGRGeometry> mega_geometry,
76233
OutputDatabase& output,
77-
unsigned int* warnings, unsigned int* errors) {
234+
unsigned int* warnings, unsigned int* errors,
235+
unsigned int* turned_around) {
78236
// This isn't an owning pointer on purpose. We are going to "steal" parts
79237
// of the geometry a few lines below but only mark them as unowned farther
80238
// 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,
87245
assert(geom);
88246
assert(geom->getGeometryType() == wkbPolygon);
89247
std::unique_ptr<OGRPolygon> p{static_cast<OGRPolygon*>(geom)};
90-
if (p->IsValid()) {
91-
polygons->push_back(std::move(p));
92-
} else {
93-
auto* ring = p->getExteriorRing()->clone();
94-
auto ls = std::unique_ptr<OGRLineString>(OGRGeometryFactory::forceToLineString(ring)->toLineString());
95-
output.add_error_line(std::move(ls), "invalid");
96-
std::unique_ptr<OGRGeometry> buf0{p->Buffer(0)};
97-
if (buf0 && buf0->getGeometryType() == wkbPolygon && buf0->IsValid()) {
98-
buf0->assignSpatialReference(srs.wgs84());
99-
polygons->push_back(static_cast_unique_ptr<OGRPolygon>(std::move(buf0)));
100-
(*warnings)++;
101-
} else {
102-
std::cerr << "Ignoring invalid polygon geometry.\n";
103-
(*errors)++;
104-
}
105-
}
248+
add_polygon_to(polygons, std::move(p), output, warnings, errors, turned_around);
106249
}
107250

108251
mega_multipolygon->removeGeometry(-1, FALSE);
@@ -112,7 +255,7 @@ void add_polygons_in_multi_to(polygon_vector_type *polygons,
112255
/**
113256
* This function assembles all the coastline rings into one huge multipolygon.
114257
*/
115-
polygon_vector_type create_polygons(CoastlineRingCollection& coastline_rings, OutputDatabase& output, unsigned int* warnings, unsigned int* errors) {
258+
polygon_vector_type create_polygons(CoastlineRingCollection& coastline_rings, OutputDatabase& output, unsigned int* warnings, unsigned int* errors, unsigned int* turned_around) {
116259
std::vector<OGRGeometry*> all_polygons = coastline_rings.add_polygons_to_vector();
117260

118261
if (all_polygons.empty()) {
@@ -136,16 +279,11 @@ polygon_vector_type create_polygons(CoastlineRingCollection& coastline_rings, Ou
136279
polygon_vector_type polygons;
137280

138281
if (mega_geometry->getGeometryType() == wkbPolygon) {
139-
if (mega_geometry->IsValid()) {
140-
polygons.push_back(static_cast_unique_ptr<OGRPolygon>(std::move(mega_geometry)));
141-
} else {
142-
std::cerr << "Ignoring invalid polygon geometry.\n";
143-
(*errors)++;
144-
}
282+
add_polygon_to(&polygons, static_cast_unique_ptr<OGRPolygon>(std::move(mega_geometry)), output, warnings, errors, turned_around);
145283
} else if (mega_geometry->getGeometryType() != wkbMultiPolygon) {
146284
throw std::runtime_error{"mega geometry isn't a (multi)polygon. Something is very wrong!"};
147285
} else {
148-
add_polygons_in_multi_to(&polygons, std::move(mega_geometry), output, warnings, errors);
286+
add_polygons_in_multi_to(&polygons, std::move(mega_geometry), output, warnings, errors, turned_around);
149287
}
150288

151289
return polygons;
@@ -354,15 +492,16 @@ int main(int argc, char *argv[]) {
354492
if (options.output_polygons != output_polygon_type::none || options.output_lines) {
355493
try {
356494
vout << "Create polygons...\n";
357-
CoastlinePolygons coastline_polygons{create_polygons(coastline_rings, *output_database, &warnings, &errors), \
495+
unsigned int turned_around = 0;
496+
CoastlinePolygons coastline_polygons{create_polygons(coastline_rings, *output_database, &warnings, &errors, &turned_around), \
358497
*output_database, \
359498
options.bbox_overlap, \
360499
options.max_points_in_polygon};
361500

362501
stats.land_polygons_before_split = coastline_polygons.num_polygons();
363502

364503
vout << "Fixing coastlines going the wrong way...\n";
365-
stats.rings_turned_around = coastline_polygons.fix_direction();
504+
stats.rings_turned_around = turned_around + coastline_polygons.fix_direction();
366505
vout << " Turned " << stats.rings_turned_around << " polygons around.\n";
367506
warnings += stats.rings_turned_around;
368507

‎test/t/invalid-direction-island-in-bbox-of-many.sh‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
#!/bin/sh
2+
#-----------------------------------------------------------------------------
3+
#
4+
# Same as invalid-direction-island-in-bbox, but with a second, correctly
5+
# mapped island. This means organizePolygons() returns a multipolygon which
6+
# is handled by a different code path.
7+
# See https://github.com/osmcode/osmcoastline/issues/41
8+
#
9+
#-----------------------------------------------------------------------------
10+
11+
# shellcheck source=test/init.sh
12+
. "$1/test/init.sh"
13+
14+
set -x
15+
16+
#-----------------------------------------------------------------------------
17+
18+
"$BIN_DIR/src/nodegrid2opl" << 'NODES' >"$INPUT"
19+
20+
0---------------1
21+
| |
22+
| | a---b
23+
| 3--------2 | |
24+
| | d---c
25+
| | 6---7
26+
| | | |
27+
| | 9---8
28+
| |
29+
5------4
30+
31+
NODES
32+
33+
cat <<'OSM' >>"$INPUT"
34+
w200 v1 Tnatural=coastline Nn100,n105,n104,n103,n102,n101,n100
35+
w201 v1 Tnatural=coastline Nn106,n107,n108,n109,n106
36+
w202 v1 Tnatural=coastline Nn110,n113,n112,n111,n110
37+
OSM
38+
39+
#-----------------------------------------------------------------------------
40+
41+
"$OSMC" --verbose --overwrite --srs="$SRID" --output-database="$DB" "$INPUT" >"$LOG" 2>&1
42+
RC=$?
43+
set -e
44+
45+
test $RC -eq 1
46+
47+
grep 'Turned 1 polygons around.$' "$LOG"
48+
49+
grep '^There were 1 warnings.$' "$LOG"
50+
grep '^There were 0 errors.$' "$LOG"
51+
52+
check_count land_polygons 3;
53+
check_count error_points 0;
54+
check_count error_lines 1;
55+
56+
echo "SELECT InsertEpsgSrid(4326);" | $SQL
57+
58+
# the island that was turned around
59+
echo "SELECT AsText(Transform(geometry, 4326)) FROM land_polygons;" | $SQL \
60+
| grep -F 'POLYGON((1.15 1.94, 1.19 1.94, 1.19 1.92, 1.15 1.92, 1.15 1.94))'
61+
62+
echo "SELECT AsText(Transform(geometry, 4326)), osm_id, error FROM error_lines;" | $SQL \
63+
| grep -F '|0|direction'
64+
65+
#-----------------------------------------------------------------------------

‎test/t/invalid-direction-island-in-bbox.sh‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
#!/bin/sh
2+
#-----------------------------------------------------------------------------
3+
#
4+
# Invalid island with coastline going the wrong direction inside the
5+
# bounding box of (but outside) a larger land polygon. Used to result in a
6+
# "Hole lies outside shell" error from GEOS which made the larger polygon
7+
# invalid. See https://github.com/osmcode/osmcoastline/issues/41
8+
#
9+
#-----------------------------------------------------------------------------
10+
11+
# shellcheck source=test/init.sh
12+
. "$1/test/init.sh"
13+
14+
set -x
15+
16+
#-----------------------------------------------------------------------------
17+
18+
"$BIN_DIR/src/nodegrid2opl" << 'NODES' >"$INPUT"
19+
20+
0---------------1
21+
| |
22+
| |
23+
| 3--------2
24+
| |
25+
| | 6---7
26+
| | | |
27+
| | 9---8
28+
| |
29+
5------4
30+
31+
NODES
32+
33+
cat <<'OSM' >>"$INPUT"
34+
w200 v1 Tnatural=coastline Nn100,n105,n104,n103,n102,n101,n100
35+
w201 v1 Tnatural=coastline Nn106,n107,n108,n109,n106
36+
OSM
37+
38+
#-----------------------------------------------------------------------------
39+
40+
"$OSMC" --verbose --overwrite --srs="$SRID" --output-database="$DB" "$INPUT" >"$LOG" 2>&1
41+
RC=$?
42+
set -e
43+
44+
test $RC -eq 1
45+
46+
grep 'Turned 1 polygons around.$' "$LOG"
47+
48+
grep '^There were 1 warnings.$' "$LOG"
49+
grep '^There were 0 errors.$' "$LOG"
50+
51+
check_count land_polygons 2;
52+
check_count error_points 0;
53+
check_count error_lines 1;
54+
55+
echo "SELECT InsertEpsgSrid(4326);" | $SQL
56+
57+
echo "SELECT AsText(Transform(geometry, 4326)) FROM land_polygons;" | $SQL \
58+
| grep -F 'POLYGON((1.15 1.94, 1.19 1.94, 1.19 1.92, 1.15 1.92, 1.15 1.94))'
59+
60+
echo "SELECT AsText(Transform(geometry, 4326)), osm_id, error FROM error_lines;" | $SQL \
61+
| grep -F '|0|direction'
62+
63+
#-----------------------------------------------------------------------------

0 commit comments

Comments
 (0)