[postgis-tickets] [SCM] PostGIS branch master updated. 3.4.0beta1-81-gfab44d035
git at osgeo.org
git at osgeo.org
Fri Jul 28 12:27:35 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 fab44d035a2445ffbd916536ba2ee884cd336ffd (commit)
from 09a0c95620008fa32f8ca8cca246ed3071161ce2 (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 fab44d035a2445ffbd916536ba2ee884cd336ffd
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Fri Jul 28 12:27:26 2023 -0700
Reflect updates from MobilityDB on use of spheroid and line-locate-point, closes #5460
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index a27529e6c..69ed8f354 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -1793,6 +1793,7 @@ extern int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwge
extern LWGEOM * geography_substring(const LWLINE *line,
+ const SPHEROID *s,
double from, double to,
double tolerance);
diff --git a/liblwgeom/lwgeodetic_measures.c b/liblwgeom/lwgeodetic_measures.c
index e1fad3a7f..1bca0a795 100644
--- a/liblwgeom/lwgeodetic_measures.c
+++ b/liblwgeom/lwgeodetic_measures.c
@@ -55,72 +55,55 @@
* Find interpolation point p between geography points p1 and p2
* so that the len(p1,p) == len(p1,p2)
* f and p falls on p1,p2 segment
- *
- * @param[in] p1,p2 3D-space points we are interpolating between
- * @param[in] v1,v2 real values and z/m coordinates
- * @param[in] f Fraction
- * @param[out] p Result
*/
static void
-interpolate_point4d_sphere(
- const POINT3D *p1, const POINT3D *p2,
- const POINT4D *v1, const POINT4D *v2,
- double f, POINT4D *p)
-{
- /* Calculate interpolated point */
- POINT3D mid;
- mid.x = p1->x + ((p2->x - p1->x) * f);
- mid.y = p1->y + ((p2->y - p1->y) * f);
- mid.z = p1->z + ((p2->z - p1->z) * f);
- normalize(&mid);
-
- /* Calculate z/m values */
- GEOGRAPHIC_POINT g;
- cart2geog(&mid, &g);
- p->x = rad2deg(g.lon);
- p->y = rad2deg(g.lat);
- p->z = v1->z + ((v2->z - v1->z) * f);
- p->m = v1->m + ((v2->m - v1->m) * f);
-}
-
+interpolate_point4d_spheroid(
+ const POINT4D *p1, const POINT4D *p2, POINT4D *p,
+ const SPHEROID *s, double f)
-
-/**
- * @brief Returns the length of the point array wrt the sphere
- */
-static double
-ptarray_length_sphere(const POINTARRAY *pa)
{
- GEOGRAPHIC_POINT a, b;
- POINT4D p;
- uint32_t i;
- double length = 0.0;
+ GEOGRAPHIC_POINT g, g1, g2;
+ geographic_point_init(p1->x, p1->y, &g1);
+ geographic_point_init(p2->x, p2->y, &g2);
+ int success;
+ double dist, dir;
- /* Return zero on non-sensical inputs */
- if ( ! pa || pa->npoints < 2 )
- return 0.0;
-
- /* Initialize first point */
- getPoint4d_p(pa, 0, &p);
- geographic_point_init(p.x, p.y, &a);
+ /* Special sphere case */
+ if ( s == NULL || s->a == s->b )
+ {
+ /* Calculate distance and direction between g1 and g2 */
+ dist = sphere_distance(&g1, &g2);
+ dir = sphere_direction(&g1, &g2, dist);
+ /* Compute interpolation point */
+ success = sphere_project(&g1, dist*f, dir, &g);
+ }
+ /* Spheroid case */
+ else
+ {
+ /* Calculate distance and direction between g1 and g2 */
+ dist = spheroid_distance(&g1, &g2, s);
+ dir = spheroid_direction(&g1, &g2, s);
+ /* Compute interpolation point */
+ success = spheroid_project(&g1, s, dist*f, dir, &g);
+ }
- /* Loop and sum the length for each segment */
- for ( i = 1; i < pa->npoints; i++ )
+ /* If success, use newly computed lat and lon,
+ * otherwise return precomputed cartesian result */
+ if (success == LW_SUCCESS)
{
- getPoint4d_p(pa, i, &p);
- geographic_point_init(p.x, p.y, &b);
- /* Add this segment length to the total */
- length += sphere_distance(&a, &b);
+ p->x = rad2deg(longitude_radians_normalize(g.lon));
+ p->y = rad2deg(latitude_radians_normalize(g.lat));
}
- return length;
}
+
/**
* @brief Return the part of a line between two fractional locations.
*/
LWGEOM *
geography_substring(
const LWLINE *lwline,
+ const SPHEROID *s,
double from, double to,
double tolerance)
{
@@ -129,7 +112,6 @@ geography_substring(
LWGEOM *lwresult;
POINT4D pt;
POINT4D p1, p2;
- POINT3D q1, q2;
GEOGRAPHIC_POINT g1, g2;
int nsegs, i;
double length, slength, tlength;
@@ -146,7 +128,7 @@ geography_substring(
ipa->npoints);
/* Compute total line length */
- length = ptarray_length_sphere(ipa);
+ length = ptarray_length_spheroid(ipa, s);
/* Get 'from' and 'to' lengths */
from = length * from;
@@ -162,7 +144,12 @@ geography_substring(
geographic_point_init(p2.x, p2.y, &g2);
/* Find the length of this segment */
- slength = sphere_distance(&g1, &g2);
+ /* Special sphere case */
+ if ( s->a == s->b )
+ slength = s->radius * sphere_distance(&g1, &g2);
+ /* Spheroid case */
+ else
+ slength = spheroid_distance(&g1, &g2, s);
/* We are before requested start. */
if (state == 0) /* before */
@@ -196,9 +183,7 @@ geography_substring(
{
/* Our start is between first and second point */
dseg = (from - tlength) / slength;
- geog2cart(&g1, &q1);
- geog2cart(&g2, &q2);
- interpolate_point4d_sphere(&q1, &q2, &p1, &p2, dseg, &pt);
+ interpolate_point4d_spheroid(&p1, &p2, &pt, s, dseg);
ptarray_append_point(dpa, &pt, LW_FALSE);
/* We're inside now, but will check 'to' point as well */
state = 1;
@@ -238,9 +223,7 @@ geography_substring(
else if (to < tlength + slength )
{
dseg = (to - tlength) / slength;
- geog2cart(&g1, &q1);
- geog2cart(&g2, &q2);
- interpolate_point4d_sphere(&q1, &q2, &p1, &p2, dseg, &pt);
+ interpolate_point4d_spheroid(&p1, &p2, &pt, s, dseg);
ptarray_append_point(dpa, &pt, LW_FALSE);
break;
}
@@ -322,9 +305,16 @@ geography_interpolate_points(
geographic_point_init(p1.x, p1.y, &g1);
for ( i = 0; i < ipa->npoints - 1 && points_found < points_to_interpolate; i++ )
{
+ double segment_length_frac;
getPoint4d_p(ipa, i+1, &p2);
geographic_point_init(p2.x, p2.y, &g2);
- double segment_length_frac = spheroid_distance(&g1, &g2, s) / length;
+
+ /* Special sphere case */
+ if ( s->a == s->b )
+ segment_length_frac = s->radius * sphere_distance(&g1, &g2) / length;
+ /* Spheroid case */
+ else
+ segment_length_frac = spheroid_distance(&g1, &g2, s) / length;
/* If our target distance is before the total length we've seen
* so far. create a new point some distance down the current
@@ -335,7 +325,7 @@ geography_interpolate_points(
geog2cart(&g1, &q1);
geog2cart(&g2, &q2);
double segment_fraction = (length_fraction - length_fraction_consumed) / segment_length_frac;
- interpolate_point4d_sphere(&q1, &q2, &p1, &p2, segment_fraction, &pt);
+ interpolate_point4d_spheroid(&p1, &p2, &pt, s, segment_fraction);
ptarray_set_point4d(opa, points_found++, &pt);
length_fraction += length_fraction_increment;
}
diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index e52671f4f..dc8bde060 100644
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -1218,7 +1218,12 @@ Datum geography_line_substring(PG_FUNCTION_ARGS)
double to_fraction = PG_GETARG_FLOAT8(2);
LWLINE *lwline;
LWGEOM *lwresult;
+ SPHEROID s;
GSERIALIZED *result;
+ bool use_spheroid = true;
+
+ if ( PG_NARGS() > 3 && ! PG_ARGISNULL(3) )
+ use_spheroid = PG_GETARG_BOOL(3);
/* Return NULL on empty argument. */
if ( gserialized_is_empty(gs) )
@@ -1253,7 +1258,13 @@ Datum geography_line_substring(PG_FUNCTION_ARGS)
PG_RETURN_NULL();
}
- lwresult = geography_substring(lwline,
+ /* Initialize spheroid */
+ spheroid_init_from_srid(gserialized_get_srid(gs), &s);
+ /* Set to sphere if requested */
+ if ( ! use_spheroid )
+ s.a = s.b = s.radius;
+
+ lwresult = geography_substring(lwline, &s,
from_fraction, to_fraction, FP_TOLERANCE);
lwline_free(lwline);
@@ -1311,8 +1322,7 @@ Datum geography_line_interpolate_point(PG_FUNCTION_ARGS)
}
/* Initialize spheroid */
- // spheroid_init_from_srid(fcinfo, srid, &s);
- spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
+ spheroid_init_from_srid(gserialized_get_srid(gs), &s);
/* Set to sphere if requested */
if ( ! use_spheroid )
@@ -1376,8 +1386,7 @@ Datum geography_line_locate_point(PG_FUNCTION_ARGS)
}
else {
/* Initialize spheroid */
- // spheroid_init_from_srid(fcinfo, srid, &s);
- spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
+ spheroid_init_from_srid(gserialized_get_srid(gs1), &s);
}
lwline = lwgeom_as_lwline(lwgeom_from_gserialized(gs1));
@@ -1457,8 +1466,7 @@ Datum geography_shortestline(PG_FUNCTION_ARGS)
}
/* Initialize spheroid */
- // spheroid_init_from_srid(fcinfo, srid, &s);
- spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
+ spheroid_init_from_srid(gserialized_get_srid(g1), &s);
/* Set to sphere if requested */
if ( ! use_spheroid )
diff --git a/regress/core/geography.sql b/regress/core/geography.sql
index 1ad7e4e04..1c07be981 100644
--- a/regress/core/geography.sql
+++ b/regress/core/geography.sql
@@ -139,16 +139,17 @@ SELECT 'lrs_lip_1', ST_AsText(ST_LineInterpolatePoint(geography 'Linestring(4.35
SELECT 'lrs_lip_2', ST_AsText(ST_LineInterpolatePoints(geography 'Linestring(4.35 50.85, 37.617222 55.755833)', 0.0, true), 2);
SELECT 'lrs_lip_3', ST_AsText(ST_LineInterpolatePoints(geography 'Linestring(4.35 50.85, 37.617222 55.755833)', 1.0, false), 2);
SELECT 'lrs_lip_4', ST_AsText(ST_LineInterpolatePoints(geography 'Linestring(4.35 50.85, 37.617222 55.755833)', 0.1, true), 2);
+SELECT 'lrs_lip_5', ST_AsText(ST_LineInterpolatePoints(geography 'Linestring(4.35 50.85, 37.617222 55.755833)', 0.1, false), 2);
SELECT 'lrs_llp_1', round(ST_LineLocatePoint(geography 'linestring(0 1, 50 1)', geography 'Point(25 0)')::numeric, 2);
SELECT 'lrs_llp_2', round(ST_LineLocatePoint(geography 'linestring(0 1, 50 1)', geography 'Point(-5 0)')::numeric, 2);
SELECT 'lrs_llp_3', round(ST_LineLocatePoint(geography 'linestring(0 1, 50 1)', geography 'Point(55 0)')::numeric, 2);
SELECT 'lrs_llp_4', round(ST_LineLocatePoint(geography 'linestring(0 1, 50 1)', geography 'Linestring(25 0,26 1)')::numeric, 2);
-SELECT 'lrs_substr_1', ST_AsText(ST_LineSubstring(geography 'linestring(0 1, 100 1)', 0.1, 0.2),2);
+SELECT 'lrs_substr_1', ST_AsText(ST_LineSubstring(geography 'linestring(0 20, 100 20)', 0.1, 0.2),2);
---SELECT 'lrs_cp_1', ST_AsText(ST_ClosestPoint(geography 'linestring(0 1, 50 1)', 'Point(25 0)'), 3);
---SELECT 'lrs_cp_2', ST_AsText(ST_ClosestPoint(geography 'Point(25 0)', geography 'linestring(0 1, 50 1)'), 3);
+--SELECT 'lrs_cp_1', ST_AsText(ST_ClosestPoint(geography 'Linestring(0 20, 50 20)', 'Point(25 20)'), 3);
+--SELECT 'lrs_cp_2', ST_AsText(ST_ClosestPoint(geography 'Point(25 20)', geography 'Linestring(0 20, 50 20)'), 3);
-SELECT 'lrs_sl_1', ST_AsText(ST_ShortestLine(geography 'linestring(0 1, 50 1)', 'Point(25 0)', true), 2);
+SELECT 'lrs_sl_1', ST_AsText(ST_ShortestLine(geography 'linestring(0 40, 50 40)', 'Point(25 40)', true), 2);
diff --git a/regress/core/geography_expected b/regress/core/geography_expected
index 4fc2e88c5..ae740cd52 100644
--- a/regress/core/geography_expected
+++ b/regress/core/geography_expected
@@ -58,10 +58,11 @@ lrs_empty_6|
lrs_lip_1|POINT(4.35 50.85)
lrs_lip_2|POINT(4.35 50.85)
lrs_lip_3|POINT(37.62 55.76)
-lrs_lip_4|MULTIPOINT((7.22 51.72),(10.23 52.53),(13.37 53.26),(16.63 53.91),(20 54.46),(23.45 54.93),(26.96 55.29),(30.51 55.55),(34.07 55.7),(37.62 55.76))
+lrs_lip_4|MULTIPOINT((7.27 51.73),(10.3 52.54),(13.43 53.27),(16.67 53.91),(20 54.47),(23.42 54.93),(26.9 55.29),(30.44 55.55),(34.02 55.7),(37.62 55.76))
+lrs_lip_5|MULTIPOINT((7.27 51.73),(10.29 52.54),(13.43 53.27),(16.67 53.91),(20 54.46),(23.42 54.92),(26.9 55.29),(30.44 55.55),(34.02 55.7),(37.62 55.76))
lrs_llp_1|0.50
lrs_llp_2|0.00
lrs_llp_3|1.00
ERROR: geography_line_locate_point: 2st arg is not a point
-lrs_substr_1|LINESTRING(6.37 1.13,14.43 1.27)
-lrs_sl_1|LINESTRING(25 0,25 1.1)
+lrs_substr_1|LINESTRING(9.28 23.25,18.97 25.93)
+lrs_sl_1|LINESTRING(25 40,25 42.79)
-----------------------------------------------------------------------
Summary of changes:
liblwgeom/liblwgeom.h.in | 1 +
liblwgeom/lwgeodetic_measures.c | 112 ++++++++++++++++++----------------------
postgis/geography_measurement.c | 22 +++++---
regress/core/geography.sql | 9 ++--
regress/core/geography_expected | 7 +--
5 files changed, 76 insertions(+), 75 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list