[geos-commits] [SCM] GEOS branch main updated. fc3eca9063096a665700651d6523f56de4dae415
git at osgeo.org
git at osgeo.org
Thu Nov 5 06:10:37 PST 2026
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, main has been updated
via fc3eca9063096a665700651d6523f56de4dae415 (commit)
from 8bace6161f6593643e45875d7972247699357720 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit fc3eca9063096a665700651d6523f56de4dae415
Author: Daniel Baston <dbaston at gmail.com>
Date: Thu Nov 5 09:10:08 2026 -0500
GridIntersection: robustness fixes (#1589)
Co-authored-by: Alejandro Pacheco <ale at alepacheco.dev>
diff --git a/include/geos/operation/grid/Grid.h b/include/geos/operation/grid/Grid.h
index d94d80916..e56d7dd84 100644
--- a/include/geos/operation/grid/Grid.h
+++ b/include/geos/operation/grid/Grid.h
@@ -19,6 +19,7 @@
#include <vector>
#include <geos/geom/Envelope.h>
+#include <geos/util/IllegalArgumentException.h>
namespace geos::operation::grid {
struct infinite_extent
@@ -69,6 +70,12 @@ class GEOS_DLL Grid
, m_num_cols{ 2 * extent_tag::padding + (extent.getMaxX() > extent.getMinX() ? static_cast<size_t>(std::round(extent.getWidth() / dx)) : 0) }
{
static_assert(std::is_same_v<extent_tag, infinite_extent>);
+ if (!isEmpty()) {
+ if (m_domain.getMinX() >= m_extent.getMinX() || m_domain.getMaxX() <= m_extent.getMaxX() ||
+ m_domain.getMinY() >= m_extent.getMinY() || m_domain.getMaxY() <= m_extent.getMaxY()) {
+ throw geos::util::IllegalArgumentException("Domain of infinite grid must be larger than extent");
+ }
+ }
}
diff --git a/src/operation/grid/Grid.cpp b/src/operation/grid/Grid.cpp
index d4db6659b..bd63dff71 100644
--- a/src/operation/grid/Grid.cpp
+++ b/src/operation/grid/Grid.cpp
@@ -41,11 +41,9 @@ Grid<infinite_extent>::getCellEnvelope(size_t row, size_t col) const
{
double cell_xmin, cell_xmax, cell_ymin, cell_ymax;
- constexpr double PADDED_CELL_SIZE = 1.0;
-
if (col == 0) {
// Left-side padding
- cell_xmin = std::min(xmin() - PADDED_CELL_SIZE, m_domain.getMinX());
+ cell_xmin = m_domain.getMinX();
} else if (col == getNumCols() - 1) {
// Right-side padding
cell_xmin = xmax(); // because rightmost col of regular may have different width from others
@@ -57,7 +55,7 @@ Grid<infinite_extent>::getCellEnvelope(size_t row, size_t col) const
switch (getNumCols() - col) {
case 1:
// Right-side padding
- cell_xmax = std::max(xmax() + PADDED_CELL_SIZE, m_domain.getMaxX());
+ cell_xmax = m_domain.getMaxX();
break;
case 2:
// Rightmost internal cell
@@ -69,7 +67,7 @@ Grid<infinite_extent>::getCellEnvelope(size_t row, size_t col) const
if (row == 0) {
// Top row padding
- cell_ymax = std::max(ymax() + PADDED_CELL_SIZE, m_domain.getMaxY());
+ cell_ymax = m_domain.getMaxY();
} else if (row == getNumRows() - 1) {
// Bottom row padding
cell_ymax = ymin(); // because bottom row of regular may have different height from others
@@ -80,7 +78,7 @@ Grid<infinite_extent>::getCellEnvelope(size_t row, size_t col) const
switch (getNumRows() - row) {
case 1:
// Bottom row padding
- cell_ymin = std::min(xmin() - PADDED_CELL_SIZE, m_domain.getMinY());
+ cell_ymin = m_domain.getMinY();
break;
case 2:
// Bottom-most internal row
diff --git a/src/operation/grid/GridIntersection.cpp b/src/operation/grid/GridIntersection.cpp
index c60c3aaec..ef94ce5ad 100644
--- a/src/operation/grid/GridIntersection.cpp
+++ b/src/operation/grid/GridIntersection.cpp
@@ -34,6 +34,24 @@ using geos::geom::CoordinateXY;
namespace geos::operation::grid {
+static Envelope
+getDomain(const Envelope& gridExtent, const Envelope& geomExtent)
+{
+ constexpr double PADDED_CELL_SIZE = 1.0;
+ constexpr double PADDED_CELL_RELATIVE_SIZE = 1e-12;
+
+ Envelope padded(gridExtent);
+ padded.expandToInclude(geomExtent);
+
+ const double xmag = std::max(std::abs(padded.getMinX()), std::abs(padded.getMaxX()));
+ const double ymag = std::max(std::abs(padded.getMinY()), std::abs(padded.getMaxY()));
+
+ padded.expandBy(std::max(PADDED_CELL_SIZE, PADDED_CELL_RELATIVE_SIZE * xmag),
+ std::max(PADDED_CELL_SIZE, PADDED_CELL_RELATIVE_SIZE * ymag));
+
+ return padded;
+}
+
std::shared_ptr<Matrix<float>>
GridIntersection::getIntersectionFractions(const Grid<bounded_extent>& raster_grid, const Geometry& g)
{
@@ -111,7 +129,7 @@ GridIntersection::processingRegion(const Envelope& raster_extent, const Geometry
}
GridIntersection::GridIntersection(const Grid<bounded_extent>& raster_grid, const Geometry& g, const std::shared_ptr<Matrix<float>>& cov)
- : m_geometry_grid{ make_infinite(raster_grid, *g.getEnvelopeInternal()) }
+ : m_geometry_grid{ make_infinite(raster_grid, getDomain(raster_grid.getExtent(), *g.getEnvelopeInternal())) }
, m_results{ cov ? cov : std::make_shared<Matrix<float>>(raster_grid.getNumRows(), raster_grid.getNumCols()) }
, m_first_geom{ true }
, m_areal{ false }
@@ -130,7 +148,7 @@ GridIntersection::GridIntersection(const Grid<bounded_extent>& raster_grid, cons
}
GridIntersection::GridIntersection(const Grid<bounded_extent>& raster_grid, const Envelope& box, const std::shared_ptr<Matrix<float>>& cov)
- : m_geometry_grid{ make_infinite(raster_grid, box) }
+ : m_geometry_grid{ make_infinite(raster_grid, getDomain(raster_grid.getExtent(), box)) }
, m_results{ cov ? cov : std::make_shared<Matrix<float>>(raster_grid.getNumRows(), raster_grid.getNumCols()) }
, m_first_geom { true }
, m_areal{ false }
@@ -182,6 +200,10 @@ GridIntersection::processRectangularRing(const Envelope& box, bool exterior_ring
auto ring_grid = getSubgrid(box, m_geometry_grid);
+ if (ring_grid.isEmpty()) {
+ return;
+ }
+
auto row_min = ring_grid.getRow(box.getMaxY());
auto row_max = ring_grid.getRow(box.getMinY());
auto col_min = ring_grid.getColumn(box.getMinX());
@@ -590,8 +612,7 @@ GridIntersection::subdividePolygon(const Grid<bounded_extent>& p_grid, const Geo
const auto cropped_grid = p_grid.shrinkToFit(p_grid.getExtent().intersection(*g.getEnvelopeInternal()), false);
- geom::Envelope gridExtent = *g.getEnvelopeInternal();
- gridExtent.expandBy(1);
+ geom::Envelope gridExtent = getDomain(cropped_grid.getExtent(), *g.getEnvelopeInternal());
const Grid<infinite_extent> cell_grid = make_infinite(cropped_grid, gridExtent);
Matrix<std::unique_ptr<Cell>> cells(cell_grid.getNumRows(), cell_grid.getNumCols());
diff --git a/tests/unit/operation/grid/GridIntersectionTest.cpp b/tests/unit/operation/grid/GridIntersectionTest.cpp
index 856bbe60d..940d54f3a 100644
--- a/tests/unit/operation/grid/GridIntersectionTest.cpp
+++ b/tests/unit/operation/grid/GridIntersectionTest.cpp
@@ -1013,4 +1013,83 @@ void object::test<54>()
ensure_equals("check area", tot/100, g->getArea(), 1e-6);
}
+template<>
+template<>
+void object::test<55>()
+{
+ set_test_name("getIntersectionFractions: traversal stays in grid when a vertex is the top corner of the geometry envelope");
+
+ // The vertex at (-259.05369950585896 811.4225290921494) is the upper-left corner of the polygon envelope
+ // and lies on a cell boundary, to within rounding.
+ Envelope e(-300.91149168825694, -175.33811514106293, 225.41343853857734, 602.13356818015927);
+
+ Grid<bounded_extent> ext(e, e.getWidth() / 6, e.getHeight() / 18);
+
+ auto g = wkt_reader_.read("POLYGON ((-70.69363468506793 685.8491525449554, -133.48032295866494 748.6358408185523, -175.33811514106293 790.4936330009504, -259.05369950585896 811.4225290921494, -238.12480341465994 706.7780486361544, -259.05369950585896 664.9202564537563, -259.05369950585896 581.2046720889604, -196.26701123226195 643.9913603625573, -112.55142686746595 623.0624642713583, -70.69363468506793 685.8491525449554))");
+
+ auto fracs = GridIntersection::getIntersectionFractions(ext, *g);
+
+ auto inside = g->intersection(g->getFactory()->toGeometry(&e).get());
+ check_area(*fracs, ext, *inside);
+}
+
+template<>
+template<>
+void object::test<56>()
+{
+ set_test_name("getIntersectionFractions: traversal stays in grid when a vertex is the bottom corner of the geometry envelope");
+
+ // The vertex at (68.44482903503193 -90.56515438144179) is the lower-right corner of the polygon envelope
+ // and lies on a cell boundary, to within rounding.
+ Envelope e(-10.004484448704343, 152.12409675101728, 35.385086543444757, 73.717768564062396);
+ Grid<bounded_extent> ext(e, e.getWidth() / 31, e.getHeight() / 14);
+
+ auto g = wkt_reader_.read("POLYGON ((173.04391368001362 60.02752498527039, 78.9047374995301 109.31240186892164, 89.36464596402827 257.1670325198754, -25.694347145451598 213.35825306774097, -93.68375216468971 161.3353274683313, -104.14366062918788 90.14606085861283, -119.83352332593513 24.432891680411146, -51.84411830669703 5.266550670102323, -20.464392913202513 -52.23247236082415, 68.44482903503193 -90.56515438144179, 173.04391368001362 -54.97052107658256, 173.04391368001362 60.02752498527039))");
+
+ auto fracs = GridIntersection::getIntersectionFractions(ext, *g);
+
+ auto inside = g->intersection(g->getFactory()->toGeometry(&e).get());
+ check_area(*fracs, ext, *inside);
+}
+
+template<>
+template<>
+void object::test<57>()
+{
+ set_test_name("subdividePolygon: traversal stays in grid when padding is lost to rounding at large coordinates");
+
+ // At 1e22, adding a fixed padding of 1 to the geometry envelope has no effect,
+ // so the vertex at (1e22 0.9) would lie on the outer edge of a padding cell.
+ Envelope e(-0.375, 2.1875, -0.375, 2.1875);
+ Grid<bounded_extent> ext(e, e.getWidth() / 64, e.getHeight() / 2);
+
+ auto g = wkt_reader_.read("POLYGON ((0.1 0.1, 0.1 0.2, 0.95 0.95, 1e22 0.9, 0.1 0.1))");
+
+ auto subd = GridIntersection::subdividePolygon(ext, *g, false);
+
+ auto inside = g->intersection(g->getFactory()->toGeometry(&e).get());
+ check_subdivided_polygon(*inside, *subd);
+}
+
+template<>
+template<>
+void object::test<58>()
+{
+ set_test_name("getIntersectionFractions: box sharing an edge with the grid");
+
+ Envelope e(0, 10, 0, 10);
+ Grid<bounded_extent> ext(e, 1, 1);
+
+ for (const auto& box : {
+ Envelope(10, 20, 0, 10),
+ Envelope(0, 10, 10, 20),
+ Envelope(-10, 0, -10, 0) })
+ {
+ auto fracs = GridIntersection::getIntersectionFractions(ext, box);
+ auto tot = std::accumulate(fracs->begin(), fracs->end(), 0.0);
+
+ ensure_equals(tot, 0.0);
+ }
+}
+
}
diff --git a/tests/unit/operation/grid/GridTest.cpp b/tests/unit/operation/grid/GridTest.cpp
index 0df5ea8a7..c142c4911 100644
--- a/tests/unit/operation/grid/GridTest.cpp
+++ b/tests/unit/operation/grid/GridTest.cpp
@@ -314,8 +314,10 @@ void object::test<17>()
set_test_name("Crop infinite grid while calculating cell boundaries relative to parent");
double res = 1.0/3601;
+ Envelope domain = global;
+ domain.expandBy(1);
- Grid<infinite_extent> grid{global, res, res, global};
+ Grid<infinite_extent> grid{global, res, res, domain};
const auto env0 = grid.getCellEnvelope(30, 30);
-----------------------------------------------------------------------
Summary of changes:
include/geos/operation/grid/Grid.h | 7 ++
src/operation/grid/Grid.cpp | 10 ++-
src/operation/grid/GridIntersection.cpp | 29 ++++++--
tests/unit/operation/grid/GridIntersectionTest.cpp | 79 ++++++++++++++++++++++
tests/unit/operation/grid/GridTest.cpp | 4 +-
5 files changed, 118 insertions(+), 11 deletions(-)
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list