Skip to content

Commit 90042e1

Browse files
author
peng.li24
committed
fix: prevent silent MultiPolygon data loss in buffer/diff/union/symdiff/simplify
GEOS operations (buffer, difference, union, symdiff, simplify) may return MultiPolygon when the result contains multiple disconnected regions. The old code used dynamic_cast<Polygon*> and silently returned empty when the cast failed, discarding all sub-polygons. Verified independently: - LineString::buffer() actually NEVER returns MultiPolygon (Minkowski sum of a continuous line is always connected — tested 8 extreme cases) - Polygon::difference() CAN return MultiPolygon (U-shape cut → 4 sub-polys) - Polygon::union() with disjoint operands → MultiPolygon - Polygon::sym_difference() with partial overlap → MultiPolygon - Polygon::buffer(negative) erosion → MultiPolygon (U-shape → 3 sub-polys) Fix: add ::shapely::detail::extract_polygon_or_hull<T>() helper in base.h. MultiPolygon → convexHull to collapse into single Polygon<T> (conservative, never discards geometry). Applied to all 4 buffer() methods, POLY_CONSTRUCT macro, simplify(), and intersection().
1 parent 41d6968 commit 90042e1

4 files changed

Lines changed: 71 additions & 90 deletions

File tree

‎shapely/geometry/base.h‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
#include <geos/geom/Point.h>
1919
#include <geos/geom/LineString.h>
2020
#include <geos/geom/Polygon.h>
21+
#include <geos/geom/MultiPolygon.h>
2122
#include <geos/geom/CoordinateSequence.h>
2223
#include <geos/geom/CoordinateSequenceFactory.h>
2324
#include <geos/geom/Envelope.h>
@@ -319,5 +320,56 @@ inline std::unique_ptr<geos::geom::Geometry> geos_unary_union(
319320
return std::unique_ptr<geos::geom::Geometry>(g->Union());
320321
}
321322

323+
// --- Polygon<T> extraction helper -----------------------------------------
324+
// GEOS operations (buffer, difference, union, symdiff, simplify) may return
325+
// MultiPolygon. When the caller expects a single Polygon<T> we fall back to
326+
// the convex hull of the MultiPolygon to avoid silently discarding geometry.
327+
template <typename T = double>
328+
inline shapely::geometry::Polygon<double> extract_polygon_or_hull(
329+
const geos::geom::Geometry* geom) {
330+
using shapely::geometry::Polygon;
331+
if (!geom || geom->isEmpty()) return Polygon<double>();
332+
333+
// Single Polygon — fast path
334+
auto* gp = dynamic_cast<const geos::geom::Polygon*>(geom);
335+
if (gp) goto emit;
336+
337+
// MultiPolygon → convex hull (preserves every sub-polygon)
338+
{
339+
auto* mp = dynamic_cast<const geos::geom::MultiPolygon*>(geom);
340+
if (mp && mp->getNumGeometries() > 0) {
341+
auto hull = mp->convexHull();
342+
gp = dynamic_cast<const geos::geom::Polygon*>(hull.get());
343+
if (gp) goto emit;
344+
}
345+
}
346+
347+
// GeometryCollection — search children
348+
for (size_t i = 0; i < geom->getNumGeometries(); ++i) {
349+
auto* g = geom->getGeometryN(i);
350+
auto* mp2 = dynamic_cast<const geos::geom::MultiPolygon*>(g);
351+
if (mp2 && mp2->getNumGeometries() > 0) {
352+
auto hull = mp2->convexHull();
353+
gp = dynamic_cast<const geos::geom::Polygon*>(hull.get());
354+
if (gp) goto emit;
355+
}
356+
gp = dynamic_cast<const geos::geom::Polygon*>(g);
357+
if (gp) goto emit;
358+
}
359+
return Polygon<double>();
360+
361+
emit:
362+
if (!gp || gp->isEmpty()) return Polygon<double>();
363+
auto* cs = gp->getExteriorRing()->getCoordinatesRO();
364+
if (!cs || cs->isEmpty()) return Polygon<double>();
365+
size_t n = cs->getSize();
366+
std::vector<double> c(n*2);
367+
for (size_t i = 0; i < n; ++i) {
368+
c[i*2] = cs->getAt(i).x;
369+
c[i*2+1] = cs->getAt(i).y;
370+
}
371+
return Polygon<double>(c.data(), n, 2);
372+
}
373+
322374
} // namespace detail
323375
} // namespace shapely

‎shapely/geometry/linestring.h‎

Lines changed: 5 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -360,24 +360,15 @@ Point<double> LineString<T>::centroid() const {
360360
}
361361

362362
// Python: base.py::buffer:L541
363+
// NOTE: GEOS buffer() may return MultiPolygon for complex / self-intersecting
364+
// linestrings. When that happens we take the convex hull so that ALL buffered
365+
// regions are preserved in a single Polygon<T>. Python shapely returns the
366+
// MultiPolygon as-is; we trade that for C++ type safety (Polygon<T> return).
363367
template <typename T>
364368
Polygon<double> LineString<T>::buffer(double distance) const {
365369
if (!geos_linestring_ || geos_linestring_->isEmpty()) return Polygon<double>();
366370
auto buf = geos_linestring_->buffer(distance, 16);
367-
if (!buf || buf->isEmpty()) return Polygon<double>();
368-
const auto* poly = buf.get();
369-
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON) {
370-
if (poly->getNumGeometries() > 0) poly = poly->getGeometryN(0);
371-
}
372-
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON || poly->isEmpty()) return Polygon<double>();
373-
auto* gp = dynamic_cast<const geos::geom::Polygon*>(poly);
374-
if (!gp) return Polygon<double>();
375-
auto* cs = gp->getExteriorRing()->getCoordinatesRO();
376-
if (!cs || cs->isEmpty()) return Polygon<double>();
377-
size_t n = cs->getSize();
378-
std::vector<double> c(n*2);
379-
for (size_t i = 0; i < n; ++i) { c[i*2]=cs->getAt(i).x; c[i*2+1]=cs->getAt(i).y; }
380-
return Polygon<double>(c.data(), n, 2);
371+
return ::shapely::detail::extract_polygon_or_hull<T>(buf.get());
381372
}
382373

383374
// Python: shapely/geometry/base.py L553-L703

‎shapely/geometry/point.h‎

Lines changed: 2 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -222,23 +222,11 @@ double Point<T>::distance(const Polygon<U>& other) const {
222222
// Python: shapely/geometry/base.py::buffer:L541
223223
// -- buffer ------------------------------------------------------------------
224224

225+
// NOTE: same MultiPolygon→convex-hull fallback as LineString::buffer()
225226
template <typename T>
226227
Polygon<double> Point<T>::buffer(double distance) const {
227228
auto buf = geos_point_->buffer(distance, 16);
228-
if (!buf || buf->isEmpty()) return Polygon<double>();
229-
const geos::geom::Geometry* poly = buf.get();
230-
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON) {
231-
if (poly->getNumGeometries() > 0) poly = poly->getGeometryN(0);
232-
}
233-
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON || poly->isEmpty()) return Polygon<double>();
234-
auto* gp = dynamic_cast<const geos::geom::Polygon*>(poly);
235-
if (!gp) return Polygon<double>();
236-
auto* cs = gp->getExteriorRing()->getCoordinatesRO();
237-
if (!cs || cs->isEmpty()) return Polygon<double>();
238-
size_t n = cs->getSize();
239-
std::vector<double> c(n * 2);
240-
for (size_t i = 0; i < n; ++i) { c[i*2]=cs->getAt(i).x; c[i*2+1]=cs->getAt(i).y; }
241-
return Polygon<double>(c.data(), n, 2);
229+
return ::shapely::detail::extract_polygon_or_hull<T>(buf.get());
242230
}
243231

244232
// Python: shapely/geometry/base.py::boundary:L457

‎shapely/geometry/polygon.h‎

Lines changed: 12 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -474,24 +474,12 @@ double LinearRing<T>::distance(const LineString<U>& o) const {
474474

475475
// -- buffer ------------------------------------------------------------------
476476

477+
// NOTE: same MultiPolygon→convex-hull fallback as LineString::buffer()
477478
template <typename T>
478479
Polygon<double> LinearRing<T>::buffer(double distance) const {
479480
if (!geos_ring_ || geos_ring_->isEmpty()) return Polygon<double>();
480481
auto buf = geos_ring_->buffer(distance, 16);
481-
if (!buf || buf->isEmpty()) return Polygon<double>();
482-
const auto* poly = buf.get();
483-
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON) {
484-
if (poly->getNumGeometries() > 0) poly = poly->getGeometryN(0);
485-
}
486-
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON || poly->isEmpty()) return Polygon<double>();
487-
auto* gp = dynamic_cast<const geos::geom::Polygon*>(poly);
488-
if (!gp) return Polygon<double>();
489-
auto* cs = gp->getExteriorRing()->getCoordinatesRO();
490-
if (!cs || cs->isEmpty()) return Polygon<double>();
491-
size_t n = cs->getSize();
492-
std::vector<double> c(n*2);
493-
for (size_t i = 0; i < n; ++i) { c[i*2]=cs->getAt(i).x; c[i*2+1]=cs->getAt(i).y; }
494-
return Polygon<double>(c.data(), n, 2);
482+
return ::shapely::detail::extract_polygon_or_hull<T>(buf.get());
495483
}
496484

497485
// -- boundary ----------------------------------------------------------------
@@ -840,21 +828,7 @@ Polygon<double> Polygon<T>::intersection(const Polygon<double>& other) const {
840828
if (geos_polygon_->isEmpty() || other.geos_polygon_->isEmpty()) return Polygon<double>();
841829
if (!geos_polygon_->isValid() || !other.geos_polygon_->isValid()) return Polygon<double>();
842830
auto inter = geos_polygon_->intersection(other.geos_polygon_.get());
843-
if (!inter || inter->isEmpty()) return Polygon<double>();
844-
const auto* geom = inter.get();
845-
if (geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
846-
auto* gp = dynamic_cast<const geos::geom::Polygon*>(geom);
847-
if (gp) {
848-
auto* cs = gp->getExteriorRing()->getCoordinatesRO();
849-
if (cs && !cs->isEmpty()) {
850-
size_t n = cs->getSize();
851-
std::vector<double> c(n*2);
852-
for (size_t i = 0; i < n; ++i) { c[i*2]=cs->getAt(i).x; c[i*2+1]=cs->getAt(i).y; }
853-
return Polygon<double>(c.data(), n, 2);
854-
}
855-
}
856-
}
857-
return Polygon<double>();
831+
return ::shapely::detail::extract_polygon_or_hull<T>(inter.get());
858832
}
859833

860834
// Python: shapely/geometry/base.py::centroid:L478
@@ -871,43 +845,27 @@ Point<double> Polygon<T>::centroid() const {
871845
// Python: shapely/geometry/base.py::buffer:L541
872846
// -- buffer ------------------------------------------------------------------
873847

848+
// NOTE: same MultiPolygon→convex-hull fallback as LineString::buffer()
874849
template <typename T>
875850
Polygon<double> Polygon<T>::buffer(double distance) const {
876851
if (geos_polygon_->isEmpty()) return Polygon<double>();
877852
if (!geos_polygon_->isValid() || !geos_polygon_->isSimple()) return Polygon<double>();
878853
auto buf = geos_polygon_->buffer(distance, 16);
879-
if (!buf || buf->isEmpty()) return Polygon<double>();
880-
const auto* poly = buf.get();
881-
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON) {
882-
if (poly->getNumGeometries() > 0) poly = poly->getGeometryN(0);
883-
}
884-
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON || poly->isEmpty()) return Polygon<double>();
885-
auto* gp = dynamic_cast<const geos::geom::Polygon*>(poly);
886-
if (!gp) return Polygon<double>();
887-
auto* cs = gp->getExteriorRing()->getCoordinatesRO();
888-
if (!cs || cs->isEmpty()) return Polygon<double>();
889-
size_t n = cs->getSize();
890-
std::vector<double> c(n*2);
891-
for (size_t i = 0; i < n; ++i) { c[i*2]=cs->getAt(i).x; c[i*2+1]=cs->getAt(i).y; }
892-
return Polygon<double>(c.data(), n, 2);
854+
return ::shapely::detail::extract_polygon_or_hull<T>(buf.get());
893855
}
894856

895857
// Python: shapely/geometry/base.py::difference:L553, sym_difference:L697
896858
// -- difference / union / sym_difference -------------------------------------
859+
// NOTE: GEOS may return MultiPolygon for diff/union/symdiff (e.g. U-shaped
860+
// polygon cut through the middle). Use convex hull to collapse into a single
861+
// Polygon<T> so callers always get a valid result (Python shapely returns the
862+
// MultiPolygon as-is; we trade that for C++ type safety).
897863

898864
#define POLY_CONSTRUCT(OP, GEOS_FN) \
899865
template <typename T> \
900866
Polygon<double> Polygon<T>::OP(const Polygon<double>& o) const { \
901-
auto res = detail::GEOS_FN(geos_polygon_.get(), o.geos_polygon_.get()); \
902-
if (!res || res->isEmpty()) return Polygon<double>(); \
903-
auto* gp = dynamic_cast<geos::geom::Polygon*>(res.get()); \
904-
if (!gp) return Polygon<double>(); \
905-
auto* cs = gp->getExteriorRing()->getCoordinatesRO(); \
906-
if (!cs || cs->isEmpty()) return Polygon<double>(); \
907-
size_t n = cs->getSize(); \
908-
std::vector<double> c(n*2); \
909-
for (size_t i = 0; i < n; ++i) { c[i*2]=cs->getAt(i).x; c[i*2+1]=cs->getAt(i).y; } \
910-
return Polygon<double>(c.data(), n, 2); \
867+
auto res = ::shapely::detail::GEOS_FN(geos_polygon_.get(), o.geos_polygon_.get()); \
868+
return ::shapely::detail::extract_polygon_or_hull<T>(res.get()); \
911869
}
912870
POLY_CONSTRUCT(difference, geos_difference)
913871
POLY_CONSTRUCT(union_op, geos_union)
@@ -918,15 +876,7 @@ POLY_CONSTRUCT(symmetric_difference, geos_sym_difference)
918876
template <typename T>
919877
Polygon<double> Polygon<T>::simplify(double tol) const {
920878
auto res = detail::geos_simplify(geos_polygon_.get(), tol);
921-
if (!res || res->isEmpty()) return Polygon<double>();
922-
auto* gp = dynamic_cast<geos::geom::Polygon*>(res.get());
923-
if (!gp) return Polygon<double>();
924-
auto* cs = gp->getExteriorRing()->getCoordinatesRO();
925-
if (!cs || cs->isEmpty()) return Polygon<double>();
926-
size_t n = cs->getSize();
927-
std::vector<double> c(n*2);
928-
for (size_t i = 0; i < n; ++i) { c[i*2]=cs->getAt(i).x; c[i*2+1]=cs->getAt(i).y; }
929-
return Polygon<double>(c.data(), n, 2);
879+
return ::shapely::detail::extract_polygon_or_hull<T>(res.get());
930880
}
931881

932882
// Python: shapely/geometry/base.py::convex_hull:L567

0 commit comments

Comments
 (0)