Add GEOSGridIntersectionFractions, GEOSSubdivideByGrid

benchmarks/operation/CMakeLists.txt

@@ -27,4 +27,12 @@ if (benchmark_FOUND)
             $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/benchmarks>)
     target_link_libraries(perf_coverage_union PRIVATE
             benchmark::benchmark geos geos_cxx_flags)
+
+    add_executable(perf_grid_intersection GridIntersectionPerfTest.cpp)
+    target_include_directories(perf_grid_intersection PUBLIC
+            $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include>
+            $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+            $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/benchmarks>)
+    target_link_libraries(perf_grid_intersection
+            benchmark::benchmark geos geos_cxx_flags)
 endif()

benchmarks/operation/GridIntersectionPerfTest.cpp

@@ -0,0 +1,121 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2025 ISciences LLC
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <benchmark/benchmark.h>
+#include <BenchmarkUtils.h>
+
+#include <geos/geom/Envelope.h>
+#include <geos/geom/prep/PreparedGeometryFactory.h>
+#include <geos/operation/grid/Grid.h>
+
+#include <geos/operation/grid/GridIntersection.h>
+#include <geos/operation/intersection/Rectangle.h>
+#include <geos/operation/intersection/RectangleIntersection.h>
+
+using geos::geom::CoordinateXY;
+using geos::geom::Envelope;
+using geos::geom::Geometry;
+using Grid = geos::operation::grid::Grid<geos::operation::grid::bounded_extent>;
+
+template<bool AreaOnly>
+struct GridIntersection {
+    static double Intersection(const Envelope& env, int nx, int ny, const Geometry& g) {
+         Grid grid(env, env.getWidth() / nx, env.getHeight() / ny);
+         if constexpr (AreaOnly) {
+             auto result = geos::operation::grid::GridIntersection::grid_intersection(grid, g);
+             float area = 0;
+             for (std::size_t i = 0; i < result->rows(); i++) {
+                 for (std::size_t j = 0; j < result->cols(); j++) {
+                    area += (*result)(i, j);
+                 }
+             }
+            return static_cast<double>(area);
+         } else {
+             auto subdivided = geos::operation::grid::GridIntersection::subdivide_polygon(grid, g, true);
+             return subdivided->getArea();
+         }
+    }
+};
+
+using GridIntersectionAreaOnly = GridIntersection<true>;
+using GridIntersectionFull = GridIntersection<false>;
+
+template<bool UseRectangleIntersection>
+struct SingleIntersection {
+
+    static double Intersection(const Envelope& env, int nx, int ny, const Geometry& g) {
+        double dx = env.getWidth() / nx;
+        double dy = env.getHeight() / ny;
+
+        double x0 = env.getMinX();
+        double y0 = env.getMinY();
+
+        const auto& gfact = *g.getFactory();
+        auto prepGeom = geos::geom::prep::PreparedGeometryFactory::prepare(&g);
+
+        double area = 0;
+
+        for (int i = 0; i < nx; i++) {
+            for (int j = 0; j < ny; j++) {
+                Envelope subEnv(x0 + i*dx, x0 + (i+1)*dx, y0 + j*dy, y0 + (j+1)*dy);
+                auto cellGeom = gfact.toGeometry(&subEnv);
+                if (!prepGeom->intersects(cellGeom.get())) {
+                    continue;
+                }
+
+                std::unique_ptr<Geometry> isect;
+                if constexpr (UseRectangleIntersection) {
+                    geos::operation::intersection::Rectangle rect(x0 + i*dx, y0 + j*dy, x0 + (i+1)*dx, y0 + (j+1)*dy);
+                    isect = geos::operation::intersection::RectangleIntersection::clip(g, rect);
+                } else {
+                    isect = g.intersection(cellGeom.get());
+                }
+
+                area += isect->getArea();
+            }
+        }
+
+        return area;
+    }
+};
+
+using PolygonIntersection = SingleIntersection<false>;
+using RectangleIntersection = SingleIntersection<true>;
+
+template<typename Impl>
+static void BM_GridIntersection(benchmark::State& state)
+{
+    auto nCells = state.range(0);
+
+    auto nx = static_cast<int>(std::ceil(std::sqrt(nCells)));
+    auto ny = static_cast<int>(std::ceil(std::sqrt(nCells)));
+
+    CoordinateXY center;
+    Envelope env(0, nx, 0, ny);
+    env.centre(center);
+
+    auto geom = geos::benchmark::createSineStar(center, env.getWidth() / 2, 500);
+
+    for (auto _ : state) {
+        Impl::Intersection(env, nx, ny, *geom);
+    }
+
+}
+
+BENCHMARK_TEMPLATE(BM_GridIntersection, GridIntersectionAreaOnly)->Range(1000, 1000000);
+BENCHMARK_TEMPLATE(BM_GridIntersection, GridIntersectionFull)->Range(1000, 1000000);
+BENCHMARK_TEMPLATE(BM_GridIntersection, RectangleIntersection)->Range(1000, 1000000);
+BENCHMARK_TEMPLATE(BM_GridIntersection, PolygonIntersection)->Range(1000, 1000000);
+
+BENCHMARK_MAIN();

capi/geos_c.cpp

@@ -713,6 +713,19 @@ extern "C" {
         return GEOSClipByRect_r(handle, g, xmin, ymin, xmax, ymax);
     }
 
+    Geometry*
+    GEOSSubdivideByGrid(const Geometry* g, double xmin, double ymin, double xmax, double ymax,
+                        unsigned nx, unsigned ny, int include_exterior)
+    {
+        return GEOSSubdivideByGrid_r(handle, g, xmin, ymin, xmax, ymax, nx, ny, include_exterior);
+    }
+
+    int
+    GEOSGridIntersectionFractions(const Geometry* g, double xmin, double ymin, double xmax, double ymax,
+                                  unsigned nx, unsigned ny, float* buf)
+    {
+        return GEOSGridIntersectionFractions_r(handle, g, xmin, ymin, xmax, ymax, nx, ny, buf);
+    }
 
     Geometry*
     GEOSGeom_transformXY(const GEOSGeometry* g, GEOSTransformXYCallback callback, void* userdata) {

capi/geos_c.h.in

@@ -1035,6 +1035,24 @@ extern GEOSGeometry GEOS_DLL *GEOSClipByRect_r(
     double xmin, double ymin,
     double xmax, double ymax);
 
+/** \see GEOSSubdivideByGrid */
+extern GEOSGeometry GEOS_DLL *GEOSSubdivideByGrid_r(
+    GEOSContextHandle_t handle,
+    const GEOSGeometry* g,
+    double xmin, double ymin,
+    double xmax, double ymax,
+    unsigned nx, unsigned ny,
+    int include_exterior);
+
+/** \see GEOSGridIntersectionFractions */
+extern int GEOS_DLL GEOSGridIntersectionFractions_r(
+    GEOSContextHandle_t handle,
+    const GEOSGeometry* g,
+    double xmin, double ymin,
+    double xmax, double ymax,
+    unsigned nx, unsigned ny,
+    float* buf);
+
 /** \see GEOSPolygonize */
 extern GEOSGeometry GEOS_DLL *GEOSPolygonize_r(
     GEOSContextHandle_t handle,
@@ -3788,6 +3806,56 @@ extern GEOSGeometry GEOS_DLL *GEOSClipByRect(
     double xmin, double ymin,
     double xmax, double ymax);
 
+/**
+* Compute the intersection of a geometry with each polygon in
+* a rectangular grid.
+* \param g The input geometry to be clipped
+* \param xmin Left bound of grd
+* \param ymin Lower bound of grid
+* \param xmax Right bound of grid
+* \param ymax Upper bound of grid
+* \param nx number of columns in grid
+* \param ny number of rows in grid
+* \param include_exterior whether to include portions of the
+         input geometry that do not intersect the grid in
+         the returned result.
+* \return A geometry collection whose components represent the
+*         intersection with each cell in the grid.
+* Caller is responsible for freeing with GEOSGeom_destroy().
+* \see GEOSGetGridIntersectionFractions
+*
+* \since 3.14.0
+*/
+extern GEOSGeometry GEOS_DLL *GEOSSubdivideByGrid(
+    const GEOSGeometry* g,
+    double xmin, double ymin,
+    double xmax, double ymax,
+    unsigned nx, unsigned ny,
+    int include_exterior);
+
+/**
+* Determine the fraction of each cell in a rectangular grid
+* that is covered by a polygon.
+* \param g The input geometry
+* \param xmin Left bound of grd
+* \param ymin Lower bound of grid
+* \param xmax Right bound of grid
+* \param ymax Upper bound of grid
+* \param nx number of columns in grid
+* \param ny number of rows in grid
+* \param buf a buffer of size nx*ny to which the grid cell
+*        coverage fractions will be written in row-major order
+* \return 1 if the operation was successful, 0 on exception
+*
+* \since 3.14.0
+*/
+extern int GEOS_DLL GEOSGridIntersectionFractions(
+    const GEOSGeometry* g,
+    double xmin, double ymin,
+    double xmax, double ymax,
+    unsigned nx, unsigned ny,
+    float* buf);
+
 /**
 * Find paths shared between the two given lineal geometries.
 *

capi/geos_ts_c.cpp

@@ -80,6 +80,8 @@
 #include <geos/operation/cluster/GeometryIntersectsClusterFinder.h>
 #include <geos/operation/distance/DistanceOp.h>
 #include <geos/operation/distance/IndexedFacetDistance.h>
+#include <geos/operation/grid/Grid.h>
+#include <geos/operation/grid/GridIntersection.h>
 #include <geos/operation/linemerge/LineMerger.h>
 #include <geos/operation/intersection/Rectangle.h>
 #include <geos/operation/intersection/RectangleIntersection.h>
@@ -1725,6 +1727,41 @@ extern "C" {
         });
     }
 
+    Geometry*
+    GEOSSubdivideByGrid_r(GEOSContextHandle_t extHandle, const Geometry* g, double xmin, double ymin,
+                          double xmax, double ymax, unsigned nx, unsigned ny, int include_exterior)
+    {
+        return execute(extHandle, [&]() {
+            Envelope env(xmin, xmax, ymin, ymax);
+            double dx = env.getWidth() / static_cast<double>(nx);
+            double dy = env.getHeight() / static_cast<double>(ny);
+            geos::operation::grid::Grid<geos::operation::grid::bounded_extent> grid(env, dx, dy);
+
+            return geos::operation::grid::GridIntersection::subdivide_polygon(grid, *g, include_exterior).release();
+        });
+    }
+
+    int
+    GEOSGridIntersectionFractions_r(GEOSContextHandle_t extHandle, const Geometry* g, double xmin, double ymin,
+                                    double xmax, double ymax, unsigned nx, unsigned ny, float* buf)
+    {
+        return execute(extHandle, 0, [&]() {
+            Envelope env(xmin, xmax, ymin, ymax);
+            double dx = env.getWidth() / static_cast<double>(nx);
+            double dy = env.getHeight() / static_cast<double>(ny);
+            geos::operation::grid::Grid<geos::operation::grid::bounded_extent> grid(env, dx, dy);
+
+            // Matrix wants a shared_ptr, but we don't actually want anything to be freed because
+            // buf is externally owned. So we give it an empty deleter.
+            std::shared_ptr<float[]> bufPtr(buf, [](float*) {});
+
+            auto cov = std::make_shared<geos::operation::grid::Matrix<float>>(ny, nx, bufPtr);
+            geos::operation::grid::GridIntersection isect(grid, *g, cov);
+
+            return 1;
+        });
+    }
+
     Geometry*
     GEOSGeom_transformXY_r(GEOSContextHandle_t handle, const GEOSGeometry* g, GEOSTransformXYCallback callback, void* userdata) {
 

include/geos/operation/grid/Cell.h

@@ -0,0 +1,107 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2018-2025 ISciences, LLC
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+
+#include <geos/operation/grid/Side.h>
+#include <geos/operation/grid/Traversal.h>
+
+namespace geos::operation::grid {
+
+/**
+ * @brief The Cell class stores information about the spatial extent of a `Grid` cell and
+ *        any cases where a line crosses that cell (recorded in a `Traversal`).
+ */
+class Cell
+{
+
+  public:
+    Cell(double xmin, double ymin, double xmax, double ymax)
+      : m_box{ xmin, ymin, xmax, ymax }
+    {
+    }
+
+    explicit Cell(const geom::Envelope& b)
+      : m_box{ b }
+    {
+    }
+
+    const geom::Envelope& box() const { return m_box; }
+
+    double width() const;
+
+    double height() const;
+
+    double area() const;
+
+    /// Force the last Coordinate processed (via `take`) to be considered as an
+    /// exit point, provided that it lies on the boundary of this Cell.
+    void force_exit();
+
+    /// Returns whether the cell can be determined to be wholly or partially
+    /// covered by a polygon.
+    bool determined() const;
+
+    /// Return the total length of a linear geometry within this Cell
+    double traversal_length() const;
+
+    /// Return the fraction of this Cell that is covered by a polygon
+    double covered_fraction() const;
+
+    /// Return a newly constructed geometry representing the portion of this Cell
+    /// that is covered by a polygon
+    std::unique_ptr<geom::Geometry> covered_polygons(const geom::GeometryFactory&) const;
+
+    /// Return the last (most recent) traversal to which a coordinate has been
+    /// added. The traversal may or may not be completed.
+    Traversal& last_traversal();
+
+    /**
+     * Attempt to take a coordinate and add it to a Traversal in progress, or start a new Traversal
+     *
+     * @param c             Coordinate to process
+     * @param prev_original The last *uninterpolated* coordinate preceding `c` in the
+     *                      boundary being processed
+     *
+     * @return `true` if the Coordinate was inside this cell, `false` otherwise
+     */
+    bool take(const geom::CoordinateXY& c, const geom::CoordinateXY* prev_original = nullptr);
+
+  private:
+    std::vector<const std::vector<geom::CoordinateXY>*> get_coord_lists() const;
+
+    enum class Location
+    {
+        INSIDE,
+        OUTSIDE,
+        BOUNDARY
+    };
+
+    geom::Envelope m_box;
+
+    std::vector<Traversal> m_traversals;
+
+    Side side(const geom::CoordinateXY& c) const;
+
+    Location location(const geom::CoordinateXY& c) const;
+
+    /// If no Traversals have been started or the most recent Traversal has been completed,
+    /// return a new Traversal. Otherwise, return the most recent Traversal.
+    Traversal& traversal_in_progress();
+};
+
+}