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) +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; } 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)