[postgis-tickets] [SCM] PostGIS branch master updated. 3.3.0rc2-1205-gf44fd54f9
git at osgeo.org
git at osgeo.org
Tue Jul 11 08:29:56 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 f44fd54f97146dd3cef22cbb9901ff887438c7f4 (commit)
from e7b473841608110b1263f553612405657e041cb4 (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 f44fd54f97146dd3cef22cbb9901ff887438c7f4
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Tue Jul 11 08:25:13 2023 -0700
Add linear referencing support functions for the geography
type. Submitted by MobilityDB, @estebanzimanyi.
* ST_LineInterpolatePoint(geography, float)
* ST_LineInterpolatePoints(geography, float)
* ST_LineLocatePoint(geography, geography)
* ST_LineSubstring(geography, float, float)
* ST_ClosestPoint(geography, geography)
* ST_ShortestLine(geography, geography)
diff --git a/doc/reference_lrs.xml b/doc/reference_lrs.xml
index 5ce26bfe2..e28f23fa5 100644
--- a/doc/reference_lrs.xml
+++ b/doc/reference_lrs.xml
@@ -16,6 +16,14 @@
<paramdef><type>geometry </type> <parameter>a_linestring</parameter></paramdef>
<paramdef><type>float8 </type> <parameter>a_fraction</parameter></paramdef>
</funcprototype>
+
+ <funcprototype>
+ <funcdef>geography <function>ST_LineInterpolatePoint</function></funcdef>
+ <paramdef><type>geography </type> <parameter>a_linestring</parameter></paramdef>
+ <paramdef><type>float8 </type> <parameter>a_fraction</parameter></paramdef>
+ <paramdef choice="opt"><type>boolean </type><parameter>use_spheroid = true</parameter></paramdef>
+ </funcprototype>
+
</funcsynopsis>
</refsynopsisdiv>
@@ -183,6 +191,15 @@ SELECT ST_AsText(
<paramdef><type>float8 </type> <parameter>a_fraction</parameter></paramdef>
<paramdef><type>boolean </type> <parameter>repeat</parameter></paramdef>
</funcprototype>
+
+ <funcprototype>
+ <funcdef>geography <function>ST_LineInterpolatePoints</function></funcdef>
+ <paramdef><type>geography </type> <parameter>a_linestring</parameter></paramdef>
+ <paramdef><type>float8 </type> <parameter>a_fraction</parameter></paramdef>
+ <paramdef choice="opt"><type>boolean </type><parameter>use_spheroid = true</parameter></paramdef>
+ <paramdef choice="opt"><type>boolean </type> <parameter>repeat = true</parameter></paramdef>
+ </funcprototype>
+
</funcsynopsis>
</refsynopsisdiv>
@@ -245,11 +262,20 @@ SELECT ST_AsText(ST_LineInterpolatePoints('LINESTRING(25 50, 100 125, 150 190)',
<refsynopsisdiv>
<funcsynopsis>
+
<funcprototype>
<funcdef>float8 <function>ST_LineLocatePoint</function></funcdef>
<paramdef><type>geometry </type> <parameter>a_linestring</parameter></paramdef>
<paramdef><type>geometry </type> <parameter>a_point</parameter></paramdef>
</funcprototype>
+
+ <funcprototype>
+ <funcdef>float8 <function>ST_LineLocatePoint</function></funcdef>
+ <paramdef><type>geography </type> <parameter>a_linestring</parameter></paramdef>
+ <paramdef><type>geography </type> <parameter>a_point</parameter></paramdef>
+ <paramdef choice="opt"><type>boolean </type><parameter>use_spheroid = true</parameter></paramdef>
+ </funcprototype>
+
</funcsynopsis>
</refsynopsisdiv>
@@ -330,6 +356,14 @@ FROM (SELECT ST_GeomFromText('LINESTRING(1 2, 4 5, 6 7)') As the_line) As foo;
<paramdef><type>float8 </type> <parameter>startfraction</parameter></paramdef>
<paramdef><type>float8 </type> <parameter>endfraction</parameter></paramdef>
</funcprototype>
+
+ <funcprototype>
+ <funcdef>geography <function>ST_LineSubstring</function></funcdef>
+ <paramdef><type>geography </type> <parameter>a_linestring</parameter></paramdef>
+ <paramdef><type>float8 </type> <parameter>startfraction</parameter></paramdef>
+ <paramdef><type>float8 </type> <parameter>endfraction</parameter></paramdef>
+ </funcprototype>
+
</funcsynopsis>
</refsynopsisdiv>
diff --git a/doc/reference_measure.xml b/doc/reference_measure.xml
index 4bd94c94d..21d5c9811 100644
--- a/doc/reference_measure.xml
+++ b/doc/reference_measure.xml
@@ -294,15 +294,25 @@ SELECT degrees( ST_Angle('LINESTRING(0 0, 0.3 0.7, 1 1)', 'LINESTRING(0 0, 0.2 0
<refsynopsisdiv>
<funcsynopsis>
+
<funcprototype>
<funcdef>geometry <function>ST_ClosestPoint</function></funcdef>
-
<paramdef><type>geometry </type>
<parameter>geom1</parameter></paramdef>
-
<paramdef><type>geometry </type>
<parameter>geom2</parameter></paramdef>
</funcprototype>
+
+ <funcprototype>
+ <funcdef>geography <function>ST_ClosestPoint</function></funcdef>
+ <paramdef><type>geography </type>
+ <parameter>geom1</parameter></paramdef>
+ <paramdef><type>geography </type>
+ <parameter>geom2</parameter></paramdef>
+ <paramdef choice="opt"><type>boolean </type>
+ <parameter>repeat = true</parameter></paramdef>
+ </funcprototype>
+
</funcsynopsis>
</refsynopsisdiv>
@@ -1879,13 +1889,21 @@ FROM ST_GeogFromText('MULTIPOLYGON(((-71.1044543107478 42.340674480411,-71.10445
<funcsynopsis>
<funcprototype>
<funcdef>geometry <function>ST_ShortestLine</function></funcdef>
-
<paramdef><type>geometry </type>
<parameter>geom1</parameter></paramdef>
-
<paramdef><type>geometry </type>
<parameter>geom2</parameter></paramdef>
</funcprototype>
+
+ <funcprototype>
+ <funcdef>geography <function>ST_ShortestLine</function></funcdef>
+ <paramdef><type>geography </type>
+ <parameter>geom1</parameter></paramdef>
+ <paramdef><type>geography </type>
+ <parameter>geom2</parameter></paramdef>
+ <paramdef choice="opt"><type>boolean </type><parameter>use_spheroid = true</parameter></paramdef>
+ </funcprototype>
+
</funcsynopsis>
</refsynopsisdiv>
diff --git a/liblwgeom/Makefile.in b/liblwgeom/Makefile.in
index 9f0530c8f..3094837af 100644
--- a/liblwgeom/Makefile.in
+++ b/liblwgeom/Makefile.in
@@ -105,6 +105,7 @@ SA_OBJS = \
gserialized1.o \
gserialized2.o \
lwgeodetic.o \
+ lwgeodetic_measures.o \
lwgeodetic_tree.o \
lwrandom.o \
lwtree.o \
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index cbac60f3c..a27529e6c 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -1791,6 +1791,26 @@ extern double lwgeom_length_spheroid(const LWGEOM *geom, const SPHEROID *s);
*/
extern int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2);
+
+extern LWGEOM * geography_substring(const LWLINE *line,
+ double from, double to,
+ double tolerance);
+
+extern LWGEOM * geography_interpolate_points(
+ const LWLINE *line,
+ double length_fraction,
+ const SPHEROID *s, char repeat);
+
+extern double
+ptarray_locate_point_spheroid(
+ const POINTARRAY *pa,
+ const POINT4D *p4d,
+ const SPHEROID *s,
+ double tolerance,
+ double *mindistout,
+ POINT4D *proj4d);
+
+
typedef struct {
POINT2D* center;
double radius;
diff --git a/liblwgeom/lwgeodetic_measures.c b/liblwgeom/lwgeodetic_measures.c
new file mode 100644
index 000000000..d8e245657
--- /dev/null
+++ b/liblwgeom/lwgeodetic_measures.c
@@ -0,0 +1,561 @@
+/*****************************************************************************
+ *
+ * This MobilityDB code is provided under The PostgreSQL License.
+ * Copyright (c) 2016-2023, Université libre de Bruxelles and MobilityDB
+ * contributors
+ *
+ * MobilityDB includes portions of PostGIS version 3 source code released
+ * under the GNU General Public License (GPLv2 or later).
+ * Copyright (c) 2001-2023, PostGIS contributors
+ *
+ * Permission to use, copy, modify, and distribute this software and its
+ * documentation for any purpose, without fee, and without a written
+ * agreement is hereby granted, provided that the above copyright notice and
+ * this paragraph and the following two paragraphs appear in all copies.
+ *
+ * IN NO EVENT SHALL UNIVERSITE LIBRE DE BRUXELLES BE LIABLE TO ANY PARTY FOR
+ * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING
+ * LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION,
+ * EVEN IF UNIVERSITE LIBRE DE BRUXELLES HAS BEEN ADVISED OF THE POSSIBILITY
+ * OF SUCH DAMAGE.
+ *
+ * UNIVERSITE LIBRE DE BRUXELLES SPECIFICALLY DISCLAIMS ANY WARRANTIES,
+ * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
+ * AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON
+ * AN "AS IS" BASIS, AND UNIVERSITE LIBRE DE BRUXELLES HAS NO OBLIGATIONS TO
+ * PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+ *
+ *****************************************************************************/
+
+/**
+ * @file
+ * @brief Spatial functions for PostGIS geography.
+ *
+ * These functions are supposed to be included in a forthcoming version of
+ * PostGIS, to be proposed as a PR. This still remains to be done.
+ * These functions are not needed in MobilityDB.
+ */
+
+// #include "point/geography_funcs.h"
+
+/* C */
+#include <float.h>
+/* PostGIS */
+#include <liblwgeom_internal.h>
+#include <lwgeodetic_tree.h>
+
+/*****************************************************************************/
+
+
+/***********************************************************************
+ * ST_LineSubstring for geographies
+ ***********************************************************************/
+
+/**
+ * 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);
+}
+
+
+
+/**
+ * @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++ )
+ {
+ 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);
+ }
+ return length;
+}
+
+/**
+ * @brief Return the part of a line between two fractional locations.
+ */
+LWGEOM *
+geography_substring(
+ const LWLINE *lwline,
+ double from, double to,
+ double tolerance)
+{
+ POINTARRAY *dpa;
+ POINTARRAY *ipa = lwline->points;
+ LWGEOM *lwresult;
+ POINT4D pt;
+ POINT4D p1, p2;
+ POINT3D q1, q2;
+ GEOGRAPHIC_POINT g1, g2;
+ int nsegs, i;
+ double length, slength, tlength;
+ int state = 0; /* 0 = before, 1 = inside */
+ uint32_t srid = lwline->srid;
+
+ /*
+ * Create a dynamic pointarray with an initial capacity
+ * equal to full copy of input points
+ */
+ dpa = ptarray_construct_empty(
+ (char) FLAGS_GET_Z(ipa->flags),
+ (char) FLAGS_GET_M(ipa->flags),
+ ipa->npoints);
+
+ /* Compute total line length */
+ length = ptarray_length_sphere(ipa);
+
+ /* Get 'from' and 'to' lengths */
+ from = length * from;
+ to = length * to;
+ tlength = 0;
+ getPoint4d_p(ipa, 0, &p1);
+ geographic_point_init(p1.x, p1.y, &g1);
+ nsegs = ipa->npoints - 1;
+ for (i = 0; i < nsegs; i++)
+ {
+ double dseg;
+ getPoint4d_p(ipa, (uint32_t) i+1, &p2);
+ geographic_point_init(p2.x, p2.y, &g2);
+
+ /* Find the length of this segment */
+ slength = sphere_distance(&g1, &g2);
+
+ /* We are before requested start. */
+ if (state == 0) /* before */
+ {
+ if (fabs ( from - ( tlength + slength ) ) <= tolerance)
+ {
+ /* Second point is our start */
+ ptarray_append_point(dpa, &p2, LW_FALSE);
+ state = 1; /* we're inside now */
+ goto END;
+ }
+ else if (fabs(from - tlength) <= tolerance)
+ {
+ /* First point is our start */
+ ptarray_append_point(dpa, &p1, LW_FALSE);
+ /*
+ * We're inside now, but will check
+ * 'to' point as well
+ */
+ state = 1;
+ }
+ /*
+ * Didn't reach the 'from' point,
+ * nothing to do
+ */
+ else if (from > tlength + slength)
+ {
+ goto END;
+ }
+ else /* tlength < from < tlength+slength */
+ {
+ /* 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);
+ ptarray_append_point(dpa, &pt, LW_FALSE);
+ /* We're inside now, but will check 'to' point as well */
+ state = 1;
+ }
+ }
+
+ if (state == 1) /* inside */
+ {
+ /* 'to' point is our second point. */
+ if (fabs(to - ( tlength + slength ) ) <= tolerance )
+ {
+ ptarray_append_point(dpa, &p2, LW_FALSE);
+ break; /* substring complete */
+ }
+ /*
+ * 'to' point is our first point.
+ * (should only happen if 'to' is 0)
+ */
+ else if (fabs(to - tlength) <= tolerance)
+ {
+ ptarray_append_point(dpa, &p1, LW_FALSE);
+ break; /* substring complete */
+ }
+ /*
+ * Didn't reach the 'end' point,
+ * just copy second point
+ */
+ else if (to > tlength + slength)
+ {
+ ptarray_append_point(dpa, &p2, LW_FALSE);
+ goto END;
+ }
+ /*
+ * 'to' point falls on this segment
+ * Interpolate and break.
+ */
+ else if (to < tlength + slength )
+ {
+ dseg = (to - tlength) / slength;
+ geog2cart(&g1, &q1);
+ geog2cart(&g2, &q2);
+ interpolate_point4d_sphere(&q1, &q2, &p1, &p2, dseg, &pt);
+ ptarray_append_point(dpa, &pt, LW_FALSE);
+ break;
+ }
+ }
+
+ END:
+ tlength += slength;
+ memcpy(&p1, &p2, sizeof(POINT4D));
+ }
+
+ if (dpa->npoints <= 1) {
+ lwresult = lwpoint_as_lwgeom(lwpoint_construct(srid, NULL, dpa));
+ }
+ else {
+ lwresult = lwline_as_lwgeom(lwline_construct(srid, NULL, dpa));
+ }
+
+ return lwresult;
+}
+
+
+/***********************************************************************
+ * Interpolate a point along a geographic line.
+ ***********************************************************************/
+
+/**
+ * @brief Interpolate a point along a geographic line.
+ */
+LWGEOM *
+geography_interpolate_points(
+ const LWLINE *line, double length_fraction,
+ const SPHEROID *s, char repeat)
+{
+ POINT4D pt;
+ uint32_t i;
+ uint32_t points_to_interpolate;
+ uint32_t points_found = 0;
+ double length;
+ double length_fraction_increment = length_fraction;
+ double length_fraction_consumed = 0;
+ char has_z = (char) lwgeom_has_z(lwline_as_lwgeom(line));
+ char has_m = (char) lwgeom_has_m(lwline_as_lwgeom(line));
+ const POINTARRAY* ipa = line->points;
+ POINTARRAY *opa;
+ POINT4D p1, p2;
+ POINT3D q1, q2;
+ LWGEOM *lwresult;
+ GEOGRAPHIC_POINT g1, g2;
+ uint32_t srid = line->srid;
+
+ /* Empty.InterpolatePoint == Point Empty */
+ if ( lwline_is_empty(line) )
+ {
+ return lwgeom_clone_deep(lwline_as_lwgeom(line));
+ }
+
+ /*
+ * If distance is one of the two extremes, return the point on that
+ * end rather than doing any computations
+ */
+ if ( length_fraction == 0.0 || length_fraction == 1.0 )
+ {
+ if ( length_fraction == 0.0 )
+ getPoint4d_p(ipa, 0, &pt);
+ else
+ getPoint4d_p(ipa, ipa->npoints-1, &pt);
+
+ opa = ptarray_construct(has_z, has_m, 1);
+ ptarray_set_point4d(opa, 0, &pt);
+ return lwpoint_as_lwgeom(lwpoint_construct(srid, NULL, opa));
+ }
+
+ /* Interpolate points along the line */
+ length = ptarray_length_spheroid(ipa, s);
+ points_to_interpolate = repeat ? (uint32_t) floor(1 / length_fraction) : 1;
+ opa = ptarray_construct(has_z, has_m, points_to_interpolate);
+
+ getPoint4d_p(ipa, 0, &p1);
+ geographic_point_init(p1.x, p1.y, &g1);
+ for ( i = 0; i < ipa->npoints - 1 && points_found < points_to_interpolate; i++ )
+ {
+ getPoint4d_p(ipa, i+1, &p2);
+ geographic_point_init(p2.x, p2.y, &g2);
+ double 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
+ * segment.
+ */
+ while ( length_fraction < length_fraction_consumed + segment_length_frac && points_found < points_to_interpolate )
+ {
+ 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);
+ ptarray_set_point4d(opa, points_found++, &pt);
+ length_fraction += length_fraction_increment;
+ }
+
+ length_fraction_consumed += segment_length_frac;
+
+ p1 = p2;
+ g1 = g2;
+ }
+
+ /* Return the last point on the line. This shouldn't happen, but
+ * could if there's some floating point rounding errors. */
+ if (points_found < points_to_interpolate)
+ {
+ getPoint4d_p(ipa, ipa->npoints - 1, &pt);
+ ptarray_set_point4d(opa, points_found, &pt);
+ }
+
+ if (opa->npoints <= 1)
+ {
+ lwresult = lwpoint_as_lwgeom(lwpoint_construct(srid, NULL, opa));
+ } else {
+ lwresult = lwmpoint_as_lwgeom(lwmpoint_construct(srid, opa));
+ }
+
+ return lwresult;
+}
+
+/**
+ * @brief Locate a point along the point array defining a geographic line.
+ */
+double
+ptarray_locate_point_spheroid(
+ const POINTARRAY *pa,
+ const POINT4D *p4d,
+ const SPHEROID *s,
+ double tolerance,
+ double *mindistout,
+ POINT4D *proj4d)
+{
+ GEOGRAPHIC_EDGE e;
+ GEOGRAPHIC_POINT a, b, nearest = {0}; /* make compiler quiet */
+ POINT4D p1, p2;
+ const POINT2D *p;
+ POINT2D proj;
+ uint32_t i, seg = 0;
+ int use_sphere = (s->a == s->b ? 1 : 0);
+ int hasz;
+ double za = 0.0, zb = 0.0;
+ double distance,
+ length, /* Used for computing lengths */
+ seglength = 0.0, /* length of the segment where the closest point is located */
+ partlength = 0.0, /* length from the beginning of the point array to the closest point */
+ totlength = 0.0; /* length of the point array */
+
+ /* Initialize our point */
+ geographic_point_init(p4d->x, p4d->y, &a);
+
+ /* Handle point/point case here */
+ if ( pa->npoints <= 1)
+ {
+ if ( pa->npoints == 1 )
+ {
+ p = getPoint2d_cp(pa, 0);
+ geographic_point_init(p->x, p->y, &b);
+ /* Sphere special case, axes equal */
+ *mindistout = s->radius * sphere_distance(&a, &b);
+ /* If close or greater than tolerance, get the real answer to be sure */
+ if ( ! use_sphere || *mindistout > 0.95 * tolerance )
+ {
+ *mindistout = spheroid_distance(&a, &b, s);
+ }
+ }
+ return 0.0;
+ }
+
+ /* Make result really big, so that everything will be smaller than it */
+ distance = FLT_MAX;
+
+ /* Initialize start of line */
+ p = getPoint2d_cp(pa, 0);
+ geographic_point_init(p->x, p->y, &(e.start));
+
+ /* Iterate through the edges in our line */
+ for ( i = 1; i < pa->npoints; i++ )
+ {
+ double d;
+ p = getPoint2d_cp(pa, i);
+ geographic_point_init(p->x, p->y, &(e.end));
+ /* Get the spherical distance between point and edge */
+ d = s->radius * edge_distance_to_point(&e, &a, &b);
+ /* New shortest distance! Record this distance / location / segment */
+ if ( d < distance )
+ {
+ distance = d;
+ nearest = b;
+ seg = i - 1;
+ }
+ /* We've gotten closer than the tolerance... */
+ if ( d < tolerance )
+ {
+ /* Working on a sphere? The answer is correct, return */
+ if ( use_sphere )
+ {
+ break;
+ }
+ /* Far enough past the tolerance that the spheroid calculation won't change things */
+ else if ( d < tolerance * 0.95 )
+ {
+ break;
+ }
+ /* On a spheroid and near the tolerance? Confirm that we are *actually* closer than tolerance */
+ else
+ {
+ d = spheroid_distance(&a, &nearest, s);
+ /* Yes, closer than tolerance, return! */
+ if ( d < tolerance )
+ break;
+ }
+ }
+ e.start = e.end;
+ }
+
+ if ( mindistout ) *mindistout = distance;
+
+ /* See if we have a third dimension */
+ hasz = (bool) FLAGS_GET_Z(pa->flags);
+
+ /* Initialize first point of array */
+ getPoint4d_p(pa, 0, &p1);
+ geographic_point_init(p1.x, p1.y, &a);
+ if ( hasz )
+ za = p1.z;
+
+ /* Loop and sum the length for each segment */
+ for ( i = 1; i < pa->npoints; i++ )
+ {
+ getPoint4d_p(pa, i, &p1);
+ geographic_point_init(p1.x, p1.y, &b);
+ if ( hasz )
+ zb = p1.z;
+
+ /* Special sphere case */
+ if ( s->a == s->b )
+ length = s->radius * sphere_distance(&a, &b);
+ /* Spheroid case */
+ else
+ length = spheroid_distance(&a, &b, s);
+
+ /* Add in the vertical displacement if we're in 3D */
+ if ( hasz )
+ length = sqrt( (zb-za)*(zb-za) + length*length );
+
+ /* Add this segment length to the total length */
+ totlength += length;
+
+ /* Add this segment length to the partial length */
+ if (i - 1 < seg)
+ partlength += length;
+ else if (i - 1 == seg)
+ /* Save segment length for computing the final value of partlength */
+ seglength = length;
+
+ /* B gets incremented in the next loop, so we save the value here */
+ a = b;
+ za = zb;
+ }
+
+ /* Copy nearest into 2D/4D holder */
+ proj4d->x = proj.x = rad2deg(nearest.lon);
+ proj4d->y = proj.y = rad2deg(nearest.lat);
+
+ /* Compute distance from beginning of the segment to closest point */
+
+ /* Start of the segment */
+ getPoint4d_p(pa, seg, &p1);
+ geographic_point_init(p1.x, p1.y, &a);
+
+ /* Closest point */
+ geographic_point_init(proj4d->x, proj4d->y, &b);
+
+ /* Special sphere case */
+ if ( s->a == s->b )
+ length = s->radius * sphere_distance(&a, &b);
+ /* Spheroid case */
+ else
+ length = spheroid_distance(&a, &b, s);
+
+ if ( hasz )
+ {
+ /* Compute Z and M values for closest point */
+ double f = length / seglength;
+ getPoint4d_p(pa, seg + 1, &p2);
+ proj4d->z = p1.z + ((p2.z - p1.z) * f);
+ proj4d->m = p1.m + ((p2.m - p1.m) * f);
+ /* Add in the vertical displacement if we're in 3D */
+ za = p1.z;
+ zb = proj4d->z;
+ length = sqrt( (zb-za)*(zb-za) + length*length );
+ }
+
+ /* Add this segment length to the total */
+ partlength += length;
+
+ /* Location of any point on a zero-length line is 0 */
+ /* See http://trac.osgeo.org/postgis/ticket/1772#comment:2 */
+ if ( partlength == 0 || totlength == 0 )
+ return 0.0;
+
+ /* For robustness, force 0 when closest point == startpoint */
+ p = getPoint2d_cp(pa, 0);
+ if ( seg == 0 && p2d_same(&proj, p) )
+ return 0.0;
+
+ /* For robustness, force 1 when closest point == endpoint */
+ p = getPoint2d_cp(pa, pa->npoints - 1);
+ if ( (seg >= (pa->npoints-2)) && p2d_same(&proj, p) )
+ return 1.0;
+
+ return partlength / totlength;
+}
+
+/*****************************************************************************/
diff --git a/liblwgeom/lwgeodetic_tree.c b/liblwgeom/lwgeodetic_tree.c
index 9baa8cd3b..a425046a0 100644
--- a/liblwgeom/lwgeodetic_tree.c
+++ b/liblwgeom/lwgeodetic_tree.c
@@ -674,7 +674,7 @@ circ_internal_nodes_sort(CIRC_NODE **nodes, uint32_t num_nodes, const CIRC_NODE
/***********************************************************************/
-static double
+double
circ_tree_distance_tree_internal(const CIRC_NODE* n1, const CIRC_NODE* n2, double threshold, double* min_dist, double* max_dist, GEOGRAPHIC_POINT* closest1, GEOGRAPHIC_POINT* closest2)
{
double max;
@@ -1026,3 +1026,83 @@ lwgeom_calculate_circ_tree(const LWGEOM* lwgeom)
}
}
+
+
+/***********************************************************************
+ * Closest point and closest line functions for geographies.
+ ***********************************************************************/
+
+LWGEOM *
+geography_tree_closestpoint(const LWGEOM* lwgeom1, const LWGEOM* lwgeom2, double threshold)
+{
+ CIRC_NODE* circ_tree1 = NULL;
+ CIRC_NODE* circ_tree2 = NULL;
+ double min_dist = FLT_MAX;
+ double max_dist = FLT_MAX;
+ GEOGRAPHIC_POINT closest1, closest2;
+ LWGEOM *result;
+ POINT4D p;
+
+ circ_tree1 = lwgeom_calculate_circ_tree(lwgeom1);
+ circ_tree2 = lwgeom_calculate_circ_tree(lwgeom2);
+
+ /* Quietly decrease the threshold just a little to avoid cases where */
+ /* the actual spheroid distance is larger than the sphere distance */
+ /* causing the return value to be larger than the threshold value */
+ // double threshold_radians = 0.95 * threshold / spheroid->radius;
+ double threshold_radians = threshold / WGS84_RADIUS;
+
+ circ_tree_distance_tree_internal(
+ circ_tree1, circ_tree2, threshold_radians,
+ &min_dist, &max_dist, &closest1, &closest2);
+
+ p.x = rad2deg(closest1.lon);
+ p.y = rad2deg(closest1.lat);
+ result = (LWGEOM *)lwpoint_make2d(lwgeom_get_srid(lwgeom1), p.x, p.y);
+
+ circ_tree_free(circ_tree1);
+ circ_tree_free(circ_tree2);
+ return result;
+}
+
+
+LWGEOM *
+geography_tree_shortestline(const LWGEOM* lwgeom1, const LWGEOM* lwgeom2, double threshold, const SPHEROID *spheroid)
+{
+ CIRC_NODE* circ_tree1 = NULL;
+ CIRC_NODE* circ_tree2 = NULL;
+ double min_dist = FLT_MAX;
+ double max_dist = FLT_MAX;
+ GEOGRAPHIC_POINT closest1, closest2;
+ LWGEOM *geoms[2];
+ LWGEOM *result;
+ POINT4D p1, p2;
+ uint32_t srid = lwgeom1->srid;
+
+ circ_tree1 = lwgeom_calculate_circ_tree(lwgeom1);
+ circ_tree2 = lwgeom_calculate_circ_tree(lwgeom2);
+
+ /* Quietly decrease the threshold just a little to avoid cases where */
+ /* the actual spheroid distance is larger than the sphere distance */
+ /* causing the return value to be larger than the threshold value */
+ // double threshold_radians = 0.95 * threshold / spheroid->radius;
+ double threshold_radians = threshold / spheroid->radius;
+
+ circ_tree_distance_tree_internal(circ_tree1, circ_tree2, threshold_radians,
+ &min_dist, &max_dist, &closest1, &closest2);
+
+ p1.x = rad2deg(closest1.lon);
+ p1.y = rad2deg(closest1.lat);
+ p2.x = rad2deg(closest2.lon);
+ p2.y = rad2deg(closest2.lat);
+
+ geoms[0] = (LWGEOM *)lwpoint_make2d(srid, p1.x, p1.y);
+ geoms[1] = (LWGEOM *)lwpoint_make2d(srid, p2.x, p2.y);
+ result = (LWGEOM *)lwline_from_lwgeom_array(srid, 2, geoms);
+
+ lwgeom_free(geoms[0]);
+ lwgeom_free(geoms[1]);
+ circ_tree_free(circ_tree1);
+ circ_tree_free(circ_tree2);
+ return result;
+}
diff --git a/liblwgeom/lwgeodetic_tree.h b/liblwgeom/lwgeodetic_tree.h
index c367ebf16..90987a6cf 100644
--- a/liblwgeom/lwgeodetic_tree.h
+++ b/liblwgeom/lwgeodetic_tree.h
@@ -56,6 +56,10 @@ CIRC_NODE* lwgeom_calculate_circ_tree(const LWGEOM* lwgeom);
int circ_tree_get_point(const CIRC_NODE* node, POINT2D* pt);
int circ_tree_get_point_outside(const CIRC_NODE* node, POINT2D* pt);
+LWGEOM * geography_tree_closestpoint(const LWGEOM* g1, const LWGEOM* g2, double threshold);
+LWGEOM * geography_tree_shortestline(const LWGEOM* g1, const LWGEOM* g2, double threshold, const SPHEROID *spheroid);
+
+
#endif /* _LWGEODETIC_TREE_H */
diff --git a/postgis/geography.h b/postgis/geography.h
index 2cc925ff6..e81547e71 100644
--- a/postgis/geography.h
+++ b/postgis/geography.h
@@ -37,3 +37,4 @@ void geography_valid_type(uint8_t type);
/* Expand the embedded bounding box in a #GSERIALIZED */
GSERIALIZED* gserialized_expand(GSERIALIZED *g, double distance);
+
diff --git a/postgis/geography.sql.in b/postgis/geography.sql.in
index 9a61d7a7a..016559adb 100644
--- a/postgis/geography.sql.in
+++ b/postgis/geography.sql.in
@@ -867,4 +867,76 @@ CREATE OR REPLACE FUNCTION ST_Intersects(text, text)
-----------------------------------------------------------------------------
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_ClosestPoint(geography, geography, use_spheroid boolean DEFAULT true)
+ RETURNS geography
+ AS 'MODULE_PATHNAME', 'geography_closestpoint'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_ClosestPoint(text, text)
+ RETURNS geometry AS
+ $$ SELECT @extschema at .ST_ClosestPoint($1::@extschema at .geometry, $2::@extschema at .geometry); $$
+ LANGUAGE 'sql' IMMUTABLE PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_ShortestLine(geography, geography, use_spheroid boolean DEFAULT true)
+ RETURNS geography
+ AS 'MODULE_PATHNAME', 'geography_shortestline'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_ShortestLine(text, text)
+ RETURNS geometry AS
+ $$ SELECT @extschema at .ST_ShortestLine($1::@extschema at .geometry, $2::@extschema at .geometry); $$
+ LANGUAGE 'sql' IMMUTABLE PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_LineSubstring(geography, float8, float8)
+ RETURNS geography
+ AS 'MODULE_PATHNAME', 'geography_line_substring'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_LineSubstring(text, float8, float8)
+ RETURNS geometry AS
+ $$ SELECT @extschema at .ST_LineSubstring($1::@extschema at .geometry, $2, $3); $$
+ LANGUAGE 'sql' IMMUTABLE PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_LineLocatePoint(geography, geography, use_spheroid boolean DEFAULT true)
+ RETURNS float
+ AS 'MODULE_PATHNAME', 'geography_line_locate_point'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_LineLocatePoint(text, text)
+ RETURNS float AS
+ $$ SELECT @extschema at .ST_LineLocatePoint($1::@extschema at .geometry, $2::@extschema at .geometry); $$
+ LANGUAGE 'sql' IMMUTABLE PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_LineInterpolatePoints(geography, float8, use_spheroid boolean DEFAULT true, repeat boolean DEFAULT true)
+ RETURNS geography
+ AS 'MODULE_PATHNAME', 'geography_line_interpolate_point'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_LineInterpolatePoints(text, float8)
+ RETURNS geometry AS
+ $$ SELECT @extschema at .ST_LineInterpolatePoints($1::@extschema at .geometry, $2); $$
+ LANGUAGE 'sql' IMMUTABLE PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_LineInterpolatePoint(geography, float8, use_spheroid boolean DEFAULT true)
+ RETURNS geography
+ AS 'MODULE_PATHNAME', 'geography_line_interpolate_point'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE;
+
+-- Availability: 3.4.0
+CREATE OR REPLACE FUNCTION ST_LineInterpolatePoint(text, float8)
+ RETURNS geometry AS
+ $$ SELECT @extschema at .ST_LineInterpolatePoint($1::@extschema at .geometry, $2); $$
+ LANGUAGE 'sql' IMMUTABLE PARALLEL SAFE;
+----------------------------------------------------------------
diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index 6bda676b3..e52671f4f 100644
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -66,6 +66,12 @@ Datum geography_project(PG_FUNCTION_ARGS);
Datum geography_azimuth(PG_FUNCTION_ARGS);
Datum geography_segmentize(PG_FUNCTION_ARGS);
+Datum geography_line_locate_point(PG_FUNCTION_ARGS);
+Datum geography_line_interpolate_point(PG_FUNCTION_ARGS);
+Datum geography_line_substring(PG_FUNCTION_ARGS);
+Datum geography_closestpoint(PG_FUNCTION_ARGS);
+Datum geography_shortestline(PG_FUNCTION_ARGS);
+
PG_FUNCTION_INFO_V1(geography_distance_knn);
Datum geography_distance_knn(PG_FUNCTION_ARGS)
@@ -135,17 +141,13 @@ Datum geography_distance_uncached(PG_FUNCTION_ARGS)
{
LWGEOM *lwgeom1 = NULL;
LWGEOM *lwgeom2 = NULL;
- GSERIALIZED *g1 = NULL;
- GSERIALIZED *g2 = NULL;
+ GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED *g2 = PG_GETARG_GSERIALIZED_P(1);
double distance;
double tolerance = FP_TOLERANCE;
bool use_spheroid = true;
SPHEROID s;
- /* Get our geometry objects loaded into memory. */
- g1 = PG_GETARG_GSERIALIZED_P(0);
- g2 = PG_GETARG_GSERIALIZED_P(1);
-
/* Read our tolerance value. */
if ( PG_NARGS() > 2 && ! PG_ARGISNULL(2) )
tolerance = PG_GETARG_FLOAT8(2);
@@ -167,7 +169,7 @@ Datum geography_distance_uncached(PG_FUNCTION_ARGS)
lwgeom2 = lwgeom_from_gserialized(g2);
/* Return NULL on empty arguments. */
- if ( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
+ if ( !lwgeom1 || !lwgeom2 || lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
{
PG_FREE_IF_COPY(g1, 0);
PG_FREE_IF_COPY(g2, 1);
@@ -180,8 +182,6 @@ Datum geography_distance_uncached(PG_FUNCTION_ARGS)
distance = lwgeom_distance_spheroid(lwgeom1, lwgeom2, &s, tolerance);
- POSTGIS_DEBUGF(2, "[GIST] '%s' got distance %g", __func__, distance);
-
/* Clean up */
lwgeom_free(lwgeom1);
lwgeom_free(lwgeom2);
@@ -213,7 +213,6 @@ Datum geography_distance(PG_FUNCTION_ARGS)
bool use_spheroid = true;
SPHEROID s;
-
if (PG_NARGS() > 2)
use_spheroid = PG_GETARG_BOOL(2);
@@ -1107,18 +1106,14 @@ Datum geography_project_geography(PG_FUNCTION_ARGS)
PG_FUNCTION_INFO_V1(geography_azimuth);
Datum geography_azimuth(PG_FUNCTION_ARGS)
{
+ GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED *g2 = PG_GETARG_GSERIALIZED_P(1);
LWGEOM *lwgeom1 = NULL;
LWGEOM *lwgeom2 = NULL;
- GSERIALIZED *g1 = NULL;
- GSERIALIZED *g2 = NULL;
double azimuth;
SPHEROID s;
uint32_t type1, type2;
- /* Get our geometry object loaded into memory. */
- g1 = PG_GETARG_GSERIALIZED_P(0);
- g2 = PG_GETARG_GSERIALIZED_P(1);
-
/* Only return for points. */
type1 = gserialized_get_type(g1);
type2 = gserialized_get_type(g2);
@@ -1154,7 +1149,7 @@ Datum geography_azimuth(PG_FUNCTION_ARGS)
PG_FREE_IF_COPY(g2, 1);
/* Return NULL for unknown (same point) azimuth */
- if( isnan(azimuth) )
+ if( !isfinite(azimuth) )
{
PG_RETURN_NULL();
}
@@ -1164,26 +1159,19 @@ Datum geography_azimuth(PG_FUNCTION_ARGS)
-/*
-** geography_segmentize(GSERIALIZED *g1, double max_seg_length)
-** returns densified geometry with no segment longer than max
-*/
+/**
+ * ST_Segmentize(geography, double max_seg_length)
+ * Returns densified geometry with no segment longer than maximum.
+ */
PG_FUNCTION_INFO_V1(geography_segmentize);
Datum geography_segmentize(PG_FUNCTION_ARGS)
{
LWGEOM *lwgeom1 = NULL;
LWGEOM *lwgeom2 = NULL;
- GSERIALIZED *g1 = NULL;
+ GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
GSERIALIZED *g2 = NULL;
- double max_seg_length;
- uint32_t type1;
-
- /* Get our geometry object loaded into memory. */
- g1 = PG_GETARG_GSERIALIZED_P(0);
- type1 = gserialized_get_type(g1);
-
- /* Convert max_seg_length from metric units to radians */
- max_seg_length = PG_GETARG_FLOAT8(1) / WGS84_RADIUS;
+ double max_seg_length = PG_GETARG_FLOAT8(1) / WGS84_RADIUS;
+ uint32_t type1 = gserialized_get_type(g1);
/* We can't densify points or points, reflect them back */
if ( type1 == POINTTYPE || type1 == MULTIPOINTTYPE || gserialized_is_empty(g1) )
@@ -1216,3 +1204,275 @@ Datum geography_segmentize(PG_FUNCTION_ARGS)
}
+/********************************************************************************/
+
+/**
+ * ST_LineSubstring(geography line, float start_fraction, float end_fraction)
+ * Return the part of a line between two fractional locations.
+ */
+PG_FUNCTION_INFO_V1(geography_line_substring);
+Datum geography_line_substring(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *gs = PG_GETARG_GSERIALIZED_P(0);
+ double from_fraction = PG_GETARG_FLOAT8(1);
+ double to_fraction = PG_GETARG_FLOAT8(2);
+ LWLINE *lwline;
+ LWGEOM *lwresult;
+ GSERIALIZED *result;
+
+ /* Return NULL on empty argument. */
+ if ( gserialized_is_empty(gs) )
+ {
+ PG_FREE_IF_COPY(gs, 0);
+ PG_RETURN_NULL();
+ }
+
+ if ( from_fraction < 0 || from_fraction > 1 )
+ {
+ elog(ERROR,"%s: second argument is not within [0,1]", __func__);
+ PG_FREE_IF_COPY(gs, 0);
+ PG_RETURN_NULL();
+ }
+ if ( to_fraction < 0 || to_fraction > 1 )
+ {
+ elog(ERROR,"%s: argument arg is not within [0,1]", __func__);
+ PG_FREE_IF_COPY(gs, 0);
+ PG_RETURN_NULL();
+ }
+ if ( from_fraction > to_fraction )
+ {
+ elog(ERROR, "%s: second argument must be smaller than third argument", __func__);
+ PG_RETURN_NULL();
+ }
+
+ lwline = lwgeom_as_lwline(lwgeom_from_gserialized(gs));
+ if ( !lwline )
+ {
+ elog(ERROR,"%s: first argument is not a line", __func__);
+ PG_FREE_IF_COPY(gs, 0);
+ PG_RETURN_NULL();
+ }
+
+ lwresult = geography_substring(lwline,
+ from_fraction, to_fraction, FP_TOLERANCE);
+
+ lwline_free(lwline);
+ PG_FREE_IF_COPY(gs, 0);
+ lwgeom_set_geodetic(lwresult, true);
+ result = geography_serialize(lwresult);
+ lwgeom_free(lwresult);
+
+ PG_RETURN_POINTER(result);
+}
+
+
+/**
+ * ST_LineInterpolatePoint(geography line, float fraction, boolean use_spheroid)
+ * Interpolate a point along a geographic line.
+ *
+ * ST_LineInterpolatePoints(geography line, float fraction, boolean use_spheroid)
+ * In-fill geographic line with multiple points using the fraction as the interval
+ * between each point.
+ */
+PG_FUNCTION_INFO_V1(geography_line_interpolate_point);
+Datum geography_line_interpolate_point(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *gs = PG_GETARG_GSERIALIZED_P(0);
+ double distance_fraction = PG_GETARG_FLOAT8(1);
+ /* Read calculation type */
+ bool use_spheroid = PG_GETARG_BOOL(2);
+ /* Read repeat mode */
+ bool repeat = (PG_NARGS() > 3) && PG_GETARG_BOOL(3);
+ LWLINE* lwline;
+ LWGEOM* lwresult;
+ SPHEROID s;
+ GSERIALIZED *result;
+
+ /* Return NULL on empty argument. */
+ if ( gserialized_is_empty(gs) )
+ {
+ PG_FREE_IF_COPY(gs, 0);
+ PG_RETURN_NULL();
+ }
+
+ if ( distance_fraction < 0 || distance_fraction > 1 )
+ {
+ elog(ERROR,"%s: second arg is not within [0,1]", __func__);
+ PG_FREE_IF_COPY(gs, 0);
+ PG_RETURN_NULL();
+ }
+
+ lwline = lwgeom_as_lwline(lwgeom_from_gserialized(gs));
+ if ( !lwline )
+ {
+ elog(ERROR,"%s: first arg is not a line", __func__);
+ PG_FREE_IF_COPY(gs, 0);
+ PG_RETURN_NULL();
+ }
+
+ /* Initialize spheroid */
+ // spheroid_init_from_srid(fcinfo, srid, &s);
+ spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
+
+ /* Set to sphere if requested */
+ if ( ! use_spheroid )
+ s.a = s.b = s.radius;
+
+ lwresult = geography_interpolate_points(lwline, distance_fraction, &s, repeat);
+
+ lwgeom_free(lwline_as_lwgeom(lwline));
+ PG_FREE_IF_COPY(gs, 0);
+
+ lwgeom_set_geodetic(lwresult, true);
+ result = geography_serialize(lwresult);
+ lwgeom_free(lwresult);
+
+ PG_RETURN_POINTER(result);
+}
+
+
+/**
+ * ST_LineLocatePoint(geography line, geography point, bool use_spheroid)
+ * Locate a point along a geographic line.
+ */
+PG_FUNCTION_INFO_V1(geography_line_locate_point);
+Datum geography_line_locate_point(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *gs1 = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED *gs2 = PG_GETARG_GSERIALIZED_P(1);
+ bool use_spheroid = PG_GETARG_BOOL(2);
+ double tolerance = FP_TOLERANCE;
+ SPHEROID s;
+ LWLINE *lwline;
+ LWPOINT *lwpoint;
+ POINTARRAY *pa;
+ POINT4D p, p_proj;
+ double ret;
+
+ gserialized_error_if_srid_mismatch(gs1, gs2, __func__);
+
+ /* Return NULL on empty argument. */
+ if ( gserialized_is_empty(gs1) || gserialized_is_empty(gs2))
+ {
+ PG_FREE_IF_COPY(gs1, 0);
+ PG_FREE_IF_COPY(gs2, 1);
+ PG_RETURN_NULL();
+ }
+
+ if ( gserialized_get_type(gs1) != LINETYPE )
+ {
+ elog(ERROR,"%s: 1st arg is not a line", __func__);
+ PG_RETURN_NULL();
+ }
+ if ( gserialized_get_type(gs2) != POINTTYPE )
+ {
+ elog(ERROR,"%s: 2st arg is not a point", __func__);
+ PG_RETURN_NULL();
+ }
+
+ /* Set to sphere if requested */
+ if ( ! use_spheroid ) {
+ s.a = s.b = s.radius;
+ }
+ else {
+ /* Initialize spheroid */
+ // spheroid_init_from_srid(fcinfo, srid, &s);
+ spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
+ }
+
+ lwline = lwgeom_as_lwline(lwgeom_from_gserialized(gs1));
+ lwpoint = lwgeom_as_lwpoint(lwgeom_from_gserialized(gs2));
+
+ pa = lwline->points;
+ lwpoint_getPoint4d_p(lwpoint, &p);
+
+ ret = ptarray_locate_point_spheroid(pa, &p, &s, tolerance, NULL, &p_proj);
+
+ PG_RETURN_FLOAT8(ret);
+}
+
+
+/**
+ * ST_ClosestPoint(geography line, geography point)
+ * Return the point in first input geography that is closest to the
+ * second input geography in 2d
+ */
+PG_FUNCTION_INFO_V1(geography_closestpoint);
+Datum geography_closestpoint(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED *g2 = PG_GETARG_GSERIALIZED_P(1);
+ LWGEOM *point, *lwg1, *lwg2;
+ GSERIALIZED *result;
+
+ gserialized_error_if_srid_mismatch(g1, g2, __func__);
+
+ lwg1 = lwgeom_from_gserialized(g1);
+ lwg2 = lwgeom_from_gserialized(g2);
+
+ /* Return NULL on empty/bad arguments. */
+ if ( !lwg1 || !lwg2 || lwgeom_is_empty(lwg1) || lwgeom_is_empty(lwg2) )
+ {
+ PG_FREE_IF_COPY(g1, 0);
+ PG_FREE_IF_COPY(g2, 1);
+ PG_RETURN_NULL();
+ }
+
+ point = geography_tree_closestpoint(lwg1, lwg2, FP_TOLERANCE);
+ result = geography_serialize(point);
+ lwgeom_free(point);
+ lwgeom_free(lwg1);
+ lwgeom_free(lwg2);
+
+ PG_FREE_IF_COPY(g1, 0);
+ PG_FREE_IF_COPY(g2, 1);
+ PG_RETURN_POINTER(result);
+}
+
+/**
+ * ST_ShortestLine(geography, geography)
+ * Return the shortest line between the first and second arguments.
+ */
+PG_FUNCTION_INFO_V1(geography_shortestline);
+Datum geography_shortestline(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED* g1 = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED* g2 = PG_GETARG_GSERIALIZED_P(1);
+ bool use_spheroid = PG_GETARG_BOOL(2);
+ LWGEOM *line, *lwgeom1, *lwgeom2;
+ GSERIALIZED* result;
+ SPHEROID s;
+
+ gserialized_error_if_srid_mismatch(g1, g2, __func__);
+
+ lwgeom1 = lwgeom_from_gserialized(g1);
+ lwgeom2 = lwgeom_from_gserialized(g2);
+
+ /* Return NULL on empty/bad arguments. */
+ if ( !lwgeom1 || !lwgeom2 || lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
+ {
+ PG_FREE_IF_COPY(g1, 0);
+ PG_FREE_IF_COPY(g2, 1);
+ PG_RETURN_NULL();
+ }
+
+ /* Initialize spheroid */
+ // spheroid_init_from_srid(fcinfo, srid, &s);
+ spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS);
+
+ /* Set to sphere if requested */
+ if ( ! use_spheroid )
+ s.a = s.b = s.radius;
+
+ line = geography_tree_shortestline(lwgeom1, lwgeom2, FP_TOLERANCE, &s);
+
+ if (lwgeom_is_empty(line))
+ PG_RETURN_NULL();
+
+ result = geography_serialize(line);
+ lwgeom_free(line);
+
+ PG_FREE_IF_COPY(g1, 0);
+ PG_FREE_IF_COPY(g2, 1);
+ PG_RETURN_POINTER(result);
+}
diff --git a/postgis/geography_measurement_trees.h b/postgis/geography_measurement_trees.h
index 501a2c6e1..6446d5d09 100644
--- a/postgis/geography_measurement_trees.h
+++ b/postgis/geography_measurement_trees.h
@@ -27,14 +27,22 @@
#include "lwgeom_cache.h"
int geography_dwithin_cache(FunctionCallInfo fcinfo,
- SHARED_GSERIALIZED *g1,
- SHARED_GSERIALIZED *g2,
- const SPHEROID *s,
- double tolerance,
- int *dwithin);
+ SHARED_GSERIALIZED *g1,
+ SHARED_GSERIALIZED *g2,
+ const SPHEROID *s,
+ double tolerance,
+ int *dwithin);
int geography_distance_cache(FunctionCallInfo fcinfo,
- SHARED_GSERIALIZED *g1,
- SHARED_GSERIALIZED *g2,
- const SPHEROID *s,
- double *distance);
-int geography_tree_distance(const GSERIALIZED* g1, const GSERIALIZED* g2, const SPHEROID* s, double tolerance, double* distance);
+ SHARED_GSERIALIZED *g1,
+ SHARED_GSERIALIZED *g2,
+ const SPHEROID *s,
+ double *distance);
+
+int geography_tree_distance(
+ const GSERIALIZED* g1,
+ const GSERIALIZED* g2,
+ const SPHEROID* s,
+ double tolerance,
+ double* distance);
+
+
diff --git a/regress/core/geography.sql b/regress/core/geography.sql
index f17bcdfec..1ad7e4e04 100644
--- a/regress/core/geography.sql
+++ b/regress/core/geography.sql
@@ -124,3 +124,31 @@ select 'dwithin_poly_line_1', ST_DWithin('POLYGON((1 1, 2 2, 3 0, 1 1))'::geogra
select 'dwithin_poly_poly_1', ST_DWithin('POLYGON((0 0, -2 -2, -3 0, 0 0))'::geography, 'POLYGON((1 1, 2 2, 3 0, 1 1))'::geography, 10);
select 'dwithin_poly_poly_2', ST_DWithin('POLYGON((0 0, -2 -2, -3 0, 0 0))'::geography, 'POLYGON((1 1, 2 2, 3 0, 1 1))'::geography, 300000);
select 'dwithin_poly_poly_3', ST_DWithin('POLYGON((1 1, -2 -2, -3 0, 1 1))'::geography, 'POLYGON((1 1, 2 2, 3 0, 1 1))'::geography, 300000);
+
+
+-- Linear Referencing functions
+
+SELECT 'lrs_empty_1', ST_LineInterpolatePoint(geography 'Linestring empty', 0.1);
+SELECT 'lrs_empty_2', ST_LineInterpolatePoints(geography 'Linestring empty', 0.1, true);
+SELECT 'lrs_empty_3', ST_LineLocatePoint(geography 'Linestring empty', 'Point(45 45)', true);
+SELECT 'lrs_empty_4', ST_LineSubstring(geography 'Linestring empty', 0.1, 0.2);
+SELECT 'lrs_empty_5', ST_ClosestPoint(geography 'Linestring empty', 'Point(45 45)', true);
+SELECT 'lrs_empty_6', ST_ShortestLine(geography 'Linestring empty', 'Point(45 45)', true);
+
+SELECT 'lrs_lip_1', ST_AsText(ST_LineInterpolatePoint(geography 'Linestring(4.35 50.85, 37.617222 55.755833)', 0.0), 2);
+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_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_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_sl_1', ST_AsText(ST_ShortestLine(geography 'linestring(0 1, 50 1)', 'Point(25 0)', true), 2);
+
diff --git a/regress/core/geography_expected b/regress/core/geography_expected
index 97846a1a0..4fc2e88c5 100644
--- a/regress/core/geography_expected
+++ b/regress/core/geography_expected
@@ -49,3 +49,19 @@ dwithin_poly_line_1|t
dwithin_poly_poly_1|f
dwithin_poly_poly_2|t
dwithin_poly_poly_3|t
+lrs_empty_1|
+lrs_empty_2|
+lrs_empty_3|
+lrs_empty_4|
+lrs_empty_5|
+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_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)
-----------------------------------------------------------------------
Summary of changes:
doc/reference_lrs.xml | 34 +++
doc/reference_measure.xml | 26 +-
liblwgeom/Makefile.in | 1 +
liblwgeom/liblwgeom.h.in | 20 ++
liblwgeom/lwgeodetic_measures.c | 561 ++++++++++++++++++++++++++++++++++
liblwgeom/lwgeodetic_tree.c | 82 ++++-
liblwgeom/lwgeodetic_tree.h | 4 +
postgis/geography.h | 1 +
postgis/geography.sql.in | 72 +++++
postgis/geography_measurement.c | 322 +++++++++++++++++--
postgis/geography_measurement_trees.h | 28 +-
regress/core/geography.sql | 28 ++
regress/core/geography_expected | 16 +
13 files changed, 1149 insertions(+), 46 deletions(-)
create mode 100644 liblwgeom/lwgeodetic_measures.c
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list