[postgis-tickets] [SCM] PostGIS branch master updated. 3.3.0rc2-746-g495b926e3
git at osgeo.org
git at osgeo.org
Wed May 3 13:55:08 PDT 2023
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 "PostGIS".
The branch, master has been updated
via 495b926e31f0bbef0980da4130f2e14fbdf911ab (commit)
from bef103b65ccb5f5ed6f4a8bb346a6e4c95c72d0e (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 495b926e31f0bbef0980da4130f2e14fbdf911ab
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Wed May 3 13:54:49 2023 -0700
Add three new functions,
ST_CoverageIsValid(geom geometry winset, tolerance float8)
ST_CoverageSimplify(geom geometry winset, tolerance float8, simplifyBoundary integer)
ST_CoverageUnion(geom geometry set)
For the purpose of these functions a "coverage" is a set of polygons with perfectly shared edges (no overlaps or gaps and exact vertex matching). For associated GEOS discussion, see libgeos/geos#867
diff --git a/doc/Makefile.in b/doc/Makefile.in
index 8175fb5c8..18f7130fe 100644
--- a/doc/Makefile.in
+++ b/doc/Makefile.in
@@ -144,6 +144,7 @@ XML_SOURCES = \
reference_output.xml \
reference_overlay.xml \
reference_processing.xml \
+ reference_coverage.xml \
reference_raster.xml \
reference_relationship.xml \
reference_srs.xml \
diff --git a/doc/postgis.xml b/doc/postgis.xml
index 4714a0fb5..6d5b713b8 100644
--- a/doc/postgis.xml
+++ b/doc/postgis.xml
@@ -46,6 +46,7 @@
<!ENTITY reference_sfcgal SYSTEM "reference_sfcgal.xml">
<!ENTITY reference_overlay SYSTEM "reference_overlay.xml">
<!ENTITY reference_processing SYSTEM "reference_processing.xml">
+<!ENTITY reference_coverage SYSTEM "reference_coverage.xml">
<!ENTITY reference_lrs SYSTEM "reference_lrs.xml">
<!ENTITY reference_srs SYSTEM "reference_srs.xml">
<!ENTITY reference_trajectory SYSTEM "reference_trajectory.xml">
diff --git a/doc/reference.xml b/doc/reference.xml
index a74281054..232fa965b 100644
--- a/doc/reference.xml
+++ b/doc/reference.xml
@@ -30,6 +30,7 @@
&reference_measure;
&reference_overlay;
&reference_processing;
+ &reference_coverage;
&reference_transformation;
&reference_cluster;
&reference_bbox;
diff --git a/doc/reference_coverage.xml b/doc/reference_coverage.xml
new file mode 100644
index 000000000..99bd0a8aa
--- /dev/null
+++ b/doc/reference_coverage.xml
@@ -0,0 +1,195 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<sect1 id="Coverage">
+ <sect1info>
+ <abstract>
+ <para>These functions operate on collections of polygonal geometry that form "implicit coverages" in which any edges shared by the polygons match exactly. These implicit coverages are fast to process, and can be operated on with window functions, which retain the implicit topology inside the window partition while altering the edges.</para>
+ </abstract>
+ </sect1info>
+
+ <title>Coverages</title>
+
+ <refentry id="ST_CoverageInvalidLocations">
+ <refnamediv>
+ <refname>ST_CoverageInvalidLocations</refname>
+
+ <refpurpose>Window function that evaluates a set of polygons in the window partition to determine if they form a valid polygonal coverage.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_CoverageInvalidLocations</function></funcdef>
+ <paramdef><type>geometry winset </type>
+ <parameter>geom</parameter></paramdef>
+ <paramdef><type>float8 </type>
+ <parameter>tolerance</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>Collects all the polygons in the window partition and then checks whether they form a "clean coverage", in which no polygons have any overlap, all vertices on shared edges are identical, and no gaps between polygons exist smaller than the tolerance value.</para>
+ <para>where polygons do not form a valid coverage, the problematic edges will be returned as MULTILINESTRING. For clean polygons the return value will be a NULL. Any non-polygonal or empty geometries on the input will be generate NULL returns.</para>
+ <para>Availability: 3.4.0 - requires GEOS >= 3.12.0</para>
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+
+ <programlisting>CREATE TABLE coverage (id integer, geom geometry);
+
+INSERT INTO coverage VALUES
+ (1, 'POLYGON((0 0, 10 0, 10.1 5, 10 10, 0 10, 0 0))'),
+ (2, 'POLYGON((10 0, 20 0, 20 10, 10 10, 10.05 5, 10 0))'),
+ (3, 'POLYGON((20 0, 30 0, 30 10, 20 10, 20 0))');
+
+SELECT id, ST_AsText(ST_CoverageInvalidLocations(geom) OVER ())
+ FROM coverage;
+
+ id | st_astext
+----+--------------------------------
+ 1 | LINESTRING(10 0,10.1 5,10 10)
+ 2 | LINESTRING(10 10,10.05 5,10 0)
+ 3 |
+ </programlisting>
+
+ <programlisting>-- Test entire table for validity
+SELECT true = ALL (
+ SELECT ST_CoverageInvalidLocations(geom) OVER () IS NULL
+ FROM coverage
+ );
+ </programlisting>
+ </refsection>
+
+ <refsection>
+ <title>See Also</title>
+ <para>
+ <xref linkend="ST_CoverageUnion" />,
+ <xref linkend="ST_CoverageSimplify" />
+ </para>
+ </refsection>
+
+ </refentry>
+
+
+ <refentry id="ST_CoverageSimplify">
+ <refnamediv>
+ <refname>ST_CoverageSimplify</refname>
+
+ <refpurpose>Window function that simplifies the edges of an implicit coverage formed by a set of polygons, applying the Visvalingam–Whyatt algorithm to the edges.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_CoverageSimplify</function></funcdef>
+ <paramdef><type>geometry winset </type>
+ <parameter>geom</parameter></paramdef>
+ <paramdef><type>float8 </type>
+ <parameter>tolerance</parameter></paramdef>
+ <paramdef choice="opt"><type>boolean </type>
+ <parameter>simplifyBoundary=true</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>Simplifies shared edges of a polygonal coverage using the <ulink url="https://en.wikipedia.org/wiki/Visvalingam%E2%80%93Whyatt_algorithm">Visvalingam–Whyatt algorithm</ulink>. For clean polygonal coverage inputs, this will result in a simplification that is consistent between polygons, and retains the coverage property for the outputs.</para>
+ <para>For invalid coverage inputs, non-shared edges are ignored, and only shared edges are simplified.</para>
+ <para>To only simplify the "internal" edges of the set of polygons (those that are shared by two polygons) set the "simplifyBoundary" parameter to false.</para>
+ <note><para>You must ensure that your coverage is valid before running simplification or you may get unexpected results in the output, such as boundary intersections or non-shared boundaries where you thought there were shared boundaries.</para></note>
+ <para>Availability: 3.4.0 - requires GEOS >= 3.12.0</para>
+
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+
+ <programlisting>CREATE TABLE coverage (id integer, geom geometry);
+
+INSERT INTO coverage VALUES
+ (1, 'POLYGON((0 0, 10 0, 10.1 5, 10 10, 0 10, 0 0))'),
+ (2, 'POLYGON((10 0, 20 0, 20 10, 10 10, 10.1 5, 10 0))'),
+ (3, 'POLYGON((20 0, 30 0, 30 10, 20 10, 20 0))');
+
+SELECT id, ST_AsText(ST_CoverageSimplify(geom, 1.0) OVER ())
+ FROM coverage;
+
+ id | st_astext
+----+---------------------------------------
+ 1 | POLYGON((10 0,10 10,0 10,0 0,10 0))
+ 2 | POLYGON((10 0,20 0,20 10,10 10,10 0))
+ 3 | POLYGON((20 0,30 0,30 10,20 10,20 0))
+ </programlisting>
+ </refsection>
+
+ <refsection>
+ <title>See Also</title>
+ <para>
+ <xref linkend="ST_CoverageUnion" />,
+ <xref linkend="ST_CoverageInvalidLocations" />
+ </para>
+ </refsection>
+
+ </refentry>
+
+ <refentry id="ST_CoverageUnion">
+ <refnamediv>
+ <refname>ST_CoverageUnion</refname>
+
+ <refpurpose>Removes adjacent edges in a polygonal coverage, thereby unioning neighboring polygons.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_CoverageUnion</function></funcdef>
+ <paramdef><type>geometry set</type>
+ <parameter>geom</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>When the inputs are a polygonal coverage, the function will provide a high speed union of adjacent polygons in the aggregation set. This function only works with inputs that are coverages: no overlaps, all shared edges identical down to the vertex level. See <xref linkend="ST_CoverageInvalidLocations" /> for a way to test that a set of polygons are a valid polygonal coverage.</para>
+
+ <note><para>You must ensure that your coverage is valid before running simplification or you may get unexpected results in the output, such as unmerged inputs returning as multipolygons.</para></note>
+
+ <para>Availability: 3.4.0 - requires GEOS >= 3.8.0</para>
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+
+ <programlisting>CREATE TABLE coverage (id integer, geom geometry);
+
+INSERT INTO coverage VALUES
+ (1, 'POLYGON((0 0, 10 0, 10.1 5, 10 10, 0 10, 0 0))'),
+ (2, 'POLYGON((10 0, 20 0, 20 10, 10 10, 10.1 5, 10 0))'),
+ (3, 'POLYGON((20 0, 30 0, 30 10, 20 10, 20 0))');
+
+SELECT ST_AsText(ST_CoverageUnion(geom))
+ FROM coverage;
+
+ st_astext
+----------------------------------------------------------
+ POLYGON((0 0,0 10,10 10,20 10,30 10,30 0,20 0,10 0,0 0))
+ </programlisting>
+ </refsection>
+
+ <refsection>
+ <title>See Also</title>
+ <para>
+ <xref linkend="ST_CoverageSimplify" />,
+ <xref linkend="ST_CoverageInvalidLocations" />
+ </para>
+ </refsection>
+
+ </refentry>
+</sect1>
diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c
index ce2505954..56b38b0be 100644
--- a/liblwgeom/lwgeom_geos.c
+++ b/liblwgeom/lwgeom_geos.c
@@ -307,16 +307,14 @@ ptarray_to_GEOSCoordSeq(const POINTARRAY* pa, uint8_t fix_ring)
sq = GEOSCoordSeq_copyFromBuffer((const double*) pa->serialized_pointlist, pa->npoints, FLAGS_GET_Z(pa->flags), FLAGS_GET_M(pa->flags));
if (!sq)
{
- lwerror("Error creating GEOS Coordinate Sequence");
- return NULL;
+ GEOS_FAIL();
}
return sq;
}
#endif
- if (!(sq = GEOSCoordSeq_create(pa->npoints + append_points, dims)))
+ if (!(sq = GEOSCoordSeq_create(pa->npoints + append_points, dims)))
{
- lwerror("Error creating GEOS Coordinate Sequence");
- return NULL;
+ GEOS_FAIL();
}
for (i = 0; i < pa->npoints; i++)
diff --git a/liblwgeom/lwgeom_topo.c b/liblwgeom/lwgeom_topo.c
index 7ecc5e430..d85ed6847 100644
--- a/liblwgeom/lwgeom_topo.c
+++ b/liblwgeom/lwgeom_topo.c
@@ -5722,7 +5722,7 @@ _lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges,
continue;
}
nearby[nearbyindex++] = g;
- ++nn;
+ nn = nn + 1;
}
LWDEBUGF(1, "Found %d isolated nodes closer than tolerance (%g)", nn, tol);
}}
diff --git a/postgis/lwgeom_accum.c b/postgis/lwgeom_accum.c
index 4ec33c852..973c40af4 100644
--- a/postgis/lwgeom_accum.c
+++ b/postgis/lwgeom_accum.c
@@ -56,6 +56,7 @@ Datum polygonize_garray(PG_FUNCTION_ARGS);
Datum clusterintersecting_garray(PG_FUNCTION_ARGS);
Datum cluster_within_distance_garray(PG_FUNCTION_ARGS);
Datum LWGEOM_makeline_garray(PG_FUNCTION_ARGS);
+Datum ST_CoverageUnion(PG_FUNCTION_ARGS);
/**
@@ -330,6 +331,34 @@ pgis_geometry_clusterwithin_finalfn(PG_FUNCTION_ARGS)
PG_RETURN_DATUM(result);
}
+
+
+/**
+* The "coverage union" final function passes the geometry[] to a
+* GEOSCoverageUnion call before returning the result.
+*/
+PG_FUNCTION_INFO_V1(pgis_geometry_coverageunion_finalfn);
+Datum
+pgis_geometry_coverageunion_finalfn(PG_FUNCTION_ARGS)
+{
+ CollectionBuildState *p;
+ Datum result = 0;
+ Datum geometry_array = 0;
+
+ if (PG_ARGISNULL(0))
+ PG_RETURN_NULL(); /* returns null iff no input values */
+
+ p = (CollectionBuildState*) PG_GETARG_POINTER(0);
+
+ geometry_array = pgis_accum_finalfn(p, CurrentMemoryContext, fcinfo);
+ result = PGISDirectFunctionCall1(ST_CoverageUnion, geometry_array);
+ if (!result)
+ PG_RETURN_NULL();
+
+ PG_RETURN_DATUM(result);
+}
+
+
/**
* A modified version of PostgreSQL's DirectFunctionCall1 which allows NULL results; this
* is required for aggregates that return NULL.
diff --git a/postgis/lwgeom_geos.c b/postgis/lwgeom_geos.c
index f4c2d919d..1e166625f 100644
--- a/postgis/lwgeom_geos.c
+++ b/postgis/lwgeom_geos.c
@@ -46,20 +46,7 @@
#include "lwgeom_accum.h"
-/* Return NULL on GEOS error
- *
- * Prints error message only if it was not for interruption, in which
- * case we let PostgreSQL deal with the error.
- */
-#define HANDLE_GEOS_ERROR(label) \
- { \
- if (strstr(lwgeom_geos_errmsg, "InterruptedException")) \
- ereport(ERROR, \
- (errcode(ERRCODE_QUERY_CANCELED), errmsg("canceling statement due to user request"))); \
- else \
- lwpgerror("%s: %s", (label), lwgeom_geos_errmsg); \
- PG_RETURN_NULL(); \
- }
+
/*
** Prototypes for SQL-bound functions
@@ -3589,3 +3576,6 @@ Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS)
PG_FREE_IF_COPY(input, 0);
PG_RETURN_POINTER(result);
}
+
+
+
diff --git a/postgis/lwgeom_geos.h b/postgis/lwgeom_geos.h
index 1a834c416..2fb5cf93c 100644
--- a/postgis/lwgeom_geos.h
+++ b/postgis/lwgeom_geos.h
@@ -49,4 +49,20 @@ Datum ST_3DDistance(PG_FUNCTION_ARGS);
uint32_t array_nelems_not_null(ArrayType* array);
+
+/* Return NULL on GEOS error
+ *
+ * Prints error message only if it was not for interruption, in which
+ * case we let PostgreSQL deal with the error.
+ */
+#define HANDLE_GEOS_ERROR(label) \
+ { \
+ if (strstr(lwgeom_geos_errmsg, "InterruptedException")) \
+ ereport(ERROR, \
+ (errcode(ERRCODE_QUERY_CANCELED), errmsg("canceling statement due to user request"))); \
+ else \
+ lwpgerror("%s: %s", (label), lwgeom_geos_errmsg); \
+ PG_RETURN_NULL(); \
+ }
+
#endif /* LWGEOM_GEOS_H_ */
diff --git a/postgis/lwgeom_window.c b/postgis/lwgeom_window.c
index 4fd14199d..af34e70fb 100644
--- a/postgis/lwgeom_window.c
+++ b/postgis/lwgeom_window.c
@@ -37,10 +37,6 @@
#include "lwgeom_log.h"
#include "lwgeom_pg.h"
-extern Datum ST_ClusterDBSCAN(PG_FUNCTION_ARGS);
-extern Datum ST_ClusterKMeans(PG_FUNCTION_ARGS);
-extern Datum ST_ClusterWithinWin(PG_FUNCTION_ARGS);
-extern Datum ST_ClusterIntersectingWin(PG_FUNCTION_ARGS);
typedef struct {
bool isdone;
@@ -117,6 +113,8 @@ read_geos_from_partition(WindowObject win_obj, uint32_t i, bool* is_null)
}
return gg;
}
+
+extern Datum ST_ClusterDBSCAN(PG_FUNCTION_ARGS);
PG_FUNCTION_INFO_V1(ST_ClusterDBSCAN);
Datum ST_ClusterDBSCAN(PG_FUNCTION_ARGS)
{
@@ -209,7 +207,7 @@ Datum ST_ClusterDBSCAN(PG_FUNCTION_ARGS)
PG_RETURN_INT32(context->clusters[row].cluster_id);
}
-
+extern Datum ST_ClusterWithinWin(PG_FUNCTION_ARGS);
PG_FUNCTION_INFO_V1(ST_ClusterWithinWin);
Datum ST_ClusterWithinWin(PG_FUNCTION_ARGS)
{
@@ -285,6 +283,7 @@ Datum ST_ClusterWithinWin(PG_FUNCTION_ARGS)
PG_RETURN_INT32(context->clusters[row].cluster_id);
}
+extern Datum ST_ClusterIntersectingWin(PG_FUNCTION_ARGS);
PG_FUNCTION_INFO_V1(ST_ClusterIntersectingWin);
Datum ST_ClusterIntersectingWin(PG_FUNCTION_ARGS)
{
@@ -350,6 +349,7 @@ Datum ST_ClusterIntersectingWin(PG_FUNCTION_ARGS)
}
+extern Datum ST_ClusterKMeans(PG_FUNCTION_ARGS);
PG_FUNCTION_INFO_V1(ST_ClusterKMeans);
Datum ST_ClusterKMeans(PG_FUNCTION_ARGS)
{
@@ -451,3 +451,423 @@ Datum ST_ClusterKMeans(PG_FUNCTION_ARGS)
curpos = WinGetCurrentPosition(winobj);
PG_RETURN_INT32(context->result[curpos]);
}
+
+
+/********************************************************************************
+ * GEOS Coverage Functions
+ *
+ * The GEOS "coverage" support is a set of functions for working with
+ * "implicit coverages". That is, collections of polygonal geometry
+ * that have perfectly shared edges. Since perfectly matched edges are
+ * fast to detect, "building a coverage" on these inputs is fast, so
+ * operations like edge simplification and so on can include the "build"
+ * step with relatively low cost.
+ *
+ */
+
+#if POSTGIS_GEOS_VERSION >= 30800
+/*
+ * For CoverageUnion and CoverageSimplify, clean up
+ * the start of a collection if things fail in mid-stream.
+ */
+static void
+coverage_destroy_geoms(GEOSGeometry **geoms, uint32 ngeoms)
+{
+ if (!geoms) return;
+ if (!ngeoms) return;
+ for (uint32 i = 0; i < ngeoms; i++)
+ {
+ if(geoms[i])
+ GEOSGeom_destroy(geoms[i]);
+ }
+}
+#endif
+
+#if POSTGIS_GEOS_VERSION >= 31200
+
+/*
+ * For CoverageIsValid and CoverageSimplify, maintain context
+ * of unioned output and the position of the resultant in that
+ * output relative to the index of the input geometries.
+ */
+typedef struct {
+ bool isdone;
+ bool isnull;
+ LWCOLLECTION *geom;
+ int64 idx[0];
+ /* variable length */
+} coverage_context;
+
+/*
+ * For CoverageIsValid and CoverageSimplify, read the context
+ * out of the WindowObject.
+ */
+static coverage_context *
+coverage_context_fetch(WindowObject winobj, int64 rowcount)
+{
+ size_t ctx_size = sizeof(coverage_context) + rowcount * sizeof(int64);
+ coverage_context *ctx = (coverage_context*) WinGetPartitionLocalMemory(
+ winobj, ctx_size);
+ return ctx;
+}
+
+/*
+ * For CoverageIsValid and CoverageSimplify, read all the
+ * geometries in this partition and form a GEOS geometry
+ * collection out of them.
+ */
+static GEOSGeometry*
+coverage_read_partition_into_collection(
+ WindowObject winobj,
+ coverage_context* context)
+{
+ int64 rowcount = WinGetPartitionRowCount(winobj);
+ GEOSGeometry* geos;
+ GEOSGeometry** geoms;
+ uint32 i, ngeoms = 0, gtype;
+
+ /* Read in all the geometries in this partition */
+ geoms = palloc(rowcount * sizeof(GEOSGeometry*));
+ for (i = 0; i < rowcount; i++)
+ {
+ GSERIALIZED* gser;
+ bool isnull, isout;
+ bool isempty, ispolygonal;
+ Datum d;
+
+ /* Read geometry in first argument */
+ d = WinGetFuncArgInPartition(winobj, 0, i,
+ WINDOW_SEEK_HEAD, false, &isnull, &isout);
+
+ /*
+ * The input to the GEOS function will be smaller than
+ * the input to this window, since we will not feed
+ * GEOS nulls or empties. So we need to maintain a
+ * map (context->idx) from the window position of the
+ * input to the GEOS position, so we can put the
+ * right result in the output stream.
+ */
+
+ /* Skip NULL inputs and move on */
+ if (isnull)
+ {
+ context->idx[i] = -1;
+ continue;
+ }
+
+ gser = (GSERIALIZED*)PG_DETOAST_DATUM(d);
+ gtype = gserialized_get_type(gser);
+ isempty = gserialized_is_empty(gser);
+ ispolygonal = (gtype == MULTIPOLYGONTYPE) || (gtype == POLYGONTYPE);
+
+ /* Skip empty or non-polygonal inputs */
+ if (isempty || !ispolygonal)
+ {
+ context->idx[i] = -1;
+ continue;
+ }
+
+ /* Skip failed inputs */
+ geos = POSTGIS2GEOS(gser);
+ if (!geos)
+ {
+ context->idx[i] = -1;
+ continue;
+ }
+
+ context->idx[i] = ngeoms;
+ geoms[ngeoms] = geos;
+ ngeoms = ngeoms + 1;
+ }
+
+ /*
+ * Create the GEOS input collection! The new
+ * collection takes ownership of the input GEOSGeometry
+ * objects, leaving just the ngeoms array, which
+ * will be cleaned up on function exit.
+ */
+ geos = GEOSGeom_createCollection(
+ GEOS_GEOMETRYCOLLECTION,
+ geoms, ngeoms);
+
+ /*
+ * If the creation failed, the objects will still be
+ * hanging around, so clean them up first. Should
+ * never happen, really.
+ */
+ if (!geos)
+ {
+ coverage_destroy_geoms(geoms, ngeoms);
+ return NULL;
+ }
+
+ pfree(geoms);
+ return geos;
+}
+
+
+/*
+ * The coverage_window_calculation function is just
+ * the shared machinery of the two window functions
+ * CoverageIsValid and CoverageSimplify which are almost
+ * the same except the GEOS function they call. That
+ * operating mode is controlled with this enumeration.
+ */
+enum {
+ COVERAGE_SIMPLIFY = 0,
+ COVERAGE_ISVALID = 1
+};
+
+/*
+ * This calculation is shared by both coverage operations
+ * since they have the same pattern of "consume collection,
+ * return collection", and are both window functions.
+ */
+static Datum
+coverage_window_calculation(PG_FUNCTION_ARGS, int mode)
+{
+ WindowObject winobj = PG_WINDOW_OBJECT();
+ int64 curpos = WinGetCurrentPosition(winobj);
+ int64 rowcount = WinGetPartitionRowCount(winobj);
+ coverage_context *context = coverage_context_fetch(winobj, rowcount);
+ GSERIALIZED* result;
+ MemoryContext uppercontext = fcinfo->flinfo->fn_mcxt;
+ MemoryContext oldcontext;
+ const LWGEOM* subgeom;
+
+ if (!context->isdone)
+ {
+ bool isnull;
+ Datum d;
+ double tolerance = 0.0;
+ bool simplifyBoundary = true;
+ GEOSGeometry *output = NULL;
+ GEOSGeometry *input = NULL;
+
+ if (!fcinfo->flinfo)
+ elog(ERROR, "%s: Could not find upper context", __func__);
+
+ if (rowcount == 0)
+ {
+ context->isdone = true;
+ context->isnull = true;
+ PG_RETURN_NULL();
+ }
+
+ /* Get the tolerance argument from second postition */
+ d = WinGetFuncArgCurrent(winobj, 1, &isnull);
+ if (!isnull)
+ tolerance = DatumGetFloat8(d);
+
+ /* The third position "preserve boundary" argument */
+ /* is only for the simplify mode */
+ if (mode == COVERAGE_SIMPLIFY)
+ {
+ d = WinGetFuncArgCurrent(winobj, 2, &isnull);
+ if (!isnull)
+ simplifyBoundary = DatumGetBool(d);
+ }
+
+ initGEOS(lwnotice, lwgeom_geos_error);
+
+ input = coverage_read_partition_into_collection(winobj, context);
+ if (!input)
+ HANDLE_GEOS_ERROR("Failed to create collection");
+
+ /* Run the correct GEOS function for the calling mode */
+ if (mode == COVERAGE_SIMPLIFY)
+ {
+ /* GEOSCoverageSimplifyVW is "preserveBoundary" so we invert simplifyBoundary */
+ output = GEOSCoverageSimplifyVW(input, tolerance, !simplifyBoundary);
+ }
+ else if (mode == COVERAGE_ISVALID)
+ {
+ GEOSCoverageIsValid(input, tolerance, &output);
+ }
+ else
+ {
+ elog(ERROR, "Unknown mode, never get here");
+ }
+
+ GEOSGeom_destroy(input);
+
+ if (!output)
+ {
+ HANDLE_GEOS_ERROR("Failed to process collection");
+ }
+
+ oldcontext = MemoryContextSwitchTo(uppercontext);
+ context->geom = (LWCOLLECTION *) GEOS2LWGEOM(output, GEOSHasZ(output));
+ MemoryContextSwitchTo(oldcontext);
+ GEOSGeom_destroy(output);
+
+ context->isdone = true;
+ }
+
+ /* Bomb out of any errors */
+ if (context->isnull)
+ PG_RETURN_NULL();
+
+ /* Propogate the null entries */
+ if (context->idx[curpos] < 0)
+ PG_RETURN_NULL();
+
+ /*
+ * Geometry serialization is not const-safe! (we
+ * calculate bounding boxes on demand) so we need
+ * to make sure we have a persistent context when
+ * we call the serialization, lest we create dangling
+ * pointers in the object. It's possible we could
+ * ignore them, skip the manual lwcollection_free
+ * and let the aggcontext deletion take
+ * care of the memory.
+ */
+ oldcontext = MemoryContextSwitchTo(uppercontext);
+ subgeom = lwcollection_getsubgeom(context->geom, context->idx[curpos]);
+ if (mode == COVERAGE_ISVALID && lwgeom_is_empty(subgeom))
+ {
+ result = NULL;
+ }
+ else
+ {
+ result = geometry_serialize((LWGEOM*)subgeom);
+ }
+ MemoryContextSwitchTo(oldcontext);
+
+ /* When at the end of the partition, release the */
+ /* simplified collection we have been reading. */
+ if (curpos == rowcount - 1)
+ lwcollection_free(context->geom);
+
+ if (!result) PG_RETURN_NULL();
+
+ PG_RETURN_POINTER(result);
+
+}
+#endif /* POSTGIS_GEOS_VERSION >= 31200 */
+
+
+extern Datum ST_CoverageSimplify(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(ST_CoverageSimplify);
+Datum ST_CoverageSimplify(PG_FUNCTION_ARGS)
+{
+#if POSTGIS_GEOS_VERSION < 31200
+
+ lwpgerror("The GEOS version this PostGIS binary "
+ "was compiled against (%d) doesn't support "
+ "'GEOSCoverageSimplifyVW' function (3.12.0+ required)",
+ POSTGIS_GEOS_VERSION);
+ PG_RETURN_NULL();
+
+#else /* POSTGIS_GEOS_VERSION >= 31200 */
+
+ return coverage_window_calculation(fcinfo, COVERAGE_SIMPLIFY);
+
+#endif
+}
+
+
+extern Datum ST_CoverageInvalidLocations(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(ST_CoverageInvalidLocations);
+Datum ST_CoverageInvalidLocations(PG_FUNCTION_ARGS)
+{
+#if POSTGIS_GEOS_VERSION < 31200
+
+ lwpgerror("The GEOS version this PostGIS binary "
+ "was compiled against (%d) doesn't support "
+ "'GEOSCoverageIsValid' function (3.12.0+ required)",
+ POSTGIS_GEOS_VERSION);
+ PG_RETURN_NULL();
+
+#else /* POSTGIS_GEOS_VERSION >= 31200 */
+
+ return coverage_window_calculation(fcinfo, COVERAGE_ISVALID);
+
+#endif
+}
+
+/**********************************************************************
+ * ST_CoverageUnion(geometry[])
+ *
+ */
+
+Datum ST_CoverageUnion(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(ST_CoverageUnion);
+Datum ST_CoverageUnion(PG_FUNCTION_ARGS)
+{
+#if POSTGIS_GEOS_VERSION < 30800
+ lwpgerror("The GEOS version this PostGIS binary "
+ "was compiled against (%d) doesn't support "
+ "'GEOSCoverageUnion' function (3.8.0+ required)",
+ POSTGIS_GEOS_VERSION);
+ PG_RETURN_NULL();
+
+#else /* POSTGIS_GEOS_VERSION >= 30800 */
+
+ GSERIALIZED *result = NULL;
+
+ Datum value;
+ bool isnull;
+
+ GEOSGeometry **geoms = NULL;
+ GEOSGeometry *geos = NULL;
+ GEOSGeometry *geos_result = NULL;
+ uint32 ngeoms = 0;
+
+ ArrayType *array = PG_GETARG_ARRAYTYPE_P(0);
+ uint32 nelems = ArrayGetNItems(ARR_NDIM(array), ARR_DIMS(array));
+ ArrayIterator iterator = array_create_iterator(array, 0, NULL);
+
+ /* Return null on 0-elements input array */
+ if (nelems == 0)
+ PG_RETURN_NULL();
+
+ /* Convert all geometries into GEOSGeometry array */
+ geoms = palloc(sizeof(GEOSGeometry *) * nelems);
+
+ initGEOS(lwpgnotice, lwgeom_geos_error);
+
+ while (array_iterate(iterator, &value, &isnull))
+ {
+ GSERIALIZED *gser;
+ /* Omit nulls */
+ if (isnull) continue;
+
+ /* Omit empty */
+ gser = (GSERIALIZED *)DatumGetPointer(value);
+ if (gserialized_is_empty(gser)) continue;
+
+ /* Omit unconvertable */
+ geos = POSTGIS2GEOS(gser);
+ if (!geos) continue;
+
+ geoms[ngeoms++] = geos;
+ }
+ array_free_iterator(iterator);
+
+ if (ngeoms == 0)
+ PG_RETURN_NULL();
+
+ geos = GEOSGeom_createCollection(
+ GEOS_GEOMETRYCOLLECTION,
+ geoms, ngeoms);
+
+ if (!geos)
+ {
+ coverage_destroy_geoms(geoms, ngeoms);
+ HANDLE_GEOS_ERROR("Geometry could not be converted");
+ }
+
+ geos_result = GEOSCoverageUnion(geos);
+ GEOSGeom_destroy(geos);
+ if (!geos_result)
+ HANDLE_GEOS_ERROR("Error computing coverage union");
+
+ result = GEOS2POSTGIS(geos_result, LW_FALSE);
+ GEOSGeom_destroy(geos_result);
+
+ PG_RETURN_POINTER(result);
+#endif
+}
+
+
diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
index 082201cf9..f30d75425 100644
--- a/postgis/postgis.sql.in
+++ b/postgis/postgis.sql.in
@@ -1911,6 +1911,7 @@ CREATE OR REPLACE FUNCTION ST_LineMerge(geometry, boolean)
LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
_COST_HIGH;
+
-----------------------------------------------------------------------------
-- Affine transforms
-----------------------------------------------------------------------------
@@ -4192,6 +4193,13 @@ CREATE OR REPLACE FUNCTION pgis_geometry_makeline_finalfn(internal)
LANGUAGE 'c' PARALLEL SAFE
_COST_MEDIUM;
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION pgis_geometry_coverageunion_finalfn(internal)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME'
+ LANGUAGE 'c' PARALLEL SAFE
+ _COST_MEDIUM;
+
-- Availability: 3.3.0
CREATE OR REPLACE FUNCTION pgis_geometry_union_parallel_transfn(internal, geometry)
RETURNS internal
@@ -4319,6 +4327,40 @@ CREATE AGGREGATE ST_MakeLine (geometry) (
FINALFUNC = pgis_geometry_makeline_finalfn
);
+
+-----------------------------------------------------------------------
+-- Polygonal coverage functions
+---------------------------------------------------------------------
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_CoverageUnion (geometry[])
+ RETURNS geometry
+ AS 'MODULE_PATHNAME', 'ST_CoverageUnion'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
+ _COST_MEDIUM;
+
+-- Availability: 3.4.0
+CREATE AGGREGATE ST_CoverageUnion (geometry) (
+ SFUNC = pgis_geometry_accum_transfn,
+ STYPE = internal,
+ PARALLEL = safe,
+ FINALFUNC = pgis_geometry_coverageunion_finalfn
+ );
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_CoverageSimplify (geom geometry, tolerance float8, simplifyBoundary boolean default true)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME', 'ST_CoverageSimplify'
+ LANGUAGE 'c' IMMUTABLE STRICT WINDOW PARALLEL SAFE
+ _COST_HIGH;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_CoverageInvalidLocations (geom geometry, tolerance float8 default 0.0)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME', 'ST_CoverageInvalidLocations'
+ LANGUAGE 'c' IMMUTABLE STRICT WINDOW PARALLEL SAFE
+ _COST_HIGH;
+
--------------------------------------------------------------------------------
-- Availability: 2.3.0
diff --git a/regress/core/coverage.sql b/regress/core/coverage.sql
new file mode 100644
index 000000000..287e66c75
--- /dev/null
+++ b/regress/core/coverage.sql
@@ -0,0 +1,103 @@
+
+CREATE TABLE coverage (id integer, seq integer, geom geometry);
+
+SELECT 'empty table a' AS test,
+ id, ST_AsText(GEOM) AS input,
+ ST_AsText(ST_CoverageInvalidLocations(geom) OVER (partition by id)) AS invalid,
+ ST_AsText(ST_CoverageSimplify(geom, 1.0) OVER (partition by id)) AS simple
+FROM coverage
+ORDER BY id, seq;
+
+SELECT 'empty table b' AS test,
+ id, ST_Area(ST_CoverageUnion(GEOM))
+FROM coverage GROUP BY id;
+
+INSERT INTO coverage VALUES
+(1, 1, 'MULTIPOLYGON(((0 0,10 0,10.1 5,10 10,0 10,0 0)))'),
+(1, 2, 'POLYGON((10 0,20 0,20 10,10 10,10.1 5,10 0))');
+
+SELECT 'one partition a' AS test,
+ id, seq, ST_AsText(GEOM) AS input,
+ ST_AsText(ST_CoverageInvalidLocations(geom) OVER (partition by id)) AS invalid,
+ ST_AsText(ST_CoverageSimplify(geom, 1.0) OVER (partition by id)) AS simple
+FROM coverage
+ORDER BY id, seq;
+
+-- Add another partition
+INSERT INTO coverage VALUES
+(2, 1, 'MULTIPOLYGON(((0 0,10 0,10.1 5,10 10,0 10,0 0)))'),
+(2, 2, 'POLYGON((10 0,20 0,20 10,10 10,10.1 5,10 0))'),
+(2, 3, NULL);
+
+SELECT 'two partition a' AS test,
+ id, seq, ST_AsText(geom) AS input,
+ ST_AsText(ST_CoverageInvalidLocations(geom) OVER (partition by id)) AS invalid,
+ ST_AsText(ST_CoverageSimplify(geom, 1.0) OVER (partition by id)) AS simple
+FROM coverage
+ORDER BY id, seq;
+
+-- Add another partition
+INSERT INTO coverage VALUES
+(3, 1, 'POLYGON EMPTY'),
+(3, 2, NULL);
+
+SELECT 'three partition a' AS test,
+ id, seq, ST_AsText(geom) AS input,
+ ST_AsText(ST_CoverageInvalidLocations(geom) OVER (partition by id)) AS invalid,
+ ST_AsText(ST_CoverageSimplify(geom, 1.0) OVER (partition by id)) AS simple
+FROM coverage
+ORDER BY id, seq;
+
+SELECT 'three partition b' AS test,
+ id, ST_Area(ST_CoverageUnion(geom))
+FROM coverage GROUP BY id ORDER BY id;
+
+
+TRUNCATE coverage;
+INSERT INTO coverage VALUES
+(4, 1, 'POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))'),
+(4, 2, 'POLYGON ((1 0, 0.9 1, 2 1, 2 0, 1 0))');
+
+WITH u AS (
+ SELECT ST_CoverageUnion(geom) AS geom FROM coverage
+)
+SELECT 'union 2' AS test,
+ ST_Area(geom), ST_GeometryType(geom)
+FROM u;
+
+
+TRUNCATE coverage;
+INSERT INTO coverage VALUES
+(5, 1, 'POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))'),
+(5, 2, 'POLYGON ((1 0, 1 1, 2 1, 2 0, 1 0))');
+
+WITH u AS (
+ SELECT ST_CoverageUnion(geom) AS geom FROM coverage
+)
+SELECT 'union 3' AS test,
+ ST_Area(geom), ST_GeometryType(geom)
+FROM u;
+
+DROP TABLE coverage;
+
+
+------------------------------------------------------------------------
+
+WITH squares AS (
+ SELECT *
+ FROM ST_SquareGrid(100.0, ST_MakeEnvelope(0, 0, 1000, 1000)) AS geom
+)
+SELECT 'grid lanes', i, ST_Area(ST_CoverageUnion(geom))
+FROM squares
+GROUP BY i
+ORDER By i;
+
+WITH squares AS (
+ SELECT i/5 AS i, j/5 AS j, geom
+ FROM ST_SquareGrid(100.0, ST_MakeEnvelope(0, 0, 1000, 1000)) AS geom
+)
+SELECT 'grid squares', i, j, ST_Area(ST_CoverageUnion(geom))
+FROM squares
+GROUP BY i, j
+ORDER By i, j;
+
diff --git a/regress/core/coverage_expected b/regress/core/coverage_expected
new file mode 100644
index 000000000..5d19e968a
--- /dev/null
+++ b/regress/core/coverage_expected
@@ -0,0 +1,39 @@
+one partition a|1|1|MULTIPOLYGON(((0 0,10 0,10.1 5,10 10,0 10,0 0)))||MULTIPOLYGON(((10 0,10 10,0 10,0 0,10 0)))
+one partition a|1|2|POLYGON((10 0,20 0,20 10,10 10,10.1 5,10 0))||POLYGON((10 0,20 0,20 10,10 10,10 0))
+two partition a|1|1|MULTIPOLYGON(((0 0,10 0,10.1 5,10 10,0 10,0 0)))||MULTIPOLYGON(((10 0,10 10,0 10,0 0,10 0)))
+two partition a|1|2|POLYGON((10 0,20 0,20 10,10 10,10.1 5,10 0))||POLYGON((10 0,20 0,20 10,10 10,10 0))
+two partition a|2|1|MULTIPOLYGON(((0 0,10 0,10.1 5,10 10,0 10,0 0)))||MULTIPOLYGON(((10 0,10 10,0 10,0 0,10 0)))
+two partition a|2|2|POLYGON((10 0,20 0,20 10,10 10,10.1 5,10 0))||POLYGON((10 0,20 0,20 10,10 10,10 0))
+two partition a|2|3|||
+three partition a|1|1|MULTIPOLYGON(((0 0,10 0,10.1 5,10 10,0 10,0 0)))||MULTIPOLYGON(((10 0,10 10,0 10,0 0,10 0)))
+three partition a|1|2|POLYGON((10 0,20 0,20 10,10 10,10.1 5,10 0))||POLYGON((10 0,20 0,20 10,10 10,10 0))
+three partition a|2|1|MULTIPOLYGON(((0 0,10 0,10.1 5,10 10,0 10,0 0)))||MULTIPOLYGON(((10 0,10 10,0 10,0 0,10 0)))
+three partition a|2|2|POLYGON((10 0,20 0,20 10,10 10,10.1 5,10 0))||POLYGON((10 0,20 0,20 10,10 10,10 0))
+three partition a|2|3|||
+three partition a|3|1|POLYGON EMPTY||
+three partition a|3|2|||
+three partition b|1|200
+three partition b|2|200
+three partition b|3|
+union 2|2.05|ST_MultiPolygon
+union 3|2|ST_Polygon
+grid lanes|0|110000
+grid lanes|1|110000
+grid lanes|2|110000
+grid lanes|3|110000
+grid lanes|4|110000
+grid lanes|5|110000
+grid lanes|6|110000
+grid lanes|7|110000
+grid lanes|8|110000
+grid lanes|9|110000
+grid lanes|10|110000
+grid squares|0|0|250000
+grid squares|0|1|250000
+grid squares|0|2|50000
+grid squares|1|0|250000
+grid squares|1|1|250000
+grid squares|1|2|50000
+grid squares|2|0|50000
+grid squares|2|1|50000
+grid squares|2|2|10000
diff --git a/regress/core/tests.mk.in b/regress/core/tests.mk.in
index 60bd6d23c..e1fc968df 100644
--- a/regress/core/tests.mk.in
+++ b/regress/core/tests.mk.in
@@ -179,6 +179,11 @@ else
$(top_srcdir)/regress/core/concave_hull
endif
+ifeq ($(shell expr "$(POSTGIS_GEOS_VERSION)" ">=" 31200),1)
+ TESTS += \
+ $(top_srcdir)/regress/core/coverage
+endif
+
ifeq ($(INTERRUPTTESTS),yes)
# Allow CI servers to configure --with-interrupt-tests
TESTS += \
-----------------------------------------------------------------------
Summary of changes:
doc/Makefile.in | 1 +
doc/postgis.xml | 1 +
doc/reference.xml | 1 +
doc/reference_coverage.xml | 195 +++++++++++++++++++
liblwgeom/lwgeom_geos.c | 8 +-
liblwgeom/lwgeom_topo.c | 2 +-
postgis/lwgeom_accum.c | 29 +++
postgis/lwgeom_geos.c | 18 +-
postgis/lwgeom_geos.h | 16 ++
postgis/lwgeom_window.c | 430 ++++++++++++++++++++++++++++++++++++++++-
postgis/postgis.sql.in | 42 ++++
regress/core/coverage.sql | 103 ++++++++++
regress/core/coverage_expected | 39 ++++
regress/core/tests.mk.in | 5 +
14 files changed, 865 insertions(+), 25 deletions(-)
create mode 100644 doc/reference_coverage.xml
create mode 100644 regress/core/coverage.sql
create mode 100644 regress/core/coverage_expected
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list