diff --git a/capi/geos_c.cpp b/capi/geos_c.cpp
index 9bd1068ba4..fb0a3f0c0a 100644
--- a/capi/geos_c.cpp
+++ b/capi/geos_c.cpp
@@ -2061,5 +2061,15 @@ extern "C" {
         return GEOSCoverageSimplifyVW_r(handle, input, tolerance, preserveBoundary);
     }
 
+    Geometry*
+    GEOSCoverageSimplifyVWWithProgress(const Geometry* input, double tolerance,
+                                       int preserveBoundary,
+                                       GEOSProgressCallback_r progressFunc,
+                                       void* progressUserData)
+    {
+        return GEOSCoverageSimplifyVWWithProgress_r(
+            handle, input, tolerance, preserveBoundary, progressFunc, progressUserData);
+    }
+
 
 } /* extern "C" */
diff --git a/capi/geos_c.h.in b/capi/geos_c.h.in
index 39671cc696..e488af7f6c 100644
--- a/capi/geos_c.h.in
+++ b/capi/geos_c.h.in
@@ -117,6 +117,21 @@ typedef void (*GEOSMessageHandler)(GEOS_PRINTF_FORMAT const char *fmt, ...)
 */
 typedef void (*GEOSMessageHandler_r)(const char *message, void *userdata);
 
+
+/**
+ * A GEOS progression callback function.
+ *
+ * Such function takes a progression ratio, and an optional message.
+ *
+ * \param progressRatio Progression ratio (between 0 and 1)
+ * \param message Information message (can be NULL)
+ * \param userdata the user data pointer that was passed to GEOS together with
+ * this callback.
+ */
+typedef void (*GEOSProgressCallback_r)(double progressRatio,
+                                       const char* message, void* userdata);
+
+
 /*
 * When we're included by geos_c.cpp, these types are #defined to the
 * C++ definitions via preprocessor. We don't touch them to allow the
@@ -886,6 +901,16 @@ GEOSCoverageSimplifyVW_r(
     double tolerance,
     int preserveBoundary);
 
+/** \see GEOSCoverageSimplifyVWWithProgress */
+extern GEOSGeometry GEOS_DLL *
+GEOSCoverageSimplifyVWWithProgress_r(
+    GEOSContextHandle_t extHandle,
+    const GEOSGeometry* input,
+    double tolerance,
+    int preserveBoundary,
+    GEOSProgressCallback_r progressFunc,
+    void* progressUserData);
+
 /** \see GEOSCoverageCleanParams_create */
 extern GEOSCoverageCleanParams GEOS_DLL *
 GEOSCoverageCleanParams_create_r(
@@ -1076,6 +1101,13 @@ extern GEOSGeometry GEOS_DLL *GEOSUnaryUnion_r(
     GEOSContextHandle_t handle,
     const GEOSGeometry* g);
 
+/** \see GEOSUnaryUnion */
+extern GEOSGeometry GEOS_DLL *GEOSUnaryUnionWithProgress_r(
+    GEOSContextHandle_t handle,
+    const GEOSGeometry* g,
+    GEOSProgressCallback_r progressFunc,
+    void* progressUserData);
+
 /** \see GEOSUnaryUnionPrec */
 extern GEOSGeometry GEOS_DLL *GEOSUnaryUnionPrec_r(
     GEOSContextHandle_t handle,
@@ -4321,6 +4353,40 @@ extern GEOSGeometry GEOS_DLL * GEOSCoverageSimplifyVW(
     double tolerance,
     int preserveBoundary);
 
+/**
+* Operates on a coverage (represented as a list of polygonal geometry
+* with exactly matching edge geometry) to apply a Visvalingam–Whyatt
+* simplification to the edges, reducing complexity in proportion with
+* the provided tolerance, while retaining a valid coverage (no edges
+* will cross or touch after the simplification).
+* Geometries never disappear, but they may be simplified down to just
+* a triangle. Also, some invalid geoms (such as Polygons which have too
+* few non-repeated points) will be returned unchanged.
+* If the input dataset is not a valid coverage due to overlaps,
+* it will still be simplified, but invalid topology such as crossing
+* edges will still be invalid.
+*
+* \param input The polygonal coverage to access,
+*        stored in a geometry collection. All members must be POLYGON
+*        or MULTIPOLYGON.
+* \param tolerance A tolerance parameter in linear units.
+* \param preserveBoundary Use 1 to preserve the outside edges
+*        of the coverage without simplification,
+*        0 to allow them to be simplified.
+* \param progressFunc Progress callback (or null)
+* \param progressUserData User data passed to progress callback (can be null)
+* \return A collection containing the simplified geometries, or null
+*         on error.
+*
+* \since 3.14
+*/
+extern GEOSGeometry GEOS_DLL * GEOSCoverageSimplifyVWWithProgress(
+    const GEOSGeometry* input,
+    double tolerance,
+    int preserveBoundary,
+    GEOSProgressCallback_r progressFunc,
+    void* progressUserData);
+
 /**
 * Create a default GEOSCoverageCleanParams object for controlling
 * the way invalid polygon interactions are repaird by \ref GEOSCoverageCleanWithParams.
@@ -5328,7 +5394,7 @@ extern char GEOS_DLL GEOSIntersects(const GEOSGeometry* g1, const GEOSGeometry*
 extern char GEOS_DLL GEOSCrosses(const GEOSGeometry* g1, const GEOSGeometry* g2);
 
 /**
-* Tests if geometry g1 is completely within g2, 
+* Tests if geometry g1 is completely within g2,
 * but not wholly contained in the boundary of g2.
 * \param g1 Input geometry
 * \param g2 Input geometry
diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp
index 27ee5b8123..06131a6cce 100644
--- a/capi/geos_ts_c.cpp
+++ b/capi/geos_ts_c.cpp
@@ -114,6 +114,7 @@
 #include <geos/util/Interrupt.h>
 #include <geos/util/UniqueCoordinateArrayFilter.h>
 #include <geos/util/Machine.h>
+#include <geos/util/Progress.h>
 #include <geos/version.h>
 
 // This should go away
@@ -1666,6 +1667,24 @@ extern "C" {
         });
     }
 
+    Geometry*
+    GEOSUnaryUnionWithProgress_r(GEOSContextHandle_t extHandle, const Geometry* g,
+                                 GEOSProgressCallback_r progressFunc,
+                                 void* progressUserData)
+    {
+        geos::util::ProgressFunction progressFunction(
+            [progressFunc, progressUserData](double progress, const char* message)
+        {
+            progressFunc(progress, message, progressUserData);
+        });
+
+        return execute(extHandle, [&]() {
+            std::unique_ptr<Geometry> g3(g->Union(&progressFunction));
+            g3->setSRID(g->getSRID());
+            return g3.release();
+        });
+    }
+
     Geometry*
     GEOSUnaryUnionPrec_r(GEOSContextHandle_t extHandle, const Geometry* g1, double gridSize)
     {
@@ -4522,6 +4541,25 @@ extern "C" {
         double tolerance,
         int preserveBoundary)
     {
+        return GEOSCoverageSimplifyVWWithProgress_r(extHandle, input, tolerance,
+                                                    preserveBoundary,
+                                                    nullptr, nullptr);
+    }
+
+    Geometry*
+    GEOSCoverageSimplifyVWWithProgress_r(GEOSContextHandle_t extHandle,
+        const Geometry* input,
+        double tolerance,
+        int preserveBoundary,
+        GEOSProgressCallback_r progressFunc,
+        void* progressUserData)
+    {
+        geos::util::ProgressFunction progressFunction(
+            [progressFunc, progressUserData](double progress, const char* message)
+        {
+            progressFunc(progress, message, progressUserData);
+        });
+
         using geos::coverage::CoverageSimplifier;
 
         return execute(extHandle, [&]() -> Geometry* {
@@ -4536,10 +4574,10 @@ extern "C" {
             CoverageSimplifier cov(coverage);
             std::vector<std::unique_ptr<Geometry>> simple;
             if (preserveBoundary == 1) {
-                simple = cov.simplifyInner(tolerance);
+                simple = cov.simplifyInner(tolerance, progressFunc ? progressFunction : geos::util::defaultProgress);
             }
             else if (preserveBoundary == 0) {
-                simple = cov.simplify(tolerance);
+                simple = cov.simplify(tolerance, progressFunc ? progressFunction : geos::util::defaultProgress);
             }
             else return nullptr;
 
diff --git a/include/geos/coverage/CoverageBoundarySegmentFinder.h b/include/geos/coverage/CoverageBoundarySegmentFinder.h
index 359ea7aecb..3639561100 100644
--- a/include/geos/coverage/CoverageBoundarySegmentFinder.h
+++ b/include/geos/coverage/CoverageBoundarySegmentFinder.h
@@ -20,6 +20,7 @@
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/LineSegment.h>
 #include <geos/export.h>
+#include <geos/util/Progress.h>
 
 namespace geos {
 namespace geom {
@@ -57,7 +58,9 @@ class CoverageBoundarySegmentFinder : public geos::geom::CoordinateSequenceFilte
 
 
     static LineSegment::UnorderedSet
-        findBoundarySegments(const std::vector<const Geometry*>& geoms);
+        findBoundarySegments(const std::vector<const Geometry*>& geoms,
+                             const util::ProgressFunction& progressFunction = util::defaultProgress);
+
 
     static bool isBoundarySegment(
         const LineSegment::UnorderedSet& boundarySegs,
diff --git a/include/geos/coverage/CoverageEdge.h b/include/geos/coverage/CoverageEdge.h
index 60afd1c7ca..b695648396 100644
--- a/include/geos/coverage/CoverageEdge.h
+++ b/include/geos/coverage/CoverageEdge.h
@@ -20,6 +20,7 @@
 #include <geos/geom/CoordinateSequence.h>
 #include <geos/geom/LineSegment.h>
 #include <geos/util.h>
+#include <geos/util/Progress.h>
 
 // Forward declarations
 namespace geos {
@@ -116,7 +117,8 @@ class GEOS_DLL CoverageEdge {
 
     static std::unique_ptr<MultiLineString> createLines(
         const std::vector<CoverageEdge*>& edges,
-        const GeometryFactory* geomFactory);
+        const GeometryFactory* geomFactory,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     std::unique_ptr<LineString> toLineString(
         const GeometryFactory* geomFactory);
diff --git a/include/geos/coverage/CoverageRingEdges.h b/include/geos/coverage/CoverageRingEdges.h
index 5a0f4ac1b9..6709672a22 100644
--- a/include/geos/coverage/CoverageRingEdges.h
+++ b/include/geos/coverage/CoverageRingEdges.h
@@ -18,6 +18,7 @@
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/LineSegment.h>
 #include <geos/coverage/CoverageEdge.h> // to materialize CoverageEdge
+#include <geos/util/Progress.h>
 
 #include <set>
 #include <map>
@@ -72,10 +73,11 @@ class GEOS_DLL CoverageRingEdges {
 
 public:
 
-    CoverageRingEdges(const std::vector<const Geometry*>& coverage)
+    CoverageRingEdges(const std::vector<const Geometry*>& coverage,
+                      const util::ProgressFunction& progressFunction = util::defaultProgress)
         : m_coverage(coverage)
     {
-        build();
+        build(progressFunction);
     };
 
 
@@ -96,14 +98,15 @@ class GEOS_DLL CoverageRingEdges {
     /**
     * Recreates the polygon coverage from the current edge values.
     *
+    * @param progressFunction Progress function or null
     * @return an array of polygonal geometries representing the coverage
     */
-    std::vector<std::unique_ptr<Geometry>> buildCoverage() const;
+    std::vector<std::unique_ptr<Geometry>> buildCoverage(const util::ProgressFunction& progressFunction = util::defaultProgress) const;
 
 
 private:
 
-    void build();
+    void build(const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     void addRingEdges(
         const LinearRing* ring,
@@ -140,10 +143,12 @@ class GEOS_DLL CoverageRingEdges {
         const CoordinateSequence& ring);
 
     Coordinate::UnorderedSet findMultiRingNodes(
-        const std::vector<const Geometry*>& coverage);
+        const std::vector<const Geometry*>& coverage,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     Coordinate::UnorderedSet findBoundaryNodes(
-        LineSegment::UnorderedSet& lineSegments);
+        LineSegment::UnorderedSet& lineSegments,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     std::unique_ptr<Geometry> buildPolygonal(
         const Geometry* geom) const;
diff --git a/include/geos/coverage/CoverageSimplifier.h b/include/geos/coverage/CoverageSimplifier.h
index 74df234dda..bf788f283f 100644
--- a/include/geos/coverage/CoverageSimplifier.h
+++ b/include/geos/coverage/CoverageSimplifier.h
@@ -19,6 +19,7 @@
 #include <vector>
 #include <memory>
 #include <geos/export.h>
+#include <geos/util/Progress.h>
 
 
 namespace geos {
@@ -87,15 +88,18 @@ class GEOS_DLL CoverageSimplifier {
     *
     * @param coverage a set of polygonal geometries forming a coverage
     * @param tolerance the simplification tolerance
+    * @param progressFunction Progress function
     * @return the simplified polygons
     */
     static std::vector<std::unique_ptr<Geometry>> simplify(
         std::vector<const Geometry*>& coverage,
-        double tolerance);
+        double tolerance,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     static std::vector<std::unique_ptr<Geometry>> simplify(
         const std::vector<std::unique_ptr<Geometry>>& coverage,
-        double tolerance);
+        double tolerance,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     /**
     * Simplifies the inner boundaries of a set of polygonal geometries forming a coverage,
@@ -104,24 +108,28 @@ class GEOS_DLL CoverageSimplifier {
     *
     * @param coverage a set of polygonal geometries forming a coverage
     * @param tolerance the simplification tolerance
+    * @param progressFunction Progress function
     * @return the simplified polygons
     */
     static std::vector<std::unique_ptr<Geometry>> simplifyInner(
         std::vector<const Geometry*>& coverage,
-        double tolerance);
+        double tolerance,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     static std::vector<std::unique_ptr<Geometry>> simplifyInner(
         const std::vector<std::unique_ptr<Geometry>>& coverage,
-        double tolerance);
+        double tolerance,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     /**
     * Computes the simplified coverage, preserving the coverage topology.
     *
     * @param tolerance the simplification tolerance
+    * @param progressFunction Progress function
     * @return the simplified polygons
     */
     std::vector<std::unique_ptr<Geometry>> simplify(
-        double tolerance);
+        double tolerance, const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     /**
     * Computes the inner-boundary simplified coverage,
@@ -129,10 +137,11 @@ class GEOS_DLL CoverageSimplifier {
     * and leaving outer boundary edges unchanged.
     *
     * @param tolerance the simplification tolerance
+    * @param progressFunction Progress function
     * @return the simplified polygons
     */
     std::vector<std::unique_ptr<Geometry>> simplifyInner(
-        double tolerance);
+        double tolerance, const util::ProgressFunction& progressFunction = util::defaultProgress);
 
 
 private:
@@ -145,7 +154,8 @@ class GEOS_DLL CoverageSimplifier {
     void simplifyEdges(
         std::vector<CoverageEdge*> edges,
         const MultiLineString* constraints,
-        double tolerance);
+        double tolerance,
+        const util::ProgressFunction& progressFunction);
 
     void setCoordinates(
         std::vector<CoverageEdge*>& edges,
diff --git a/include/geos/coverage/TPVWSimplifier.h b/include/geos/coverage/TPVWSimplifier.h
index 0b904c8a83..fbe3cc5095 100644
--- a/include/geos/coverage/TPVWSimplifier.h
+++ b/include/geos/coverage/TPVWSimplifier.h
@@ -21,6 +21,7 @@
 #include <geos/simplify/LinkedLine.h>
 #include <geos/coverage/Corner.h>
 #include <geos/export.h>
+#include <geos/util/Progress.h>
 
 
 namespace geos {
@@ -164,11 +165,13 @@ class GEOS_DLL TPVWSimplifier
     *
     * @param lines the lines to simplify
     * @param distanceTolerance the simplification tolerance
+    * @param progressFunction Progress function
     * @return the simplified lines
     */
     static std::unique_ptr<MultiLineString> simplify(
         const MultiLineString* lines,
-        double distanceTolerance);
+        double distanceTolerance,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     /**
     * Simplifies a set of lines, preserving the topology of the lines between
@@ -181,13 +184,15 @@ class GEOS_DLL TPVWSimplifier
     * @param freeRings flags indicating which ring edges do not have node endpoints
     * @param constraintLines the linear constraints
     * @param distanceTolerance the simplification tolerance
+    * @param progressFunction Progress function
     * @return the simplified lines
     */
     static std::unique_ptr<MultiLineString> simplify(
         const MultiLineString* lines,
         std::vector<bool>& freeRings,
         const MultiLineString* constraintLines,
-        double distanceTolerance);
+        double distanceTolerance,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     // Constructor
     TPVWSimplifier(const MultiLineString* lines,
@@ -209,11 +214,12 @@ class GEOS_DLL TPVWSimplifier
 
     void setConstraints(const MultiLineString* constraints);
 
-    std::unique_ptr<MultiLineString> simplify();
+    std::unique_ptr<MultiLineString> simplify(const util::ProgressFunction& progressFunction = util::defaultProgress);
 
     std::vector<Edge> createEdges(
         const MultiLineString* lines,
-        std::vector<bool>& freeRing);
+        std::vector<bool>& freeRing,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
 
 }; // TPVWSimplifier
diff --git a/include/geos/coverage/VertexRingCounter.h b/include/geos/coverage/VertexRingCounter.h
index c0982ad2b4..5cc242ca70 100644
--- a/include/geos/coverage/VertexRingCounter.h
+++ b/include/geos/coverage/VertexRingCounter.h
@@ -19,6 +19,7 @@
 #include <geos/geom/CoordinateSequenceFilter.h>
 #include <geos/geom/Coordinate.h>
 #include <geos/export.h>
+#include <geos/util/Progress.h>
 
 namespace geos {
 namespace geom {
@@ -54,7 +55,8 @@ class VertexRingCounter : public geos::geom::CoordinateSequenceFilter
 
     static void count(
         const std::vector<const Geometry*>& geoms,
-        std::map<Coordinate, std::size_t>& counts);
+        std::map<Coordinate, std::size_t>& counts,
+        const util::ProgressFunction& progressFunction = util::defaultProgress);
 
 private:
 
diff --git a/include/geos/geom/Geometry.h b/include/geos/geom/Geometry.h
index 36e64f0421..4325767a60 100644
--- a/include/geos/geom/Geometry.h
+++ b/include/geos/geom/Geometry.h
@@ -35,6 +35,7 @@
 #include <geos/geom/Dimension.h> // for Dimension::DimensionType
 #include <geos/geom/GeometryComponentFilter.h> // for inheritance
 #include <geos/geom/CoordinateSequence.h> // to materialize CoordinateSequence
+#include <geos/util/Progress.h>
 
 #include <algorithm>
 #include <string>
@@ -726,9 +727,15 @@ class GEOS_DLL Geometry {
      *
      * @see UnaryUnionOp
      */
-    Ptr Union() const;
+    Ptr Union(geos::util::ProgressFunction* progressFunction) const;
     // throw(IllegalArgumentException *, TopologyException *);
 
+    Ptr Union() const
+    {
+        geos::util::ProgressFunction* progressFunction = nullptr;
+        return Union(progressFunction);
+    }
+
     /**
      * \brief
      * Returns a Geometry representing the points making up this
diff --git a/include/geos/operation/union/CascadedPolygonUnion.h b/include/geos/operation/union/CascadedPolygonUnion.h
index 2b1e292670..fcc5ba9a76 100644
--- a/include/geos/operation/union/CascadedPolygonUnion.h
+++ b/include/geos/operation/union/CascadedPolygonUnion.h
@@ -25,6 +25,7 @@
 #include <geos/export.h>
 
 #include <geos/operation/union/UnionStrategy.h>
+#include <geos/util/Progress.h>
 
 // Forward declarations
 namespace geos {
@@ -140,7 +141,7 @@ class GEOS_DLL CascadedPolygonUnion {
      *              ownership of elements *and* vector are left to caller.
      */
     static std::unique_ptr<geom::Geometry> Union(std::vector<geom::Polygon*>* polys);
-    static std::unique_ptr<geom::Geometry> Union(std::vector<geom::Polygon*>* polys, UnionStrategy* unionFun);
+    static std::unique_ptr<geom::Geometry> Union(std::vector<geom::Polygon*>* polys, UnionStrategy* unionFun, geos::util::ProgressFunction* progressFunction);
 
     /** \brief
      * Computes the union of a set of polygonal [Geometrys](@ref geom::Geometry).
@@ -149,17 +150,18 @@ class GEOS_DLL CascadedPolygonUnion {
      * @param start start iterator
      * @param end end iterator
      * @param unionStrategy strategy to apply
+     * @param progressFunction progress function
      */
     template <class T>
     static std::unique_ptr<geom::Geometry>
-    Union(T start, T end, UnionStrategy *unionStrategy)
+    Union(T start, T end, UnionStrategy *unionStrategy, geos::util::ProgressFunction* progressFunction)
     {
         std::vector<geom::Polygon*> polys;
         for(T i = start; i != end; ++i) {
             const geom::Polygon* p = dynamic_cast<const geom::Polygon*>(*i);
             polys.push_back(const_cast<geom::Polygon*>(p));
         }
-        return Union(&polys, unionStrategy);
+        return Union(&polys, unionStrategy, progressFunction);
     }
 
     /** \brief
@@ -167,8 +169,9 @@ class GEOS_DLL CascadedPolygonUnion {
      *
      * @param polys a collection of polygonal [Geometrys](@ref geom::Geometry).
      *              Ownership of elements *and* vector are left to caller.
+     * @param progressFunction progress function
      */
-    static std::unique_ptr<geom::Geometry> Union(const geom::MultiPolygon* polys);
+    static std::unique_ptr<geom::Geometry> Union(const geom::MultiPolygon* polys, geos::util::ProgressFunction* progressFunction);
 
     /** \brief
      * Creates a new instance to union the given collection of
@@ -192,10 +195,11 @@ class GEOS_DLL CascadedPolygonUnion {
     /** \brief
      * Computes the union of the input geometries.
      *
+     * @param progressFunction progress function
      * @return the union of the input geometries
      * @return `null` if no input geometries were provided
      */
-    std::unique_ptr<geom::Geometry> Union();
+    std::unique_ptr<geom::Geometry> Union(geos::util::ProgressFunction* progressFunction);
 
 private:
 
@@ -211,7 +215,11 @@ class GEOS_DLL CascadedPolygonUnion {
      * @param end the index after the end of the section
      * @return the union of the list section
      */
-    std::unique_ptr<geom::Geometry> binaryUnion(const std::vector<const geom::Geometry*> & geoms, std::size_t start, std::size_t end);
+    std::unique_ptr<geom::Geometry> binaryUnion(
+        const std::vector<const geom::Geometry*> & geoms,
+        std::size_t start,
+        std::size_t end,
+        std::function<void()>* unitProgress);
 
     /**
      * Computes the union of two geometries,
diff --git a/include/geos/operation/union/UnaryUnionOp.h b/include/geos/operation/union/UnaryUnionOp.h
index 3fed7e32bc..98b6a9e700 100644
--- a/include/geos/operation/union/UnaryUnionOp.h
+++ b/include/geos/operation/union/UnaryUnionOp.h
@@ -28,6 +28,7 @@
 #include <geos/geom/Polygon.h>
 #include <geos/geom/util/GeometryExtracter.h>
 #include <geos/operation/union/CascadedPolygonUnion.h>
+#include <geos/util/Progress.h>
 
 #include <geos/util.h>
 
@@ -93,7 +94,8 @@ class GEOS_DLL UnaryUnionOp {
     Union(const T& geoms)
     {
         UnaryUnionOp op(geoms);
-        return op.Union();
+        geos::util::ProgressFunction* progressFunction = nullptr;
+        return op.Union(progressFunction);
     }
 
     template <class T>
@@ -102,14 +104,15 @@ class GEOS_DLL UnaryUnionOp {
           geom::GeometryFactory& geomFact)
     {
         UnaryUnionOp op(geoms, geomFact);
-        return op.Union();
+        geos::util::ProgressFunction* progressFunction = nullptr;
+        return op.Union(progressFunction);
     }
 
     static std::unique_ptr<geom::Geometry>
-    Union(const geom::Geometry& geom)
+    Union(const geom::Geometry& geom, geos::util::ProgressFunction* progressFunction)
     {
         UnaryUnionOp op(geom);
-        return op.Union();
+        return op.Union(progressFunction);
     }
 
     template <class T>
@@ -150,7 +153,7 @@ class GEOS_DLL UnaryUnionOp {
      * @return an empty GEOMETRYCOLLECTION if no geometries were provided
      *         in the input
      */
-    std::unique_ptr<geom::Geometry> Union();
+    std::unique_ptr<geom::Geometry> Union(geos::util::ProgressFunction* progressFunction);
 
 private:
 
diff --git a/include/geos/util/Progress.h b/include/geos/util/Progress.h
new file mode 100644
index 0000000000..34ff622d04
--- /dev/null
+++ b/include/geos/util/Progress.h
@@ -0,0 +1,121 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2025 Even Rouault <even.rouault@spatialys.com>
+ *
+ * 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 <functional>
+#include <optional>
+
+#include <geos/export.h>
+
+namespace geos::util {
+
+/** A ProgressFunction is a wrapper around an optional user-defined callback,
+ * taking a progress ratio (between 0 and 1) and an optional message. If not
+ * provided with a callback at construction, the ProgressFunction will do nothing
+ * when invoked.
+ */
+class GEOS_DLL ProgressFunction {
+public:
+    ProgressFunction() : m_function(std::nullopt) {}
+
+    ProgressFunction(std::function<void(double, const char*)> f) : m_function(f) {}
+
+    void operator()(double percentage, const char* message) const {
+        if (m_function.has_value()) {
+            m_function.value()(percentage, message);
+        }
+    }
+
+   /** Create a ProgressFunction to manage a subset of the work reported
+    * by this ProgressFunction.
+    *
+    * Sometimes when an operation wants to report progress, it actually
+    * invokes several subprocesses which also take a ProgressFunction,
+    * and it is desirable to map the progress of each sub operation into
+    * a portion of 0.0 to 1.0 progress of the overall process. The scaled
+    * progress function can be used for this.
+    *
+    * For each subsection a scaled progress function is created and
+    * instead of passing the overall progress func down to the sub functions,
+    * the scaled progress function is passed instead.
+    *
+    * @param from the completion fraction (0 to 1) to which 0.0 in the sub operation is mapped.
+    * @param to the completion fraction (0 to 1) to which 1.0 is the sub operation is mapped.
+    *
+    * @return scaled progress function.
+    */
+    ProgressFunction subProgress(double from, double to) const;
+
+    bool isSpecified() const {
+        return m_function.has_value();
+    }
+
+private:
+    std::optional<std::function<void(double, const char*)>> m_function;
+};
+
+static const ProgressFunction defaultProgress;
+
+/** A ProgressContext manages the invocation of a ProgressFunction at a specified frequency.
+ */
+class GEOS_DLL ProgressContext {
+
+public:
+    /** Create a ProgressContext
+     *
+     * @param pCallback the function to call with progress updates
+     * @param pIterCount the total number of expected iterations
+     */
+    ProgressContext(ProgressFunction pCallback, size_t pIterCount) :
+        callback(pCallback),
+        i(0),
+        iterCount(pIterCount),
+        notificationInterval(std::max<size_t>(1, iterCount / 100)),
+        iNotify(0) {}
+
+    /** Set the resolution of the progress reporting as a fraction from 0 to 1. By default,
+     *  the progress function will be called for each 1% change in progress.
+     */
+    void setResolution(double res) {
+        notificationInterval = std::max<size_t>(1, static_cast<std::size_t>(static_cast<double>(iterCount) * res));
+    }
+
+    /** Update the progress. This method should be called once per iteration. */
+    void update() {
+        if (callback.isSpecified()) {
+            if (iNotify + 1 == notificationInterval) {
+                callback(static_cast<double>(i + 1)/static_cast<double>(iterCount), nullptr);
+                iNotify = 0;
+            }
+            else {
+                ++iNotify;
+            }
+        }
+    }
+
+    void finish() const {
+        callback(1.0, nullptr);
+    }
+
+private:
+    const ProgressFunction callback;
+    std::size_t i;
+    std::size_t iterCount;
+    std::size_t notificationInterval;
+    std::size_t iNotify;
+};
+
+} // namespace geos::util
+
diff --git a/src/coverage/CoverageBoundarySegmentFinder.cpp b/src/coverage/CoverageBoundarySegmentFinder.cpp
index d5759be7d4..c1fd187346 100644
--- a/src/coverage/CoverageBoundarySegmentFinder.cpp
+++ b/src/coverage/CoverageBoundarySegmentFinder.cpp
@@ -28,12 +28,17 @@ namespace coverage { // geos.coverage
 // public static
 LineSegment::UnorderedSet
 CoverageBoundarySegmentFinder::findBoundarySegments(
-    const std::vector<const Geometry*>& geoms)
+    const std::vector<const Geometry*>& geoms,
+    const util::ProgressFunction& progressFunction)
 {
     LineSegment::UnorderedSet segs;
     CoverageBoundarySegmentFinder finder(segs);
-    for (const Geometry* geom : geoms) {
+
+    util::ProgressContext progress(progressFunction, geoms.size());
+
+    for (const auto& geom : geoms) {
         geom->apply_ro(finder);
+        progress.update();
     }
     return segs;
 }
diff --git a/src/coverage/CoverageEdge.cpp b/src/coverage/CoverageEdge.cpp
index 65971fc3ad..dd5bfe1eff 100644
--- a/src/coverage/CoverageEdge.cpp
+++ b/src/coverage/CoverageEdge.cpp
@@ -55,14 +55,22 @@ CoverageEdge::createEdge(const CoordinateSequence& ring,
 std::unique_ptr<MultiLineString>
 CoverageEdge::createLines(
     const std::vector<CoverageEdge*>& edges,
-    const GeometryFactory* geomFactory)
+    const GeometryFactory* geomFactory,
+    const util::ProgressFunction& progressFunction)
 {
     std::vector<std::unique_ptr<LineString>> lines;
+
+    util::ProgressContext progress(progressFunction, edges.size());
+
     for (const CoverageEdge* edge : edges) {
         auto cs = edge->getCoordinates()->clone();
         auto ls = geomFactory->createLineString(std::move(cs));
         lines.push_back(std::move(ls));
+        progress.update();
     }
+
+    progress.finish();
+
     return geomFactory->createMultiLineString(std::move(lines));
 }
 
diff --git a/src/coverage/CoverageRingEdges.cpp b/src/coverage/CoverageRingEdges.cpp
index 9d5eba07fd..9801f381e5 100644
--- a/src/coverage/CoverageRingEdges.cpp
+++ b/src/coverage/CoverageRingEdges.cpp
@@ -61,14 +61,26 @@ CoverageRingEdges::selectEdges(std::size_t ringCount) const
 
 /* private */
 void
-CoverageRingEdges::build()
+CoverageRingEdges::build(const util::ProgressFunction& progressFunction)
 {
-    Coordinate::UnorderedSet nodes = findMultiRingNodes(m_coverage);
-    LineSegment::UnorderedSet boundarySegs = CoverageBoundarySegmentFinder::findBoundarySegments(m_coverage);
-    Coordinate::UnorderedSet boundaryNodes = findBoundaryNodes(boundarySegs);
+    constexpr double RATIO_FIRST_PASS = 0.1;
+    constexpr double RATIO_SECOND_PASS = 0.2;
+    constexpr double RATIO_THIRD_PASS = 0.9;
+    constexpr double RATIO_FOURTH_PASS = 1.0;
+
+    Coordinate::UnorderedSet nodes = findMultiRingNodes(m_coverage, progressFunction.subProgress(0, RATIO_FIRST_PASS));
+
+    LineSegment::UnorderedSet boundarySegs = CoverageBoundarySegmentFinder::findBoundarySegments(
+        m_coverage, progressFunction.subProgress(RATIO_FIRST_PASS, RATIO_SECOND_PASS));
+
+    Coordinate::UnorderedSet boundaryNodes = findBoundaryNodes(
+        boundarySegs, progressFunction.subProgress(RATIO_SECOND_PASS, RATIO_THIRD_PASS));
     nodes.insert(boundaryNodes.begin(), boundaryNodes.end());
 
     std::map<LineSegment, CoverageEdge*> uniqueEdgeMap;
+
+    util::ProgressContext progress(progressFunction.subProgress(RATIO_THIRD_PASS, RATIO_FOURTH_PASS), m_coverage.size());
+
     for (const Geometry* geom : m_coverage) {
         for (std::size_t ipoly = 0; ipoly < geom->getNumGeometries(); ipoly++) {
             util::ensureNoCurvedComponents(geom->getGeometryN(ipoly));
@@ -91,7 +103,11 @@ CoverageRingEdges::build()
                 addRingEdges(hole, nodes, boundarySegs, uniqueEdgeMap);
             }
         }
+
+        progress.update();
     }
+
+    progress.finish();
 }
 
 /* private */
@@ -254,55 +270,84 @@ CoverageRingEdges::next(std::size_t index, const CoordinateSequence& ring)
 
 /* private */
 Coordinate::UnorderedSet
-CoverageRingEdges::findMultiRingNodes(const std::vector<const Geometry*>& coverage)
+CoverageRingEdges::findMultiRingNodes(const std::vector<const Geometry*>& coverage,
+                                      const util::ProgressFunction& progressFunction)
 {
     std::map<Coordinate, std::size_t> vertexRingCount;
-    VertexRingCounter::count(coverage, vertexRingCount);
+
+    VertexRingCounter::count(coverage, vertexRingCount, progressFunction.subProgress(0, 0.5));
     Coordinate::UnorderedSet nodes;
     // for (Coordinate v : vertexCount.keySet()) {
     //     if (vertexCount.get(v) > 2) {
     //         nodes.add(v);
     //     }
     // }
+
+    util::ProgressContext progress(progressFunction.subProgress(0.5, 1.0), vertexRingCount.size());
+
     for (const auto &mapPair : vertexRingCount) {
         const Coordinate& v = mapPair.first;
         std::size_t count = mapPair.second;
         if (count > 2)
             nodes.insert(v);
+
+        progress.update();
     }
+    progress.finish();
     return nodes;
 }
 
 
 /* private */
 Coordinate::UnorderedSet
-CoverageRingEdges::findBoundaryNodes(LineSegment::UnorderedSet& boundarySegments)
+CoverageRingEdges::findBoundaryNodes(LineSegment::UnorderedSet& boundarySegments,
+                                     const util::ProgressFunction& progressFunction)
 {
     std::unordered_map<Coordinate, std::size_t, Coordinate::HashCode> counter;
-    for (const LineSegment& seg : boundarySegments) {
-        counter[seg.p0] += 1;
-        counter[seg.p1] += 1;
+
+    util::ProgressContext progress1(progressFunction.subProgress(0, 0.5), boundarySegments.size());
+
+    {
+        for (const LineSegment& seg : boundarySegments) {
+            counter[seg.p0] += 1;
+            counter[seg.p1] += 1;
+            progress1.update();
+        }
     }
 
+    util::ProgressContext progress2(progressFunction.subProgress(0.5, 1.0), counter.size());
+
     Coordinate::UnorderedSet nodes;
     for (const auto& c : counter) {
         const Coordinate& v = c.first;
         std::size_t count = c.second;
         if (count > 2)
             nodes.insert(v);
+
+        progress2.update();
     }
+
+    progress2.finish();
+
     return nodes;
 }
 
 
 /* public */
 std::vector<std::unique_ptr<Geometry>>
-CoverageRingEdges::buildCoverage() const
+CoverageRingEdges::buildCoverage(const util::ProgressFunction& progressFunction) const
 {
     std::vector<std::unique_ptr<Geometry>> result;
+
+    util::ProgressContext progress(progressFunction, m_coverage.size());
+
     for (const Geometry* geom : m_coverage) {
         result.push_back(buildPolygonal(geom));
+        progress.update();
     }
+
+    progress.finish();
+
     return result;
 }
 
diff --git a/src/coverage/CoverageSimplifier.cpp b/src/coverage/CoverageSimplifier.cpp
index 98ba1ad877..e80f49c960 100644
--- a/src/coverage/CoverageSimplifier.cpp
+++ b/src/coverage/CoverageSimplifier.cpp
@@ -26,7 +26,7 @@
 using geos::geom::Geometry;
 using geos::geom::GeometryFactory;
 using geos::geom::MultiLineString;
-
+using geos::util::ProgressFunction;
 
 namespace geos {     // geos
 namespace coverage { // geos.coverage
@@ -35,23 +35,25 @@ namespace coverage { // geos.coverage
 std::vector<std::unique_ptr<Geometry>>
 CoverageSimplifier::simplify(
     std::vector<const Geometry*>& coverage,
-    double tolerance)
+    double tolerance,
+    const ProgressFunction& progressFunction)
 {
     CoverageSimplifier simplifier(coverage);
-    return simplifier.simplify(tolerance);
+    return simplifier.simplify(tolerance, progressFunction);
 }
 
 /* public static */
 std::vector<std::unique_ptr<Geometry>>
 CoverageSimplifier::simplify(
     const std::vector<std::unique_ptr<Geometry>>& coverage,
-    double tolerance)
+    double tolerance,
+    const ProgressFunction& progressFunction)
 {
     std::vector<const Geometry*> geoms;
     for (auto& geom : coverage) {
         geoms.push_back(geom.get());
     }
-    return simplify(geoms, tolerance);
+    return simplify(geoms, tolerance, progressFunction);
 }
 
 
@@ -59,10 +61,11 @@ CoverageSimplifier::simplify(
 std::vector<std::unique_ptr<Geometry>>
 CoverageSimplifier::simplifyInner(
     std::vector<const Geometry*>& coverage,
-    double tolerance)
+    double tolerance,
+    const ProgressFunction& progressFunction)
 {
     CoverageSimplifier simplifier(coverage);
-    return simplifier.simplifyInner(tolerance);
+    return simplifier.simplifyInner(tolerance, progressFunction);
 }
 
 
@@ -70,13 +73,14 @@ CoverageSimplifier::simplifyInner(
 std::vector<std::unique_ptr<Geometry>>
 CoverageSimplifier::simplifyInner(
     const std::vector<std::unique_ptr<Geometry>>& coverage,
-    double tolerance)
+    double tolerance,
+    const ProgressFunction& progressFunction)
 {
     std::vector<const Geometry*> geoms;
     for (auto& geom : coverage) {
         geoms.push_back(geom.get());
     }
-    return simplifyInner(geoms, tolerance);
+    return simplifyInner(geoms, tolerance, progressFunction);
 }
 
 
@@ -94,24 +98,28 @@ CoverageSimplifier::CoverageSimplifier(const std::vector<const Geometry*>& cover
 
 /* public */
 std::vector<std::unique_ptr<Geometry>>
-CoverageSimplifier::simplify(double tolerance)
+CoverageSimplifier::simplify(double tolerance,
+                             const ProgressFunction& progressFunction)
 {
-    CoverageRingEdges cov(m_input);
-    simplifyEdges(cov.getEdges(), nullptr, tolerance);
-    return cov.buildCoverage();
+    CoverageRingEdges cov(m_input, progressFunction.subProgress(0, 0.8));
+    simplifyEdges(cov.getEdges(), nullptr, tolerance, progressFunction.subProgress(0.8, 0.9));
+    return cov.buildCoverage(progressFunction.subProgress(0.9, 1.0));
 }
 
 /* public */
 std::vector<std::unique_ptr<Geometry>>
-CoverageSimplifier::simplifyInner(double tolerance)
+CoverageSimplifier::simplifyInner(double tolerance,
+                                  const ProgressFunction& progressFunction)
 {
-    CoverageRingEdges cov(m_input);
+    CoverageRingEdges cov(m_input, progressFunction.subProgress(0, 0.7));
     std::vector<CoverageEdge*> innerEdges = cov.selectEdges(2);
     std::vector<CoverageEdge*> outerEdges = cov.selectEdges(1);
-    std::unique_ptr<MultiLineString> constraintEdges = CoverageEdge::createLines(outerEdges, m_geomFactory);
 
-    simplifyEdges(innerEdges, constraintEdges.get(), tolerance);
-    return cov.buildCoverage();
+    std::unique_ptr<MultiLineString> constraintEdges = CoverageEdge::createLines(
+        outerEdges, m_geomFactory, progressFunction.subProgress(0.7, 0.8));
+
+    simplifyEdges(innerEdges, constraintEdges.get(), tolerance, progressFunction.subProgress(0.8, 0.9));
+    return cov.buildCoverage(progressFunction.subProgress(0.9, 1.0));
 }
 
 /* private */
@@ -119,11 +127,17 @@ void
 CoverageSimplifier::simplifyEdges(
     std::vector<CoverageEdge*> edges,
     const MultiLineString* constraints,
-    double tolerance)
+    double tolerance,
+    const ProgressFunction& progressFunction)
 {
-    std::unique_ptr<MultiLineString> lines = CoverageEdge::createLines(edges, m_geomFactory);
+    constexpr double RATIO_FIRST_PASS = 0.5;
+
+    std::unique_ptr<MultiLineString> lines = CoverageEdge::createLines(
+        edges, m_geomFactory, progressFunction.subProgress(0, RATIO_FIRST_PASS));
     std::vector<bool> freeRings = getFreeRings(edges);
-    std::unique_ptr<MultiLineString> linesSimp = TPVWSimplifier::simplify(lines.get(), freeRings, constraints, tolerance);
+    std::unique_ptr<MultiLineString> linesSimp = TPVWSimplifier::simplify(
+        lines.get(), freeRings, constraints, tolerance,
+        progressFunction.subProgress(RATIO_FIRST_PASS, 1.0));
     //Assert: mlsSimp.getNumGeometries = edges.length
 
     setCoordinates(edges, linesSimp.get());
diff --git a/src/coverage/TPVWSimplifier.cpp b/src/coverage/TPVWSimplifier.cpp
index 0bf36ac303..a6087c4227 100644
--- a/src/coverage/TPVWSimplifier.cpp
+++ b/src/coverage/TPVWSimplifier.cpp
@@ -47,10 +47,11 @@ typedef TPVWSimplifier::EdgeIndex EdgeIndex;
 std::unique_ptr<MultiLineString>
 TPVWSimplifier::simplify(
     const MultiLineString* lines,
-    double distanceTolerance)
+    double distanceTolerance,
+    const util::ProgressFunction& progressFunction)
 {
     TPVWSimplifier simp(lines, distanceTolerance);
-    std::unique_ptr<MultiLineString> result = simp.simplify();
+    std::unique_ptr<MultiLineString> result = simp.simplify(progressFunction);
     return result;
 }
 
@@ -61,12 +62,13 @@ TPVWSimplifier::simplify(
     const MultiLineString* p_lines,
     std::vector<bool>& p_freeRings,
     const MultiLineString* p_constraintLines,
-    double distanceTolerance)
+    double distanceTolerance,
+    const util::ProgressFunction& progressFunction)
 {
     TPVWSimplifier simp(p_lines, distanceTolerance);
     simp.setFreeRingIndices(p_freeRings);
     simp.setConstraints(p_constraintLines);
-    std::unique_ptr<MultiLineString> result = simp.simplify();
+    std::unique_ptr<MultiLineString> result = simp.simplify(progressFunction);
     return result;
 }
 
@@ -99,22 +101,42 @@ TPVWSimplifier::setFreeRingIndices(std::vector<bool>& freeRing)
 
 /* private */
 std::unique_ptr<MultiLineString>
-TPVWSimplifier::simplify()
+TPVWSimplifier::simplify(const util::ProgressFunction& progressFunction)
 {
     std::vector<bool> emptyList;
-    std::vector<Edge> edges = createEdges(inputLines, isFreeRing);
-    std::vector<Edge> constraintEdges = createEdges(constraintLines, emptyList);
+
+    constexpr double RATIO_FIRST_PASS = 0.8;
+
+    const double ratioInputLinesOverInputAndConstraint =
+        RATIO_FIRST_PASS *
+        static_cast<double>(inputLines ? inputLines->getNumGeometries() : 0) /
+        static_cast<double>(std::max<size_t>(1,
+            (inputLines ? inputLines->getNumGeometries() : 0) +
+            (constraintLines ? constraintLines->getNumGeometries() : 0)));
+
+    std::vector<Edge> edges = createEdges(inputLines, isFreeRing,
+        progressFunction.subProgress(0, ratioInputLinesOverInputAndConstraint));
+
+    std::vector<Edge> constraintEdges = createEdges(
+        constraintLines, emptyList,
+        progressFunction.subProgress(ratioInputLinesOverInputAndConstraint, RATIO_FIRST_PASS));
 
     EdgeIndex edgeIndex;
     edgeIndex.add(edges);
     edgeIndex.add(constraintEdges);
 
     std::vector<std::unique_ptr<LineString>> result;
+
+    util::ProgressContext progress(progressFunction.subProgress(RATIO_FIRST_PASS, 1.0), edges.size());
+
     for (auto& edge : edges) {
         std::unique_ptr<CoordinateSequence> ptsSimp = edge.simplify(edgeIndex);
         auto ls = geomFactory->createLineString(std::move(ptsSimp));
         result.emplace_back(ls.release());
+        progress.update();
     }
+
+    progress.finish();
     return geomFactory->createMultiLineString(std::move(result));
 }
 
@@ -122,18 +144,24 @@ TPVWSimplifier::simplify()
 std::vector<Edge>
 TPVWSimplifier::createEdges(
     const MultiLineString* lines,
-    std::vector<bool>& freeRing)
+    std::vector<bool>& freeRing,
+    const util::ProgressFunction& progressFunction)
 {
     std::vector<Edge> edges;
 
     if (lines == nullptr)
         return edges;
+    const size_t iterCount = lines->getNumGeometries();
+
+    util::ProgressContext progress(progressFunction, lines->getNumGeometries());
 
-    for (std::size_t i = 0; i < lines->getNumGeometries(); i++) {
+    for (size_t i = 0; i < iterCount; ++i) {
         const LineString* line = lines->getGeometryN(i);
         bool isFree = freeRing.empty() ? false : freeRing[i];
         edges.emplace_back(line, isFree, areaTolerance);
+        progress.update();
     }
+    progress.finish();
     return edges;
 }
 
diff --git a/src/coverage/VertexRingCounter.cpp b/src/coverage/VertexRingCounter.cpp
index fe742889f6..bccfde3383 100644
--- a/src/coverage/VertexRingCounter.cpp
+++ b/src/coverage/VertexRingCounter.cpp
@@ -34,12 +34,19 @@ namespace coverage { // geos.coverage
 void
 VertexRingCounter::count(
     const std::vector<const Geometry*>& geoms,
-    std::map<Coordinate, std::size_t>& counts)
+    std::map<Coordinate, std::size_t>& counts,
+    const util::ProgressFunction& progressFunction)
 {
     VertexRingCounter vertextCounter(counts);
+
+    util::ProgressContext progress(progressFunction, geoms.size());
+
     for (const Geometry* geom : geoms) {
         geom->apply_ro(vertextCounter);
+        progress.update();
     }
+
+    progress.finish();
 }
 
 
diff --git a/src/geom/Geometry.cpp b/src/geom/Geometry.cpp
index a34c7bf352..b833a04990 100644
--- a/src/geom/Geometry.cpp
+++ b/src/geom/Geometry.cpp
@@ -624,10 +624,10 @@ Geometry::Union(const Geometry* other) const
 
 /* public */
 Geometry::Ptr
-Geometry::Union() const
+Geometry::Union(geos::util::ProgressFunction* progressFunction) const
 {
     using geos::operation::geounion::UnaryUnionOp;
-    return UnaryUnionOp::Union(*this);
+    return UnaryUnionOp::Union(*this, progressFunction);
 }
 
 std::unique_ptr<Geometry>
diff --git a/src/geom/util/GeometryFixer.cpp b/src/geom/util/GeometryFixer.cpp
index 6d73ecbc9f..7523978f46 100644
--- a/src/geom/util/GeometryFixer.cpp
+++ b/src/geom/util/GeometryFixer.cpp
@@ -341,7 +341,8 @@ GeometryFixer::unionGeometry(std::vector<const Geometry*>& polys) const
     }
 
     UnaryUnionOp op(polys);
-    return op.Union();
+    geos::util::ProgressFunction* progressFunction = nullptr;
+    return op.Union(progressFunction);
     // return OverlayNGRobust::Union(polys);
 }
 
diff --git a/src/operation/overlayng/OverlayNGRobust.cpp b/src/operation/overlayng/OverlayNGRobust.cpp
index e9af20896a..fa12b9340b 100644
--- a/src/operation/overlayng/OverlayNGRobust.cpp
+++ b/src/operation/overlayng/OverlayNGRobust.cpp
@@ -76,7 +76,8 @@ OverlayNGRobust::Union(const Geometry* a)
     geounion::UnaryUnionOp op(*a);
     SRUnionStrategy unionSRFun;
     op.setUnionFunction(&unionSRFun);
-    return op.Union();
+    geos::util::ProgressFunction* progressFunction = nullptr;
+    return op.Union(progressFunction);
 }
 
 /*public static*/
diff --git a/src/operation/overlayng/UnaryUnionNG.cpp b/src/operation/overlayng/UnaryUnionNG.cpp
index f09b987b93..609fa93ee9 100644
--- a/src/operation/overlayng/UnaryUnionNG.cpp
+++ b/src/operation/overlayng/UnaryUnionNG.cpp
@@ -37,7 +37,8 @@ UnaryUnionNG::Union(const Geometry* geom, const PrecisionModel& pm)
     NGUnionStrategy ngUnionStrat(pm);
     geounion::UnaryUnionOp op(*geom);
     op.setUnionFunction(&ngUnionStrat);
-    return op.Union();
+    geos::util::ProgressFunction* progressFunction = nullptr;
+    return op.Union(progressFunction);
 }
 
 /*public static*/
diff --git a/src/operation/polygonize/BuildArea.cpp b/src/operation/polygonize/BuildArea.cpp
index 37ccd371eb..9f16ac6476 100644
--- a/src/operation/polygonize/BuildArea.cpp
+++ b/src/operation/polygonize/BuildArea.cpp
@@ -270,7 +270,7 @@ std::unique_ptr<geom::Geometry> BuildArea::build(const geom::Geometry* geom) {
 
     /* Run a single overlay operation to dissolve shared edges */
     auto shp = std::unique_ptr<geom::Geometry>(
-        geos::operation::geounion::CascadedPolygonUnion::CascadedPolygonUnion::Union(tmp.get()));
+        geos::operation::geounion::CascadedPolygonUnion::CascadedPolygonUnion::Union(tmp.get(), nullptr));
     if( shp ) {
         shp->setSRID(geom->getSRID());
     }
diff --git a/src/operation/union/CascadedPolygonUnion.cpp b/src/operation/union/CascadedPolygonUnion.cpp
index 82cd91bcc7..bc1b12d290 100644
--- a/src/operation/union/CascadedPolygonUnion.cpp
+++ b/src/operation/union/CascadedPolygonUnion.cpp
@@ -47,18 +47,19 @@ std::unique_ptr<geom::Geometry>
 CascadedPolygonUnion::Union(std::vector<geom::Polygon*>* polys)
 {
     CascadedPolygonUnion op(polys);
-    return op.Union();
+    geos::util::ProgressFunction* progressFunction = nullptr;
+    return op.Union(progressFunction);
 }
 
 std::unique_ptr<geom::Geometry>
-CascadedPolygonUnion::Union(std::vector<geom::Polygon*>* polys, UnionStrategy* unionFun)
+CascadedPolygonUnion::Union(std::vector<geom::Polygon*>* polys, UnionStrategy* unionFun, geos::util::ProgressFunction* progressFunction)
 {
     CascadedPolygonUnion op(polys, unionFun);
-    return op.Union();
+    return op.Union(progressFunction);
 }
 
 std::unique_ptr<geom::Geometry>
-CascadedPolygonUnion::Union(const geom::MultiPolygon* multipoly)
+CascadedPolygonUnion::Union(const geom::MultiPolygon* multipoly, geos::util::ProgressFunction* progressFunction)
 {
     std::vector<geom::Polygon*> polys;
 
@@ -67,11 +68,11 @@ CascadedPolygonUnion::Union(const geom::MultiPolygon* multipoly)
     }
 
     CascadedPolygonUnion op(&polys);
-    return op.Union();
+    return op.Union(progressFunction);
 }
 
 std::unique_ptr<geom::Geometry>
-CascadedPolygonUnion::Union()
+CascadedPolygonUnion::Union(geos::util::ProgressFunction* progressFunction)
 {
     if(inputPolys->empty()) {
         return nullptr;
@@ -94,28 +95,43 @@ CascadedPolygonUnion::Union()
     // TODO avoid creating this vector and run binaryUnion off the iterators directly
     std::vector<const geom::Geometry*> geoms(index.items().begin(), index.items().end());
 
-    return binaryUnion(geoms, 0, geoms.size());
+    size_t inc = 0;
+    std::function<void()> UnitProgress = [progressFunction, &inc, &geoms]()
+    {
+        ++inc;
+        (*progressFunction)(static_cast<double>(inc)/static_cast<double>(geoms.size()), "");
+    };
+
+    return binaryUnion(geoms, 0, geoms.size(), progressFunction ? &UnitProgress : nullptr);
 }
 
 
 std::unique_ptr<geom::Geometry>
 CascadedPolygonUnion::binaryUnion(const std::vector<const geom::Geometry*> & geoms,
-                                  std::size_t start, std::size_t end)
+                                  std::size_t start, std::size_t end,
+                                  std::function<void()>* unitProgress)
 {
     if(end - start == 0) {
         return nullptr;
     }
     else if(end - start == 1) {
+        if( unitProgress )
+            (*unitProgress)();
         return unionSafe(geoms[start], nullptr);
     }
     else if(end - start == 2) {
+        if( unitProgress )
+        {
+            (*unitProgress)();
+            (*unitProgress)();
+        }
         return unionSafe(geoms[start], geoms[start + 1]);
     }
     else {
         // recurse on both halves of the list
         std::size_t mid = (end + start) / 2;
-        std::unique_ptr<geom::Geometry> g0(binaryUnion(geoms, start, mid));
-        std::unique_ptr<geom::Geometry> g1(binaryUnion(geoms, mid, end));
+        std::unique_ptr<geom::Geometry> g0(binaryUnion(geoms, start, mid, unitProgress));
+        std::unique_ptr<geom::Geometry> g1(binaryUnion(geoms, mid, end, unitProgress));
         return unionSafe(std::move(g0), std::move(g1));
     }
 }
diff --git a/src/operation/union/UnaryUnionOp.cpp b/src/operation/union/UnaryUnionOp.cpp
index a801e56582..c894d7cfd3 100644
--- a/src/operation/union/UnaryUnionOp.cpp
+++ b/src/operation/union/UnaryUnionOp.cpp
@@ -64,7 +64,7 @@ UnaryUnionOp::unionWithNull(std::unique_ptr<geom::Geometry> g0,
 
 /*public*/
 std::unique_ptr<geom::Geometry>
-UnaryUnionOp::Union()
+UnaryUnionOp::Union(geos::util::ProgressFunction* progressFunction)
 {
     typedef std::unique_ptr<geom::Geometry> GeomPtr;
 
@@ -95,7 +95,7 @@ UnaryUnionOp::Union()
 
     GeomPtr unionPolygons;
     if(!polygons.empty()) {
-        unionPolygons = CascadedPolygonUnion::Union(polygons.begin(), polygons.end(), unionFunction);
+        unionPolygons = CascadedPolygonUnion::Union(polygons.begin(), polygons.end(), unionFunction, progressFunction);
     }
 
     /*
diff --git a/src/util/Progress.cpp b/src/util/Progress.cpp
new file mode 100644
index 0000000000..4a89c62d4f
--- /dev/null
+++ b/src/util/Progress.cpp
@@ -0,0 +1,31 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2025 Even Rouault <even.rouault@spatialys.com>
+ *
+ * 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 <geos/util/Progress.h>
+
+namespace geos::util {
+
+ProgressFunction
+ProgressFunction::subProgress(double from, double to) const {
+    if (m_function.has_value()) {
+        return ProgressFunction([from, to, this](double ratio, const char* msg)
+        {
+            m_function.value()(from + (to - from) * ratio, msg);
+        });
+    }
+
+    return ProgressFunction();
+}
+
+}
\ No newline at end of file
diff --git a/tests/unit/capi/GEOSCoverageSimplifyTest.cpp b/tests/unit/capi/GEOSCoverageSimplifyTest.cpp
index 4ca7e9c394..b3c5703161 100644
--- a/tests/unit/capi/GEOSCoverageSimplifyTest.cpp
+++ b/tests/unit/capi/GEOSCoverageSimplifyTest.cpp
@@ -84,8 +84,21 @@ template<> void object::test<3>
     const char* inputWKT = "GEOMETRYCOLLECTION(POLYGON(( 0 0,10 0,10.1 5,10 10,0 10,0 0)),POLYGON((10 0,20 0,20 10,10 10,10.1 5,10 0)))";
 
     input_ = fromWKT(inputWKT);
-    result_ = GEOSCoverageSimplifyVW(input_, 1.0, 0);
 
+    struct Cbk
+    {
+        double lastRatio = 0;
+
+        static void Func(double progressRatio, const char* /* msg */, void* userdata)
+        {
+            Cbk* self = static_cast<Cbk*>(userdata);
+            self->lastRatio = progressRatio;
+        }
+    };
+
+    Cbk cbk;
+    result_ = GEOSCoverageSimplifyVWWithProgress(input_, 1.0, 0, &Cbk::Func, &cbk);
+    ensure("cbk.lastRatio == 1.0", cbk.lastRatio == 1.0);
     ensure( result_ != nullptr );
     ensure( GEOSGeomTypeId(result_) == GEOS_GEOMETRYCOLLECTION );
 
diff --git a/tests/unit/coverage/CoverageSimplifierTest.cpp b/tests/unit/coverage/CoverageSimplifierTest.cpp
index 1973ae965e..725ddcdd89 100644
--- a/tests/unit/coverage/CoverageSimplifierTest.cpp
+++ b/tests/unit/coverage/CoverageSimplifierTest.cpp
@@ -43,8 +43,15 @@ struct test_coveragesimplifier_data {
         double tolerance,
         const std::vector<std::unique_ptr<Geometry>>& expected)
     {
-        std::vector<std::unique_ptr<Geometry>> actual = CoverageSimplifier::simplify(input, tolerance);
+        double lastRatio = 0.0;
+        geos::util::ProgressFunction myProgress([&lastRatio](double ratio, const char*)
+        {
+            ensure("ratio >= lastRatio", ratio >= lastRatio);
+            lastRatio = ratio;
+        });
+        std::vector<std::unique_ptr<Geometry>> actual = CoverageSimplifier::simplify(input, tolerance, myProgress);
         checkArrayEqual(expected, actual);
+        ensure("lastRatio == 1.0", lastRatio == 1.0);
     }
 
     void checkResultInner(