[postgis-tickets] [SCM] PostGIS branch master updated. 3.1.0alpha2-65-gc97b84f

git at osgeo.org git at osgeo.org
Fri Sep 11 01:54:54 PDT 2020


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  c97b84fb39976897dc1ef621eab5cc6dc1187efe (commit)
      from  a15f8921b434454f9264eb4de0dbe12ff79a0711 (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 c97b84fb39976897dc1ef621eab5cc6dc1187efe
Author: Sandro Santilli <strk at kbt.io>
Date:   Mon Aug 24 22:40:32 2020 +0200

    Expose fixed-precision overlay functions
    
    Adds an optional gridSize argument to:
      - ST_Intersection
      - ST_Difference
      - ST_SymDifference
      - ST_Union (including aggregate)
      - ST_UnaryUnion
      - ST_Subdivide

diff --git a/doc/reference_processing.xml b/doc/reference_processing.xml
index 205a872..024ab31 100644
--- a/doc/reference_processing.xml
+++ b/doc/reference_processing.xml
@@ -1255,6 +1255,7 @@ GEOMETRYCOLLECTION Z (POLYGON Z ((14 14 10,20 10 150,34 6 25,14 14 10))
             <funcdef>geometry <function>ST_Difference</function></funcdef>
             <paramdef><type>geometry </type> <parameter>geomA</parameter></paramdef>
             <paramdef><type>geometry </type> <parameter>geomB</parameter></paramdef>
+            <paramdef choice="opt"><type>float8 </type> <parameter>gridSize</parameter></paramdef>
           </funcprototype>
         </funcsynopsis>
       </refsynopsisdiv>
@@ -1267,8 +1268,16 @@ GEOMETRYCOLLECTION Z (POLYGON Z ((14 14 10,20 10 150,34 6 25,14 14 10))
             then an empty geometry collection is returned.</para>
         <note><para>Order matters. B - A will always return a portion of B</para></note>
 
+            <para>
+If a gridSize argument is provided, the inputs are
+snapped to a grid of the given size and all intersections are computed
+on that same grid. This method requires GEOS-3.9.0 or higher.
+            </para>
+
         <para>Performed by the GEOS module</para>
 
+        <para>Enhanced: 3.1.0 accept a gridSize parameter.</para>
+
         <para>&sfs_compliant; s2.1.1.3</para>
         <para>&sqlmm_compliant; SQL-MM 3: 5.1.20</para>
         <para>&Z_support; However it seems to only consider x y when
@@ -1631,6 +1640,10 @@ Returns a geometry that represents the shared portion of geomA and geomB.
                         <type>geometry</type>
                         <parameter>geomB</parameter>
                     </paramdef>
+                    <paramdef choice="opt">
+                        <type>float8</type>
+                        <parameter>gridSize</parameter>
+                    </paramdef>
                 </funcprototype>
                 <funcprototype>
                     <funcdef>geography <function>ST_Intersection</function></funcdef>
@@ -1658,6 +1671,12 @@ Returns a geometry that represents the shared portion of geomA and geomB.
             <para>ST_Intersection in conjunction with ST_Intersects is very useful for clipping geometries such as in bounding box, buffer, region
                 queries where you only want to return that portion of a geometry that sits in a country or region of interest.</para>
 
+            <para>
+If a gridSize argument is provided, the inputs are
+snapped to a grid of the given size and all intersections are computed
+on that same grid. This method requires GEOS-3.9.0 or higher.
+            </para>
+
             <note><para>Geography: For geography this is really a thin wrapper around the geometry implementation. It first determines the best SRID that
                     fits the bounding box of the 2 geography objects (if geography objects are within one half zone UTM but not same UTM will pick one of those) (favoring UTM or Lambert Azimuthal Equal Area (LAEA) north/south pole, and falling back on mercator in worst case scenario)  and then intersection in that best fit planar spatial ref and retransforms back to WGS84 geography.</para></note>
 
@@ -1669,6 +1688,7 @@ Returns a geometry that represents the shared portion of geomA and geomB.
 
           <para>Availability: 1.5 support for geography data type was introduced.</para>
             <para>Changed: 3.0.0 does not depend on SFCGAL.</para>
+            <para>Changed: 3.1.0 added gridSize argument.</para>
 
           <para>&sfs_compliant; s2.1.1.3</para>
           <para>&sqlmm_compliant; SQL-MM 3: 5.1.18</para>
@@ -3380,6 +3400,7 @@ GEOMETRYCOLLECTION(
             <funcdef>geometry <function>ST_SymDifference</function></funcdef>
             <paramdef><type>geometry </type> <parameter>geomA</parameter></paramdef>
             <paramdef><type>geometry </type> <parameter>geomB</parameter></paramdef>
+            <paramdef choice="opt"><type>float8 </type> <parameter>gridSize</parameter></paramdef>
           </funcprototype>
         </funcsynopsis>
       </refsynopsisdiv>
@@ -3392,7 +3413,16 @@ GEOMETRYCOLLECTION(
             ST_SymDifference(A,B) = ST_SymDifference(B,A). One can think of this as ST_Union(geomA,geomB) - ST_Intersection(A,B).
             </para>
 
+            <para>
+If a gridSize argument is provided, the inputs are
+snapped to a grid of the given size and all intersections are computed
+on that same grid. This method requires GEOS-3.9.0 or higher.
+            </para>
+
         <para>Performed by the GEOS module</para>
+
+        <para>Enhanced: 3.1.0 accept a gridSize parameter.</para>
+
         <para>&sfs_compliant; s2.1.1.3</para>
         <para>&sqlmm_compliant; SQL-MM 3: 5.1.21</para>
         <para>&Z_support; However it seems to only consider x y when
@@ -3482,6 +3512,7 @@ MULTILINESTRING((1 3 2.75,1 4 2),(1 1 3,1 2 2.25))
                 <funcdef>setof geometry <function>ST_Subdivide</function></funcdef>
                 <paramdef><type>geometry</type> <parameter>geom</parameter></paramdef>
                 <paramdef><type>integer</type> <parameter>max_vertices=256</parameter></paramdef>
+                <paramdef choice="opt"><type>float8</type> <parameter>gridSize</parameter></paramdef>
             </funcprototype>
         </funcsynopsis>
     </refsynopsisdiv>
@@ -3494,12 +3525,13 @@ MULTILINESTRING((1 3 2.75,1 4 2),(1 1 3,1 2 2.25))
             Point-in-polygon and other overlay operations are normally faster for indexed subdivided dataset:
             "miss" cases are faster to check as boxes for all parts typically cover smaller area than original geometry box,
             "hit" cases are faster because recheck operates on less points.
-            Uses the same envelope clipping as <code>ST_ClipByBox2D</code>.
             <code>max_vertices</code> must be 5 or more, as 5 points are needed to represent a closed box.
+            <code>gridSize</code> can be specified to have clipping work in fixed-precision space (requires GEOS-3.9.0+).
         </para>
         <para>Performed by the GEOS module.</para>
         <para>Availability: 2.2.0</para>
         <para>Enhanced: 2.5.0 reuses existing points on polygon split, vertex count is lowered from 8 to 5.</para>
+        <para>Enhanced: 3.1.0 accept a gridSize parameter.</para>
       </refsection>
 
       <refsection>
@@ -3601,6 +3633,12 @@ LINESTRING(44.7994523421035 82.5156766227011,85 85)</screen>
       </funcprototype>
       <funcprototype>
         <funcdef>geometry <function>ST_Union</function></funcdef>
+        <paramdef><type>geometry</type> <parameter>g1</parameter></paramdef>
+        <paramdef><type>geometry</type> <parameter>g2</parameter></paramdef>
+        <paramdef><type>float8</type> <parameter>gridSize</parameter></paramdef>
+      </funcprototype>
+      <funcprototype>
+        <funcdef>geometry <function>ST_Union</function></funcdef>
         <paramdef><type>geometry[]</type> <parameter>g1_array</parameter></paramdef>
       </funcprototype>
     </funcsynopsis>
@@ -3622,6 +3660,11 @@ LINESTRING(44.7994523421035 82.5156766227011,85 85)</screen>
         input geometries. Output type can be a MULTI*, NON-MULTI or
         GEOMETRYCOLLECTION. If any are NULL, then NULL is returned.</para>
 
+    <para>
+        For non-aggregate version, <code>gridSize</code> can be specified
+        to work in fixed-precision space (requires GEOS-3.9.0+).
+    </para>
+
     <note><para>ST_Collect and ST_Union are often interchangeable.
         ST_Union is in general orders of magnitude slower than ST_Collect
         because it tries to dissolve boundaries and reorder geometries to ensure that a constructed Multi* doesn't
@@ -3636,6 +3679,7 @@ LINESTRING(44.7994523421035 82.5156766227011,85 85)</screen>
         word.</para>
     <para>Availability: 1.4.0 - ST_Union was enhanced. ST_Union(geomarray) was introduced and also faster aggregate collection in PostgreSQL.</para>
     <para>Changed: 3.0.0 does not depend on SFCGAL.</para>
+    <para>Enhanced: 3.1.0 accept a gridSize parameter.</para>
 
     <para>&sfs_compliant; s2.1.1.3</para>
     <note><para>Aggregate version is not explicitly defined in OGC SPEC.</para></note>
@@ -3732,6 +3776,7 @@ MULTILINESTRING((3 4,4 5),(1 2,3 4))
           <funcprototype>
             <funcdef>geometry <function>ST_UnaryUnion</function></funcdef>
             <paramdef><type>geometry </type> <parameter>geom</parameter></paramdef>
+            <paramdef choice="opt"><type>float8 </type> <parameter>gridSize</parameter></paramdef>
           </funcprototype>
 
         </funcsynopsis>
@@ -3758,9 +3803,16 @@ MULTILINESTRING((3 4,4 5),(1 2,3 4))
         balance between ST_Union and ST_MemUnion.
         </para>
 
+        <para>
+If a gridSize argument is provided, the inputs are
+snapped to a grid of the given size and all intersections are computed
+on that same grid. This method requires GEOS-3.9.0 or higher.
+        </para>
+
         <para>&Z_support;</para>
 
         <para>Availability: 2.0.0</para>
+        <para>Enhanced: 3.1.0 accept a gridSize parameter.</para>
       </refsection>
 
 
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index 59a94ce..f44588f 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -2364,15 +2364,21 @@ LWGEOM* lwgeom_geos_noop(const LWGEOM *geom) ;
 
 LWGEOM *lwgeom_normalize(const LWGEOM *geom);
 LWGEOM *lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2);
+LWGEOM *lwgeom_intersection_prec(const LWGEOM *geom1, const LWGEOM *geom2, double gridSize);
 LWGEOM *lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2);
+LWGEOM *lwgeom_difference_prec(const LWGEOM *geom1, const LWGEOM *geom2, double gridSize);
 LWGEOM *lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2);
+LWGEOM *lwgeom_symdifference_prec(const LWGEOM* geom1, const LWGEOM* geom2, double gridSize);
 LWGEOM *lwgeom_pointonsurface(const LWGEOM* geom);
 LWGEOM *lwgeom_centroid(const LWGEOM* geom);
 LWGEOM *lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2);
+LWGEOM *lwgeom_union_prec(const LWGEOM *geom1, const LWGEOM *geom2, double gridSize);
 LWGEOM *lwgeom_linemerge(const LWGEOM *geom1);
 LWGEOM *lwgeom_unaryunion(const LWGEOM *geom1);
+LWGEOM *lwgeom_unaryunion_prec(const LWGEOM *geom1, double gridSize);
 LWGEOM *lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double y1);
 LWCOLLECTION *lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices);
+LWCOLLECTION *lwgeom_subdivide_prec(const LWGEOM *geom, uint32_t maxvertices, double gridSize);
 
 /**
  * Snap vertices and segments of a geometry to another using a given tolerance.
diff --git a/liblwgeom/lwgeom.c b/liblwgeom/lwgeom.c
index d8ee087..1ca1bd5 100644
--- a/liblwgeom/lwgeom.c
+++ b/liblwgeom/lwgeom.c
@@ -2262,14 +2262,16 @@ static void lwgeom_subdivide_recursive(const LWGEOM *geom,
 				       uint8_t dimension,
 				       uint32_t maxvertices,
 				       uint32_t depth,
-				       LWCOLLECTION *col);
+				       LWCOLLECTION *col,
+				       double gridSize);
 
 static void
 lwgeom_subdivide_recursive(const LWGEOM *geom,
 			   uint8_t dimension,
 			   uint32_t maxvertices,
 			   uint32_t depth,
-			   LWCOLLECTION *col)
+			   LWCOLLECTION *col,
+			   double gridSize)
 {
 	const uint32_t maxdepth = 50;
 
@@ -2317,7 +2319,7 @@ lwgeom_subdivide_recursive(const LWGEOM *geom,
 		/* Don't increment depth yet, since we aren't actually
 		 * subdividing geometries yet */
 		for (uint32_t i = 0; i < incol->ngeoms; i++ )
-			lwgeom_subdivide_recursive(incol->geoms[i], dimension, maxvertices, depth, col);
+			lwgeom_subdivide_recursive(incol->geoms[i], dimension, maxvertices, depth, col, gridSize);
 		return;
 	}
 
@@ -2421,31 +2423,31 @@ lwgeom_subdivide_recursive(const LWGEOM *geom,
 	{
 		LWGEOM *subbox = (LWGEOM *)lwpoly_construct_envelope(
 		    geom->srid, subbox1.xmin, subbox1.ymin, subbox1.xmax, subbox1.ymax);
-		LWGEOM *clipped = lwgeom_intersection(geom, subbox);
+		LWGEOM *clipped = lwgeom_intersection_prec(geom, subbox, gridSize);
 		lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
 		lwgeom_free(subbox);
 		if (clipped && !lwgeom_is_empty(clipped))
 		{
-			lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
+			lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col, gridSize);
 			lwgeom_free(clipped);
 		}
 	}
 	{
 		LWGEOM *subbox = (LWGEOM *)lwpoly_construct_envelope(
 		    geom->srid, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
-		LWGEOM *clipped = lwgeom_intersection(geom, subbox);
+		LWGEOM *clipped = lwgeom_intersection_prec(geom, subbox, gridSize);
 		lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
 		lwgeom_free(subbox);
 		if (clipped && !lwgeom_is_empty(clipped))
 		{
-			lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
+			lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col, gridSize);
 			lwgeom_free(clipped);
 		}
 	}
 }
 
 LWCOLLECTION *
-lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices)
+lwgeom_subdivide_prec(const LWGEOM *geom, uint32_t maxvertices, double gridSize)
 {
 	static uint32_t startdepth = 0;
 	static uint32_t minmaxvertices = 5;
@@ -2462,11 +2464,19 @@ lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices)
 		lwerror("%s: cannot subdivide to fewer than %d vertices per output", __func__, minmaxvertices);
 	}
 
-	lwgeom_subdivide_recursive(geom, lwgeom_dimension(geom), maxvertices, startdepth, col);
+	lwgeom_subdivide_recursive(geom, lwgeom_dimension(geom), maxvertices, startdepth, col, gridSize);
 	lwgeom_set_srid((LWGEOM*)col, geom->srid);
 	return col;
 }
 
+LWCOLLECTION *
+lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices)
+{
+	return lwgeom_subdivide_prec(geom, maxvertices, -1);
+}
+
+
+
 int
 lwgeom_is_trajectory(const LWGEOM *geom)
 {
diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c
index f735f70..1ca3f4f 100644
--- a/liblwgeom/lwgeom_geos.c
+++ b/liblwgeom/lwgeom_geos.c
@@ -18,7 +18,7 @@
  *
  **********************************************************************
  *
- * Copyright 2011-2014 Sandro Santilli <strk at kbt.io>
+ * Copyright 2011-2020 Sandro Santilli <strk at kbt.io>
  * Copyright 2015-2018 Daniel Baston <dbaston at gmail.com>
  * Copyright 2017-2018 Darafei Praliaskouski <me at komzpa.net>
  *
@@ -673,7 +673,13 @@ lwgeom_normalize(const LWGEOM* geom)
 }
 
 LWGEOM*
-lwgeom_intersection(const LWGEOM* geom1, const LWGEOM* geom2)
+lwgeom_intersection(const LWGEOM* g1, const LWGEOM* g2)
+{
+	return lwgeom_intersection_prec(g1, g2, -1.0);
+}
+
+LWGEOM*
+lwgeom_intersection_prec(const LWGEOM* geom1, const LWGEOM* geom2, double prec)
 {
 	LWGEOM* result;
 	int32_t srid = RESULT_SRID(geom1, geom2);
@@ -695,7 +701,19 @@ lwgeom_intersection(const LWGEOM* geom1, const LWGEOM* geom2)
 	if (!(g1 = LWGEOM2GEOS(geom1, AUTOFIX))) GEOS_FAIL();
 	if (!(g2 = LWGEOM2GEOS(geom2, AUTOFIX))) GEOS_FREE_AND_FAIL(g1);
 
-	g3 = GEOSIntersection(g1, g2);
+	if ( prec >= 0) {
+#if POSTGIS_GEOS_VERSION < 39
+		lwerror("Fixed-precision intersection requires GEOS-3.9 or higher");
+		GEOS_FREE_AND_FAIL(g1, g2);
+		return NULL;
+#else
+		g3 = GEOSIntersectionPrec(g1, g2, prec);
+#endif
+	}
+	else
+	{
+		g3 = GEOSIntersection(g1, g2);
+	}
 
 	if (!g3) GEOS_FREE_AND_FAIL(g1);
 	GEOSSetSRID(g3, srid);
@@ -740,6 +758,12 @@ lwgeom_linemerge(const LWGEOM* geom)
 LWGEOM*
 lwgeom_unaryunion(const LWGEOM* geom)
 {
+	return lwgeom_unaryunion_prec(geom, -1.0);
+}
+
+LWGEOM*
+lwgeom_unaryunion_prec(const LWGEOM* geom, double prec)
+{
 	LWGEOM* result;
 	int32_t srid = RESULT_SRID(geom);
 	uint8_t is3d = FLAGS_GET_Z(geom->flags);
@@ -755,7 +779,19 @@ lwgeom_unaryunion(const LWGEOM* geom)
 
 	if (!(g1 = LWGEOM2GEOS(geom, AUTOFIX))) GEOS_FAIL();
 
-	g3 = GEOSUnaryUnion(g1);
+	if ( prec >= 0) {
+#if POSTGIS_GEOS_VERSION < 39
+		lwerror("Fixed-precision union requires GEOS-3.9 or higher");
+		GEOS_FREE_AND_FAIL(g1);
+		return NULL;
+#else
+		g3 = GEOSUnaryUnionPrec(g1, prec);
+#endif
+	}
+	else
+	{
+		g3 = GEOSUnaryUnion(g1);
+	}
 
 	if (!g3) GEOS_FREE_AND_FAIL(g1);
 	GEOSSetSRID(g3, srid);
@@ -771,6 +807,12 @@ lwgeom_unaryunion(const LWGEOM* geom)
 LWGEOM*
 lwgeom_difference(const LWGEOM* geom1, const LWGEOM* geom2)
 {
+	return lwgeom_difference_prec(geom1, geom2, -1.0);
+}
+
+LWGEOM*
+lwgeom_difference_prec(const LWGEOM* geom1, const LWGEOM* geom2, double prec)
+{
 	LWGEOM* result;
 	int32_t srid = RESULT_SRID(geom1, geom2);
 	uint8_t is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags));
@@ -789,7 +831,19 @@ lwgeom_difference(const LWGEOM* geom1, const LWGEOM* geom2)
 	if (!(g1 = LWGEOM2GEOS(geom1, AUTOFIX))) GEOS_FAIL();
 	if (!(g2 = LWGEOM2GEOS(geom2, AUTOFIX))) GEOS_FREE_AND_FAIL(g1);
 
-	g3 = GEOSDifference(g1, g2);
+	if ( prec >= 0) {
+#if POSTGIS_GEOS_VERSION < 39
+		lwerror("Fixed-precision difference requires GEOS-3.9 or higher");
+		GEOS_FREE_AND_FAIL(g1, g2);
+		return NULL;
+#else
+		g3 = GEOSDifferencePrec(g1, g2, prec);
+#endif
+	}
+	else
+	{
+		g3 = GEOSDifference(g1, g2);
+	}
 
 	if (!g3) GEOS_FREE_AND_FAIL(g1, g2);
 	GEOSSetSRID(g3, srid);
@@ -804,6 +858,12 @@ lwgeom_difference(const LWGEOM* geom1, const LWGEOM* geom2)
 LWGEOM*
 lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2)
 {
+	return lwgeom_symdifference_prec(geom1, geom2, -1.0);
+}
+
+LWGEOM*
+lwgeom_symdifference_prec(const LWGEOM* geom1, const LWGEOM* geom2, double prec)
+{
 	LWGEOM* result;
 	int32_t srid = RESULT_SRID(geom1, geom2);
 	uint8_t is3d = (FLAGS_GET_Z(geom1->flags) || FLAGS_GET_Z(geom2->flags));
@@ -822,7 +882,19 @@ lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2)
 	if (!(g1 = LWGEOM2GEOS(geom1, AUTOFIX))) GEOS_FAIL();
 	if (!(g2 = LWGEOM2GEOS(geom2, AUTOFIX))) GEOS_FREE_AND_FAIL(g1);
 
-	g3 = GEOSSymDifference(g1, g2);
+	if ( prec >= 0) {
+#if POSTGIS_GEOS_VERSION < 39
+		lwerror("Fixed-precision difference requires GEOS-3.9 or higher");
+		GEOS_FREE_AND_FAIL(g1, g2);
+		return NULL;
+#else
+		g3 = GEOSSymDifferencePrec(g1, g2, prec);
+#endif
+	}
+	else
+	{
+		g3 = GEOSSymDifference(g1, g2);
+	}
 
 	if (!g3) GEOS_FREE_AND_FAIL(g1, g2);
 	GEOSSetSRID(g3, srid);
@@ -901,7 +973,13 @@ lwgeom_pointonsurface(const LWGEOM *geom)
 }
 
 LWGEOM*
-lwgeom_union(const LWGEOM* geom1, const LWGEOM* geom2)
+lwgeom_union(const LWGEOM* g1, const LWGEOM* g2)
+{
+	return lwgeom_union_prec(g1, g2, -1.0);
+}
+
+LWGEOM*
+lwgeom_union_prec(const LWGEOM* geom1, const LWGEOM* geom2, double gridSize)
 {
 	LWGEOM* result;
 	int32_t srid = RESULT_SRID(geom1, geom2);
@@ -921,7 +999,19 @@ lwgeom_union(const LWGEOM* geom1, const LWGEOM* geom2)
 	if (!(g1 = LWGEOM2GEOS(geom1, AUTOFIX))) GEOS_FAIL();
 	if (!(g2 = LWGEOM2GEOS(geom2, AUTOFIX))) GEOS_FREE_AND_FAIL(g1);
 
-	g3 = GEOSUnion(g1, g2);
+	if ( gridSize >= 0) {
+#if POSTGIS_GEOS_VERSION < 39
+		lwerror("Fixed-precision union requires GEOS-3.9 or higher");
+		GEOS_FREE_AND_FAIL(g1, g2);
+		return NULL;
+#else
+		g3 = GEOSUnionPrec(g1, g2, gridSize);
+#endif
+	}
+	else
+	{
+		g3 = GEOSUnion(g1, g2);
+	}
 
 	if (!g3) GEOS_FREE_AND_FAIL(g1, g2);
 	GEOSSetSRID(g3, srid);
diff --git a/postgis/legacy.sql.in b/postgis/legacy.sql.in
index d750541..c90452a 100644
--- a/postgis/legacy.sql.in
+++ b/postgis/legacy.sql.in
@@ -1806,13 +1806,13 @@ CREATE OR REPLACE FUNCTION StartPoint(geometry)
 -- Deprecation in 1.2.3
 CREATE OR REPLACE FUNCTION symdifference(geometry,geometry)
 	RETURNS geometry
-	AS 'MODULE_PATHNAME','symdifference'
+	AS 'MODULE_PATHNAME','ST_SymDifference'
 	LANGUAGE 'c' IMMUTABLE STRICT;
 
 -- Deprecation in 1.2.3
 CREATE OR REPLACE FUNCTION symmetricdifference(geometry,geometry)
 	RETURNS geometry
-	AS 'MODULE_PATHNAME','symdifference'
+	AS 'MODULE_PATHNAME','ST_SymDifference'
 	LANGUAGE 'c' IMMUTABLE STRICT;
 
 -- Deprecation in 1.2.3
diff --git a/postgis/lwgeom_accum.c b/postgis/lwgeom_accum.c
index eb96fe4..4b4a688 100644
--- a/postgis/lwgeom_accum.c
+++ b/postgis/lwgeom_accum.c
@@ -72,6 +72,7 @@ pgis_geometry_accum_transfn(PG_FUNCTION_ARGS)
 	LWGEOM *geom = NULL;
 	GSERIALIZED *gser = NULL;
 	Datum argType = get_fn_expr_argtype(fcinfo->flinfo, 1);
+	double gridSize = -1.0;
 
 	if (argType == InvalidOid)
 		ereport(ERROR,
@@ -92,6 +93,7 @@ pgis_geometry_accum_transfn(PG_FUNCTION_ARGS)
 		state = MemoryContextAlloc(aggcontext, sizeof(CollectionBuildState));
 		state->geoms = NULL;
 		state->geomOid = argType;
+		state->gridSize = gridSize;
 
 		for (int i = 0; i < n; i++)
 		{
@@ -110,6 +112,13 @@ pgis_geometry_accum_transfn(PG_FUNCTION_ARGS)
 	if (!PG_ARGISNULL(1))
 		gser = PG_GETARG_GSERIALIZED_P(1);
 
+	if (PG_NARGS()>2 && !PG_ARGISNULL(2))
+	{
+		gridSize = PG_GETARG_FLOAT8(2);
+		/*lwnotice("Passed gridSize %g", gridSize);*/
+		if ( gridSize > state->gridSize ) state->gridSize = gridSize;
+	}
+
 	/* Take a copy of the geometry into the aggregate context */
 	old = MemoryContextSwitchTo(aggcontext);
 	if (gser)
diff --git a/postgis/lwgeom_accum.h b/postgis/lwgeom_accum.h
index c500d67..e1da21b 100644
--- a/postgis/lwgeom_accum.h
+++ b/postgis/lwgeom_accum.h
@@ -38,6 +38,7 @@ typedef struct CollectionBuildState
 	List *geoms;  /* collected geometries */
 	Datum data[CollectionBuildStateDataSize];
 	Oid geomOid;
+	float8 gridSize;
 } CollectionBuildState;
 
 
diff --git a/postgis/lwgeom_dump.c b/postgis/lwgeom_dump.c
index 7d6ae64..e096bf3 100644
--- a/postgis/lwgeom_dump.c
+++ b/postgis/lwgeom_dump.c
@@ -347,6 +347,7 @@ Datum ST_Subdivide(PG_FUNCTION_ARGS)
 		LWCOLLECTION *col;
 		/* default to maxvertices < page size */
 		int maxvertices = 128;
+		double gridSize = -1;
 
 		/* create a function context for cross-call persistence */
 		funcctx = SRF_FIRSTCALL_INIT();
@@ -369,9 +370,15 @@ Datum ST_Subdivide(PG_FUNCTION_ARGS)
 			maxvertices = PG_GETARG_INT32(1);
 
 		/*
+		* Get the gridSize value
+		*/
+		if ( PG_NARGS() > 2 && ! PG_ARGISNULL(2) )
+			gridSize = PG_GETARG_FLOAT8(2);
+
+		/*
 		* Compute the subdivision of the geometry
 		*/
-		col = lwgeom_subdivide(geom, maxvertices);
+		col = lwgeom_subdivide_prec(geom, maxvertices, gridSize);
 
 		if ( ! col )
 			SRF_RETURN_DONE(funcctx);
diff --git a/postgis/lwgeom_geos.c b/postgis/lwgeom_geos.c
index e7fd51a..ed5d27e 100644
--- a/postgis/lwgeom_geos.c
+++ b/postgis/lwgeom_geos.c
@@ -83,7 +83,7 @@ Datum convexhull(PG_FUNCTION_ARGS);
 Datum topologypreservesimplify(PG_FUNCTION_ARGS);
 Datum ST_Difference(PG_FUNCTION_ARGS);
 Datum boundary(PG_FUNCTION_ARGS);
-Datum symdifference(PG_FUNCTION_ARGS);
+Datum ST_SymDifference(PG_FUNCTION_ARGS);
 Datum ST_Union(PG_FUNCTION_ARGS);
 Datum issimple(PG_FUNCTION_ARGS);
 Datum isring(PG_FUNCTION_ARGS);
@@ -683,23 +683,13 @@ Datum pgis_geometry_union_finalfn(PG_FUNCTION_ARGS)
 	*/
 	if (ngeoms > 0)
 	{
-		GEOSGeometry *g = NULL;
-		GEOSGeometry *g_union = NULL;
 		LWCOLLECTION* col = lwcollection_construct(COLLECTIONTYPE, srid, NULL, ngeoms, geoms);
-
-		initGEOS(lwpgnotice, lwgeom_geos_error);
-		g = LWGEOM2GEOS((LWGEOM*)col, LW_FALSE);
-		if (!g)
-			HANDLE_GEOS_ERROR("Could not create GEOS COLLECTION from geometry array");
-
-		g_union = GEOSUnaryUnion(g);
-		GEOSGeom_destroy(g);
-		if (!g_union)
-			HANDLE_GEOS_ERROR("GEOSUnaryUnion");
-
-		GEOSSetSRID(g_union, srid);
-		gser_out = GEOS2POSTGIS(g_union, has_z);
-		GEOSGeom_destroy(g_union);
+		LWGEOM *out = lwgeom_unaryunion_prec(lwcollection_as_lwgeom(col), state->gridSize);
+		if ( ! out )
+		{
+			lwcollection_free(col);
+		}
+		gser_out = geometry_serialize(out);
 	}
 	/* No real geometries in our array, any empties? */
 	else
@@ -737,12 +727,15 @@ Datum ST_UnaryUnion(PG_FUNCTION_ARGS)
 	GSERIALIZED *geom1;
 	GSERIALIZED *result;
 	LWGEOM *lwgeom1, *lwresult ;
+	double prec = -1;
 
 	geom1 = PG_GETARG_GSERIALIZED_P(0);
+	if (PG_NARGS() > 1 && ! PG_ARGISNULL(1))
+		prec = PG_GETARG_FLOAT8(1);
 
 	lwgeom1 = lwgeom_from_gserialized(geom1) ;
 
-	lwresult = lwgeom_unaryunion(lwgeom1);
+	lwresult = lwgeom_unaryunion_prec(lwgeom1, prec);
 	result = geometry_serialize(lwresult) ;
 
 	lwgeom_free(lwgeom1) ;
@@ -760,14 +753,17 @@ Datum ST_Union(PG_FUNCTION_ARGS)
 	GSERIALIZED *geom2;
 	GSERIALIZED *result;
 	LWGEOM *lwgeom1, *lwgeom2, *lwresult;
+	double gridSize = -1;
 
 	geom1 = PG_GETARG_GSERIALIZED_P(0);
 	geom2 = PG_GETARG_GSERIALIZED_P(1);
+	if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
+		gridSize = PG_GETARG_FLOAT8(2);
 
 	lwgeom1 = lwgeom_from_gserialized(geom1);
 	lwgeom2 = lwgeom_from_gserialized(geom2);
 
-	lwresult = lwgeom_union(lwgeom1, lwgeom2);
+	lwresult = lwgeom_union_prec(lwgeom1, lwgeom2, gridSize);
 	result = geometry_serialize(lwresult);
 
 	lwgeom_free(lwgeom1);
@@ -780,26 +776,40 @@ Datum ST_Union(PG_FUNCTION_ARGS)
 	PG_RETURN_POINTER(result);
 }
 
+/* This is retained for backward ABI compatibility
+ * with PostGIS < 3.1.0 */
+PG_FUNCTION_INFO_V1(symdifference);
+Datum symdifference(PG_FUNCTION_ARGS)
+{
+  PG_RETURN_DATUM(DirectFunctionCall2(
+     ST_SymDifference,
+     PG_GETARG_DATUM(0), PG_GETARG_DATUM(1)
+  ));
+}
+
 /**
  *  @example symdifference {@link #symdifference} - SELECT ST_SymDifference(
  *      'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))',
  *      'POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))');
  */
-PG_FUNCTION_INFO_V1(symdifference);
-Datum symdifference(PG_FUNCTION_ARGS)
+PG_FUNCTION_INFO_V1(ST_SymDifference);
+Datum ST_SymDifference(PG_FUNCTION_ARGS)
 {
 	GSERIALIZED *geom1;
 	GSERIALIZED *geom2;
 	GSERIALIZED *result;
 	LWGEOM *lwgeom1, *lwgeom2, *lwresult ;
+	double prec = -1;
 
 	geom1 = PG_GETARG_GSERIALIZED_P(0);
 	geom2 = PG_GETARG_GSERIALIZED_P(1);
+	if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
+		prec = PG_GETARG_FLOAT8(2);
 
 	lwgeom1 = lwgeom_from_gserialized(geom1) ;
 	lwgeom2 = lwgeom_from_gserialized(geom2) ;
 
-	lwresult = lwgeom_symdifference(lwgeom1, lwgeom2) ;
+	lwresult = lwgeom_symdifference_prec(lwgeom1, lwgeom2, prec);
 	result = geometry_serialize(lwresult) ;
 
 	lwgeom_free(lwgeom1) ;
@@ -1406,14 +1416,17 @@ Datum ST_Intersection(PG_FUNCTION_ARGS)
 	GSERIALIZED *geom2;
 	GSERIALIZED *result;
 	LWGEOM *lwgeom1, *lwgeom2, *lwresult;
+	double prec = -1;
 
 	geom1 = PG_GETARG_GSERIALIZED_P(0);
 	geom2 = PG_GETARG_GSERIALIZED_P(1);
+	if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
+		prec = PG_GETARG_FLOAT8(2);
 
 	lwgeom1 = lwgeom_from_gserialized(geom1);
 	lwgeom2 = lwgeom_from_gserialized(geom2);
 
-	lwresult = lwgeom_intersection(lwgeom1, lwgeom2);
+	lwresult = lwgeom_intersection_prec(lwgeom1, lwgeom2, prec);
 	result = geometry_serialize(lwresult);
 
 	lwgeom_free(lwgeom1);
@@ -1433,14 +1446,17 @@ Datum ST_Difference(PG_FUNCTION_ARGS)
 	GSERIALIZED *geom2;
 	GSERIALIZED *result;
 	LWGEOM *lwgeom1, *lwgeom2, *lwresult;
+	double prec = -1;
 
 	geom1 = PG_GETARG_GSERIALIZED_P(0);
 	geom2 = PG_GETARG_GSERIALIZED_P(1);
+	if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
+		prec = PG_GETARG_FLOAT8(2);
 
 	lwgeom1 = lwgeom_from_gserialized(geom1);
 	lwgeom2 = lwgeom_from_gserialized(geom2);
 
-	lwresult = lwgeom_difference(lwgeom1, lwgeom2);
+	lwresult = lwgeom_difference_prec(lwgeom1, lwgeom2, prec);
 	result = geometry_serialize(lwresult);
 
 	lwgeom_free(lwgeom1);
diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
index 8fcd15f..9b19f05 100644
--- a/postgis/postgis.sql.in
+++ b/postgis/postgis.sql.in
@@ -3422,7 +3422,8 @@ CREATE OR REPLACE FUNCTION ST_IsValidTrajectory(geometry)
 -- GEOS
 ---------------------------------------------------------------
 
-CREATE OR REPLACE FUNCTION ST_Intersection(geom1 geometry, geom2 geometry)
+-- Changed: 3.1.0 to add gridSize default argument
+CREATE OR REPLACE FUNCTION ST_Intersection(geom1 geometry, geom2 geometry, gridSize float8 DEFAULT -1)
 	RETURNS geometry
 	AS 'MODULE_PATHNAME','ST_Intersection'
 	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
@@ -3563,8 +3564,9 @@ CREATE OR REPLACE FUNCTION ST_MaximumInscribedCircle(geometry, OUT center geomet
 	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
 	_COST_HIGH;
 
--- PostGIS equivalent function: difference(geom1 geometry, geom2 geometry)
-CREATE OR REPLACE FUNCTION ST_Difference(geom1 geometry, geom2 geometry)
+-- PostGIS equivalent function: ST_difference(geom1 geometry, geom2 geometry)
+-- Changed: 3.1.0 to add gridSize default argument
+CREATE OR REPLACE FUNCTION ST_Difference(geom1 geometry, geom2 geometry, gridSize float8 DEFAULT -1.0)
 	RETURNS geometry
 	AS 'MODULE_PATHNAME','ST_Difference'
 	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
@@ -3585,28 +3587,35 @@ CREATE OR REPLACE FUNCTION ST_Points(geometry)
 	_COST_MEDIUM;
 
 -- PostGIS equivalent function: symdifference(geom1 geometry, geom2 geometry)
-CREATE OR REPLACE FUNCTION ST_SymDifference(geom1 geometry, geom2 geometry)
+-- Changed: 3.1.0 to add gridSize default argument
+CREATE OR REPLACE FUNCTION ST_SymDifference(geom1 geometry, geom2 geometry, gridSize float8 DEFAULT -1.0)
 	RETURNS geometry
-	AS 'MODULE_PATHNAME','symdifference'
+	AS 'MODULE_PATHNAME','ST_SymDifference'
 	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
 	_COST_HIGH;
 
 -- Availability: 1.2.2
 CREATE OR REPLACE FUNCTION ST_SymmetricDifference(geom1 geometry, geom2 geometry)
 	RETURNS geometry
-	AS 'MODULE_PATHNAME','symdifference'
+	AS 'SELECT ST_SymDifference(geom1, geom2, -1.0);'
+	LANGUAGE 'sql';
+
+CREATE OR REPLACE FUNCTION ST_Union(geom1 geometry, geom2 geometry)
+	RETURNS geometry
+	AS 'MODULE_PATHNAME','ST_Union'
 	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
 	_COST_HIGH;
 
--- PostGIS equivalent function: GeomUnion(geom1 geometry, geom2 geometry)
-CREATE OR REPLACE FUNCTION ST_Union(geom1 geometry, geom2 geometry)
+-- Availability: 3.1.0
+CREATE OR REPLACE FUNCTION ST_Union(geom1 geometry, geom2 geometry, gridSize float8)
 	RETURNS geometry
 	AS 'MODULE_PATHNAME','ST_Union'
 	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
 	_COST_HIGH;
 
 -- Availability: 2.0.0
-CREATE OR REPLACE FUNCTION ST_UnaryUnion(geometry)
+-- Changed: 3.1.0 to add gridSize default argument
+CREATE OR REPLACE FUNCTION ST_UnaryUnion(geometry, gridSize float8 DEFAULT -1.0)
 	RETURNS geometry
 	AS 'MODULE_PATHNAME','ST_UnaryUnion'
 	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
@@ -3633,7 +3642,8 @@ CREATE OR REPLACE FUNCTION ST_ClipByBox2d(geom geometry, box box2d)
 	_COST_MEDIUM;
 
 -- Availability: 2.2.0
-CREATE OR REPLACE FUNCTION ST_Subdivide(geom geometry, maxvertices integer DEFAULT 256)
+-- Changed: 3.1.0 to add gridSize default argument
+CREATE OR REPLACE FUNCTION ST_Subdivide(geom geometry, maxvertices integer DEFAULT 256, gridSize float8 DEFAULT -1.0)
 	RETURNS setof geometry
 	AS 'MODULE_PATHNAME', 'ST_Subdivide'
 	LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
@@ -4006,6 +4016,15 @@ CREATE AGGREGATE ST_Union (geometry) (
 	finalfunc = pgis_geometry_union_finalfn
 	);
 
+-- Availability: 3.1.0
+CREATE AGGREGATE ST_Union (geometry, gridSize float8) (
+	sfunc = pgis_geometry_accum_transfn,
+	stype = internal,
+	parallel = safe,
+	finalfunc = pgis_geometry_union_finalfn
+	);
+
+
 -- Availability: 1.2.2
 -- Changed: 2.4.0: marked parallel safe
 -- Changed: 2.5.0 use 'internal' stype
diff --git a/postgis/postgis_after_upgrade.sql b/postgis/postgis_after_upgrade.sql
index 02f8692..a7b929f 100644
--- a/postgis/postgis_after_upgrade.sql
+++ b/postgis/postgis_after_upgrade.sql
@@ -4,7 +4,7 @@
 -- PostGIS - Spatial Types for PostgreSQL
 -- http://postgis.net
 --
--- Copyright (C) 2011-2012 Sandro Santilli <strk at kbt.io>
+-- Copyright (C) 2011-2020 Sandro Santilli <strk at kbt.io>
 -- Copyright (C) 2010-2012 Regina Obe <lr at pcorp.us>
 -- Copyright (C) 2009      Paul Ramsey <pramsey at cleverelephant.ca>
 --
diff --git a/postgis/postgis_before_upgrade.sql b/postgis/postgis_before_upgrade.sql
index fd104bc..dcd69aa 100644
--- a/postgis/postgis_before_upgrade.sql
+++ b/postgis/postgis_before_upgrade.sql
@@ -217,6 +217,12 @@ DROP FUNCTION IF EXISTS st_buffer(geometry, double precision); -- Does not confl
 DROP FUNCTION IF EXISTS ST_CurveToLine(geometry, integer); -- Does not conflict
 DROP FUNCTION IF EXISTS ST_CurveToLine(geometry); -- Does not conflict
 
+DROP FUNCTION IF EXISTS st_intersection(geometry, geometry); -- replaced in 3.1.0 by 3 args version
+DROP FUNCTION IF EXISTS st_subdivide(geometry, integer); -- replaced in 3.1.0 by 3 args version
+DROP FUNCTION IF EXISTS st_difference(geometry, geometry); -- replaced in 3.1.0 by 3 args version
+DROP FUNCTION IF EXISTS st_symdifference(geometry, geometry); -- replaced in 3.1.0 by 3 args version
+DROP FUNCTION IF EXISTS st_unaryunion(geometry); -- replaced in 3.1.0 by 3 args version
+
 -- geometry_columns changed parameter types so we verify if it needs to be dropped
 -- We check the catalog to see if the view (geometry_columns) has a column
 -- with name `f_table_schema` and type `character varying(256)` as it was
diff --git a/regress/core/subdivide_expected b/regress/core/subdivide_expected
index bf58e60..10fc0e5 100644
--- a/regress/core/subdivide_expected
+++ b/regress/core/subdivide_expected
@@ -1,7 +1,7 @@
 1|t|5|10
 2|t|4|8
 3|t|10|10
-ERROR:  lwgeom_subdivide: cannot subdivide to fewer than 5 vertices per output
+ERROR:  lwgeom_subdivide_prec: cannot subdivide to fewer than 5 vertices per output
 #3522|POINT(1 1)
 #3744|1600000000000000
 4|29321996468.6|29321996468.6|1857|256
diff --git a/regress/core/tickets_expected b/regress/core/tickets_expected
index a329953..9bbed01 100644
--- a/regress/core/tickets_expected
+++ b/regress/core/tickets_expected
@@ -327,10 +327,10 @@ NOTICE:  Too few points in geometry component at or near point 0 0
 #4011|ST_MultiLineString|MULTILINESTRING EMPTY|t|t
 #4011|ST_GeometryCollection|MULTILINESTRING((0 0,0 0))|f|f
 #4025|
-ERROR:  lwgeom_intersection: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection
-ERROR:  lwgeom_difference: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection
-ERROR:  lwgeom_symdifference: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection
-ERROR:  lwgeom_union: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection
+ERROR:  lwgeom_intersection_prec: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection
+ERROR:  lwgeom_difference_prec: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection
+ERROR:  lwgeom_symdifference_prec: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection
+ERROR:  lwgeom_union_prec: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection
 #4055a|4326
 #4055b|4326
 #4089|LINESTRING Z (1 1 1,3 3 1)

-----------------------------------------------------------------------

Summary of changes:
 doc/reference_processing.xml       |  54 ++++++++++++++++++-
 liblwgeom/liblwgeom.h.in           |   6 +++
 liblwgeom/lwgeom.c                 |  28 ++++++----
 liblwgeom/lwgeom_geos.c            | 106 ++++++++++++++++++++++++++++++++++---
 postgis/legacy.sql.in              |   4 +-
 postgis/lwgeom_accum.c             |   9 ++++
 postgis/lwgeom_accum.h             |   1 +
 postgis/lwgeom_dump.c              |   9 +++-
 postgis/lwgeom_geos.c              |  64 +++++++++++++---------
 postgis/postgis.sql.in             |  39 ++++++++++----
 postgis/postgis_after_upgrade.sql  |   2 +-
 postgis/postgis_before_upgrade.sql |   6 +++
 regress/core/subdivide_expected    |   2 +-
 regress/core/tickets_expected      |   8 +--
 14 files changed, 277 insertions(+), 61 deletions(-)


hooks/post-receive
-- 
PostGIS


More information about the postgis-tickets mailing list