[SCM] PostGIS branch master updated. 3.6.0rc2-509-gbbbe84257
git at osgeo.org
git at osgeo.org
Tue Jun 9 12:40:11 PDT 2026
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, master has been updated
via bbbe84257be29b0101088e5da69677924f65cdd0 (commit)
via 8da97650d154e7b70b0a58f5be4c2870848bd06c (commit)
from 6216b684aaff36b393d5d5707c95a896b1e89ffc (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 bbbe84257be29b0101088e5da69677924f65cdd0
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Tue Jun 9 12:40:05 2026 -0700
News item for #4560
diff --git a/NEWS b/NEWS
index aaf9f2c05..17b07b97f 100644
--- a/NEWS
+++ b/NEWS
@@ -40,6 +40,7 @@ To take advantage of all postgis_sfcgal extension features SFCGAL 2.2+ is needed
- #2614, Use GEOSPreparedDistance and caching to accelerate ST_DWithin (Paul Ramsey)
- GH-839, ST_Multi support for TIN and surfaces (Loïc Bartoletti)
- #4398, Add ST_CatmullSmoothing smoothes but retains original vertices (Paul Ramsey)
+ - #4560, ST_3DInterpolatePoint for M interpolation from XYZ inputs (Paul Ramsey)
* Enhancements *
commit 8da97650d154e7b70b0a58f5be4c2870848bd06c
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Tue Jun 9 19:32:57 2026 +0000
Add ST_3DInterpolatePoint for ZM geometry measure interpolation
Projects a query point onto a LINESTRING ZM using 3D (XYZ) distance,
returning the interpolated M value. Closes #4560.
diff --git a/doc/reference_lrs.xml b/doc/reference_lrs.xml
index 570d4daea..c0c797c65 100644
--- a/doc/reference_lrs.xml
+++ b/doc/reference_lrs.xml
@@ -744,6 +744,68 @@ SELECT ST_AsText(
<para><xref linkend="ST_AddMeasure"/>, <xref linkend="ST_LocateAlong"/>, <xref linkend="ST_LocateBetween"/></para>
</refsection>
+ </refentry>
+
+ <refentry xml:id="ST_3DInterpolatePoint">
+ <refnamediv>
+ <refname>ST_3DInterpolatePoint</refname>
+
+ <refpurpose>Returns the interpolated measure of a geometry closest to a point in 3D.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>float8 <function>ST_3DInterpolatePoint</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>linear_geom_with_measure</parameter></paramdef>
+ <paramdef><type>geometry </type> <parameter>point</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>Returns the interpolated measure value of a linear ZM geometry
+ at the location closest to the given point, using 3D (XYZ) distance
+ for the projection. Use this function when the geometry has
+ significant Z variation, such as flight trajectories, where
+ <xref linkend="ST_InterpolatePoint"/> would give incorrect results
+ by ignoring the Z dimension.</para>
+
+ <note><para>The line must have both Z and M dimensions.
+ The point should have a Z dimension.</para></note>
+
+ <para role="availability" conformance="3.7.0">Availability: 3.7.0</para>
+
+ <para>&Z_support;</para>
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+
+ <programlisting>-- Line rising diagonally in Z; point at the 3D midpoint
+SELECT ST_3DInterpolatePoint(
+ 'LINESTRING ZM (0 0 0 0, 10 0 10 20)',
+ 'POINT Z (5 0 5)');
+ ---------------------
+ 10
+
+-- Vertical line (zero XY extent); 3D projection works where 2D would not
+SELECT ST_3DInterpolatePoint(
+ 'LINESTRING ZM (0 0 0 0, 0 0 10 100)',
+ 'POINT Z (5 5 5)');
+ ---------------------
+ 50
+ </programlisting>
+ </refsection>
+
+ <refsection>
+ <title>See Also</title>
+ <para><xref linkend="ST_InterpolatePoint"/>, <xref linkend="ST_AddMeasure"/>,
+ <xref linkend="ST_LocateAlong"/>, <xref linkend="ST_LocateBetween"/></para>
+ </refsection>
+
</refentry>
diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c
index 1194785c4..6ac9da3c7 100644
--- a/liblwgeom/cunit/cu_algorithm.c
+++ b/liblwgeom/cunit/cu_algorithm.c
@@ -596,6 +596,52 @@ test_lwline_interpolate_point_3d(void)
lwline_free(line);
}
+static void
+test_lwgeom_interpolate_point_3d(void)
+{
+ LWGEOM *line;
+ LWPOINT *pt;
+ double result;
+
+ /* Simple diagonal line in XYZ, point at 3D midpoint */
+ line = lwgeom_from_wkt("LINESTRING ZM (0 0 0 0, 10 0 10 20)", LW_PARSER_CHECK_NONE);
+ pt = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT Z (5 0 5)", LW_PARSER_CHECK_NONE));
+ result = lwgeom_interpolate_point_3d(line, pt);
+ CU_ASSERT_DOUBLE_EQUAL(result, 10.0, 1e-10);
+ lwpoint_free(pt);
+
+ /* Point offset in Y but same XZ projection: same result */
+ pt = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT Z (5 99 5)", LW_PARSER_CHECK_NONE));
+ result = lwgeom_interpolate_point_3d(line, pt);
+ CU_ASSERT_DOUBLE_EQUAL(result, 10.0, 1e-10);
+ lwpoint_free(pt);
+ lwgeom_free(line);
+
+ /* Vertical line (zero XY extent): 3D projection works where 2D would not */
+ line = lwgeom_from_wkt("LINESTRING ZM (0 0 0 0, 0 0 10 100)", LW_PARSER_CHECK_NONE);
+ pt = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT Z (5 5 5)", LW_PARSER_CHECK_NONE));
+ result = lwgeom_interpolate_point_3d(line, pt);
+ CU_ASSERT_DOUBLE_EQUAL(result, 50.0, 1e-10);
+ lwpoint_free(pt);
+ lwgeom_free(line);
+
+ /* Point beyond end of line: clamped to end M */
+ line = lwgeom_from_wkt("LINESTRING ZM (0 0 0 0, 10 0 0 100)", LW_PARSER_CHECK_NONE);
+ pt = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT Z (20 0 0)", LW_PARSER_CHECK_NONE));
+ result = lwgeom_interpolate_point_3d(line, pt);
+ CU_ASSERT_DOUBLE_EQUAL(result, 100.0, 1e-10);
+ lwpoint_free(pt);
+ lwgeom_free(line);
+
+ /* Point before start of line: clamped to start M */
+ line = lwgeom_from_wkt("LINESTRING ZM (0 0 0 0, 10 0 0 100)", LW_PARSER_CHECK_NONE);
+ pt = lwgeom_as_lwpoint(lwgeom_from_wkt("POINT Z (-5 0 0)", LW_PARSER_CHECK_NONE));
+ result = lwgeom_interpolate_point_3d(line, pt);
+ CU_ASSERT_DOUBLE_EQUAL(result, 0.0, 1e-10);
+ lwpoint_free(pt);
+ lwgeom_free(line);
+}
+
static void test_lwline_clip(void)
{
LWCOLLECTION *c;
@@ -1844,6 +1890,7 @@ void algorithms_suite_setup(void)
PG_ADD_TEST(suite,test_point_interpolate);
PG_ADD_TEST(suite,test_lwline_interpolate_points);
PG_ADD_TEST(suite, test_lwline_interpolate_point_3d);
+ PG_ADD_TEST(suite, test_lwgeom_interpolate_point_3d);
PG_ADD_TEST(suite,test_lwline_clip);
PG_ADD_TEST(suite, test_lwpoly_clip);
PG_ADD_TEST(suite, test_lwtriangle_clip);
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index 3ac202217..692167c13 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -1623,6 +1623,7 @@ extern LWCOLLECTION* lwgeom_locate_between(const LWGEOM *lwin, double from, doub
* Find the measure value at the location on the line closest to the point.
*/
extern double lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt);
+extern double lwgeom_interpolate_point_3d(const LWGEOM *lwin, const LWPOINT *lwpt);
/**
* Find the time of closest point of approach
diff --git a/liblwgeom/lwlinearreferencing.c b/liblwgeom/lwlinearreferencing.c
index 610f0cf4e..ace84cc08 100644
--- a/liblwgeom/lwlinearreferencing.c
+++ b/liblwgeom/lwlinearreferencing.c
@@ -947,6 +947,86 @@ lwgeom_interpolate_point(const LWGEOM *lwin, const LWPOINT *lwpt)
return ret;
}
+double
+lwgeom_interpolate_point_3d(const LWGEOM *lwin, const LWPOINT *lwpt)
+{
+ POINT4D p, A, B, closest;
+ double mindist2 = DBL_MAX;
+ double best_m = 0.0;
+ uint32_t i;
+
+ if (!lwin)
+ lwerror("lwgeom_interpolate_point_3d: null input geometry!");
+
+ if (!lwgeom_has_m(lwin))
+ lwerror("Input geometry does not have a measure dimension");
+
+ if (!lwgeom_has_z(lwin))
+ lwerror("Input geometry does not have a Z dimension");
+
+ if (lwgeom_is_empty(lwin) || lwpoint_is_empty(lwpt))
+ lwerror("Input geometry is empty");
+
+ switch (lwin->type)
+ {
+ case LINETYPE:
+ {
+ LWLINE *lwline = lwgeom_as_lwline(lwin);
+ POINTARRAY *pa = lwline->points;
+
+ lwpoint_getPoint4d_p(lwpt, &p);
+
+ if (pa->npoints == 1)
+ {
+ getPoint4d_p(pa, 0, &closest);
+ return closest.m;
+ }
+
+ for (i = 0; i < pa->npoints - 1; i++)
+ {
+ double dx, dy, dz, len2, r;
+ double cx, cy, cz, d2;
+
+ getPoint4d_p(pa, i, &A);
+ getPoint4d_p(pa, i + 1, &B);
+
+ dx = B.x - A.x;
+ dy = B.y - A.y;
+ dz = B.z - A.z;
+ len2 = dx * dx + dy * dy + dz * dz;
+
+ if (len2 == 0.0)
+ {
+ /* Degenerate segment: distance to point A */
+ r = 0.0;
+ }
+ else
+ {
+ r = ((p.x - A.x) * dx + (p.y - A.y) * dy + (p.z - A.z) * dz) / len2;
+ if (r < 0.0) r = 0.0;
+ if (r > 1.0) r = 1.0;
+ }
+
+ cx = A.x + r * dx;
+ cy = A.y + r * dy;
+ cz = A.z + r * dz;
+
+ d2 = (p.x - cx) * (p.x - cx) + (p.y - cy) * (p.y - cy) + (p.z - cz) * (p.z - cz);
+
+ if (d2 < mindist2)
+ {
+ mindist2 = d2;
+ best_m = A.m + r * (B.m - A.m);
+ }
+ }
+ break;
+ }
+ default:
+ lwerror("This function does not accept %s geometries.", lwtype_name(lwin->type));
+ }
+ return best_m;
+}
+
/*
* Time of closest point of approach
*
diff --git a/postgis/lwgeom_functions_lrs.c b/postgis/lwgeom_functions_lrs.c
index 474d07ea9..b02afb5bb 100644
--- a/postgis/lwgeom_functions_lrs.c
+++ b/postgis/lwgeom_functions_lrs.c
@@ -214,6 +214,47 @@ Datum ST_InterpolatePoint(PG_FUNCTION_ARGS)
}
+Datum ST_3DInterpolatePoint(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(ST_3DInterpolatePoint);
+Datum ST_3DInterpolatePoint(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *gser_line = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED *gser_point = PG_GETARG_GSERIALIZED_P(1);
+ LWGEOM *lwline;
+ LWPOINT *lwpoint;
+
+ if ( gserialized_get_type(gser_line) != LINETYPE )
+ {
+ elog(ERROR,"ST_3DInterpolatePoint: 1st argument isn't a line");
+ PG_RETURN_NULL();
+ }
+ if ( gserialized_get_type(gser_point) != POINTTYPE )
+ {
+ elog(ERROR,"ST_3DInterpolatePoint: 2nd argument isn't a point");
+ PG_RETURN_NULL();
+ }
+
+ gserialized_error_if_srid_mismatch(gser_line, gser_point, __func__);
+
+ if ( ! gserialized_has_m(gser_line) )
+ {
+ elog(ERROR,"ST_3DInterpolatePoint only accepts geometries that have an M dimension");
+ PG_RETURN_NULL();
+ }
+
+ if ( ! gserialized_has_z(gser_line) )
+ {
+ elog(ERROR,"ST_3DInterpolatePoint only accepts geometries that have a Z dimension");
+ PG_RETURN_NULL();
+ }
+
+ lwpoint = lwgeom_as_lwpoint(lwgeom_from_gserialized(gser_point));
+ lwline = lwgeom_from_gserialized(gser_line);
+
+ PG_RETURN_FLOAT8(lwgeom_interpolate_point_3d(lwline, lwpoint));
+}
+
+
Datum LWGEOM_line_locate_point(PG_FUNCTION_ARGS);
PG_FUNCTION_INFO_V1(LWGEOM_line_locate_point);
Datum LWGEOM_line_locate_point(PG_FUNCTION_ARGS)
diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
index 73d41a465..bf3a094d9 100644
--- a/postgis/postgis.sql.in
+++ b/postgis/postgis.sql.in
@@ -6682,6 +6682,13 @@ CREATE OR REPLACE FUNCTION ST_InterpolatePoint(Line geometry, Point geometry)
LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
_COST_MEDIUM;
+-- Availability: 3.7.0
+CREATE OR REPLACE FUNCTION ST_3DInterpolatePoint(Line geometry, Point geometry)
+ RETURNS float8
+ AS 'MODULE_PATHNAME', 'ST_3DInterpolatePoint'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
+ _COST_MEDIUM;
+
---------------------------------------------------------------
-- Grid / Hexagon coverage functions
---
diff --git a/regress/core/regress_lrs.sql b/regress/core/regress_lrs.sql
index 06e419dae..92706eb26 100644
--- a/regress/core/regress_lrs.sql
+++ b/regress/core/regress_lrs.sql
@@ -99,3 +99,11 @@ select 'addMeasure2', ST_AsText(ST_AddMeasure('LINESTRING(0 0, 9 0, 10 0)', 10,
select 'interpolatePoint1', ST_InterpolatePoint('LINESTRINGM(0 0 0, 10 0 4)', 'POINT(5 0)');
select 'interpolatePoint2', ST_InterpolatePoint('LINESTRINGM(0 0 0, 10 0 4)', 'POINT(7.5 0)');
+
+--
+-- ST_3DInterpolatePoint
+--
+
+select '3DinterpolatePoint1', ST_3DInterpolatePoint('LINESTRING ZM (0 0 0 0, 10 0 10 20)', 'POINT Z (5 0 5)');
+select '3DinterpolatePoint2', ST_3DInterpolatePoint('LINESTRING ZM (0 0 0 0, 10 0 10 20)', 'POINT Z (5 5 5)');
+select '3DinterpolatePoint3', ST_3DInterpolatePoint('LINESTRING ZM (0 0 0 0, 0 0 10 100)', 'POINT Z (5 5 5)');
diff --git a/regress/core/regress_lrs_expected b/regress/core/regress_lrs_expected
index bdc232637..868662a6f 100644
--- a/regress/core/regress_lrs_expected
+++ b/regress/core/regress_lrs_expected
@@ -47,3 +47,6 @@ addMeasure1|LINESTRING M (0 0 10,2 0 15,4 0 20)
addMeasure2|LINESTRING M (0 0 10,9 0 19,10 0 20)
interpolatePoint1|2
interpolatePoint2|3
+3DinterpolatePoint1|10
+3DinterpolatePoint2|10
+3DinterpolatePoint3|50
-----------------------------------------------------------------------
Summary of changes:
NEWS | 1 +
doc/reference_lrs.xml | 62 ++++++++++++++++++++++++++++++
liblwgeom/cunit/cu_algorithm.c | 47 +++++++++++++++++++++++
liblwgeom/liblwgeom.h.in | 1 +
liblwgeom/lwlinearreferencing.c | 80 +++++++++++++++++++++++++++++++++++++++
postgis/lwgeom_functions_lrs.c | 41 ++++++++++++++++++++
postgis/postgis.sql.in | 7 ++++
regress/core/regress_lrs.sql | 8 ++++
regress/core/regress_lrs_expected | 3 ++
9 files changed, 250 insertions(+)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list