mirror of
https://git.osgeo.org/gitea/postgis/postgis
synced 2024-10-23 00:22:38 +00:00
Reflect updates from MobilityDB on use of spheroid and line-locate-point, closes #5460
This commit is contained in:
Paul Ramsey
parent
09a0c95620
commit
fab44d035a
|
|
@ -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);
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
interpolate_point4d_spheroid(
|
||||
const POINT4D *p1, const POINT4D *p2, POINT4D *p,
|
||||
const SPHEROID *s, double f)
|
||||
|
||||
{
|
||||
/* 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);
|
||||
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;
|
||||
|
||||
/* 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);
|
||||
}
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* @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;
|
||||
|
||||
/* 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);
|
||||
|
||||
/* Loop and sum the length for each segment */
|
||||
for ( i = 1; i < pa->npoints; i++ )
|
||||
/* Special sphere case */
|
||||
if ( s == NULL || s->a == s->b )
|
||||
{
|
||||
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);
|
||||
/* 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);
|
||||
}
|
||||
|
||||
/* If success, use newly computed lat and lon,
|
||||
* otherwise return precomputed cartesian result */
|
||||
if (success == LW_SUCCESS)
|
||||
{
|
||||
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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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 )
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
Loading…
Reference in a new issue