2018/xx/xx
* New Features *
- #3876, ST_Angle function (Rémi Cura)
+ - #3564, ST_LineInterpolatePoints (Dan Baston)
PostGIS 2.4.0
2017/09/30
../images/st_linecrossingdirection03.png \
../images/st_linecrossingdirection04.png \
../images/st_line_interpolate_point01.png \
+ ../images/st_line_interpolate_points01.png \
../images/st_line_substring01.png \
../images/st_longestline01.png \
../images/st_longestline02.png \
--- /dev/null
+Style2;LINESTRING(25 50, 100 125, 150 190)
+Style1-thinline;MULTIPOINT(51.5974135047432 76.5974135047432,78.1948270094864 103.194827009486,104.132163186446 130.37181214238,127.066081593223 160.18590607119,150 190)
+
<refsection>
<title>See Also</title>
- <para><xref linkend="ST_AsText" />, <xref linkend="ST_AsEWKT" />, <xref linkend="ST_Length" />, <xref linkend="ST_LineLocatePoint" /></para>
+ <para>
+ <xref linkend="ST_AsText" />,
+ <xref linkend="ST_AsEWKT" />,
+ <xref linkend="ST_Length" />,
+ <xref linkend="ST_LineInterpolatePoints" />
+ <xref linkend="ST_LineLocatePoint" />
+ O</para>
+ </refsection>
+ </refentry>
+
+ <refentry id="ST_LineInterpolatePoints">
+ <refnamediv>
+ <refname>ST_LineInterpolatePoints</refname>
+
+ <refpurpose>
+ Returns one or more points interpolated along a line.
+ </refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_LineInterpolatePoints</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>a_linestring</parameter></paramdef>
+ <paramdef><type>float8 </type> <parameter>a_fraction</parameter></paramdef>
+ <paramdef><type>boolean </type> <parameter>repeat</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>Returns one or more points interpolated along a line. First argument
+ must be a LINESTRING. Second argument is a float8 between 0 and 1
+ representing the spacing between the points as a fraction of total
+ LineString length. If the third argument is false, at most one point
+ will be constructed (the function will be equivalent to <xref linkend="ST_LineInterpolatePoint" />.)
+ </para>
+
+ <para>
+ If the result has zero or one points, it will be returned as a POINT.
+ If it has two or more points, it will be returned as a MULTIPOINT.
+ </para>
+
+
+ <para>Availability: 2.5.0</para>
+ <para>&Z_support;</para>
+ <para>&M_support;</para>
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+ <informalfigure>
+ <mediaobject>
+ <imageobject>
+ <imagedata fileref="images/st_line_interpolate_points01.png" />
+ </imageobject>
+ <caption><para>A linestring with the interpolated points every 20% </para></caption>
+ </mediaobject>
+ </informalfigure>
+ <programlisting>--Return points each 20% along a 2D line
+SELECT ST_AsText(ST_LineInterpolatePoints('LINESTRING(25 50, 100 125, 150 190)', 0.20))
+ st_astext
+----------------
+ MULTIPOINT(51.5974135047432 76.5974135047432,78.1948270094864 103.194827009486,104.132163186446 130.37181214238,127.066081593223 160.18590607119,150 190)
+</programlisting>
+ </refsection>
+
+ <refsection>
+ <title>See Also</title>
+
+ <para>
+ <xref linkend="ST_LineInterpolatePoint" />
+ <xref linkend="ST_LineLocatePoint" />
+ </para>
</refsection>
</refentry>
}
+static void test_lwline_interpolate_points(void)
+{
+ LWLINE* line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING ZM (0 0 3 1, 1 1 2 2, 10 10 4 3)", LW_PARSER_CHECK_NONE));
+ LWLINE* line2d = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING (1 1, 3 7, 9 12)", LW_PARSER_CHECK_NONE));
+ LWLINE* empty_line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING EMPTY", LW_PARSER_CHECK_NONE));
+
+ POINTARRAY* rpa;
+ POINT4D pta;
+ POINT4D ptb;
+ double eps = 1e-10;
+
+ /* Empty line gives empty point */
+ rpa = lwline_interpolate_points(empty_line, 0.5, LW_TRUE);
+ ASSERT_INT_EQUAL(rpa->npoints, 0);
+ ptarray_free(rpa);
+
+ /* Get first endpoint when fraction = 0 */
+ rpa = lwline_interpolate_points(line, 0, LW_TRUE);
+ ASSERT_INT_EQUAL(rpa->npoints, 1);
+ pta = getPoint4d(line->points, 0);
+ ptb = getPoint4d(rpa, 0);
+ ASSERT_POINT4D_EQUAL(pta, ptb, eps);
+ ptarray_free(rpa);
+
+ /* Get last endpoint when fraction = 0 */
+ rpa = lwline_interpolate_points(line, 1, LW_TRUE);
+ ASSERT_INT_EQUAL(rpa->npoints, 1);
+ pta = getPoint4d(line->points, line->points->npoints - 1);
+ ptb = getPoint4d(rpa, 0);
+ ASSERT_POINT4D_EQUAL(pta, ptb, eps);
+ ptarray_free(rpa);
+
+ /* Interpolate a single point */
+ /* First vertex is at 10% */
+ rpa = lwline_interpolate_points(line, 0.1, LW_FALSE);
+ ASSERT_INT_EQUAL(rpa->npoints, 1);
+ pta = getPoint4d(line->points, 1);
+ ptb = getPoint4d(rpa, 0);
+ ASSERT_POINT4D_EQUAL(pta, ptb, eps);
+ ptarray_free(rpa);
+
+ /* 5% is halfway to first vertex */
+ rpa = lwline_interpolate_points(line, 0.05, LW_FALSE);
+ ASSERT_INT_EQUAL(rpa->npoints, 1);
+ pta.x = 0.5;
+ pta.y = 0.5;
+ pta.m = 1.5;
+ pta.z = 2.5;
+ ptb = getPoint4d(rpa, 0);
+ ASSERT_POINT4D_EQUAL(pta, ptb, eps);
+ ptarray_free(rpa);
+
+ /* Now repeat points */
+ rpa = lwline_interpolate_points(line, 0.4, LW_TRUE);
+ ASSERT_INT_EQUAL(rpa->npoints, 2);
+ pta.x = 4;
+ pta.y = 4;
+ ptb = getPoint4d(rpa, 0);
+ ASSERT_POINT2D_EQUAL(pta, ptb, eps);
+
+ pta.x = 8;
+ pta.y = 8;
+ ptb = getPoint4d(rpa, 1);
+ ASSERT_POINT2D_EQUAL(pta, ptb, eps);
+ ptarray_free(rpa);
+
+ /* Make sure it works on 2D lines */
+ rpa = lwline_interpolate_points(line2d, 0.4, LW_TRUE);
+ ASSERT_INT_EQUAL(rpa->npoints, 2);
+ CU_ASSERT_FALSE(ptarray_has_z(rpa));
+ CU_ASSERT_FALSE(ptarray_has_m(rpa));
+ ptarray_free(rpa);
+
+ lwgeom_free(lwline_as_lwgeom(line));
+ lwgeom_free(lwline_as_lwgeom(line2d));
+ lwgeom_free(lwline_as_lwgeom(empty_line));
+}
+
static void test_lwline_clip(void)
{
LWCOLLECTION *c;
PG_ADD_TEST(suite,test_lwpoint_set_ordinate);
PG_ADD_TEST(suite,test_lwpoint_get_ordinate);
PG_ADD_TEST(suite,test_point_interpolate);
+ PG_ADD_TEST(suite,test_lwline_interpolate_points);
PG_ADD_TEST(suite,test_lwline_clip);
PG_ADD_TEST(suite,test_lwline_clip_big);
PG_ADD_TEST(suite,test_lwmline_clip);
CU_PASS(); \
} while(0);
+#define ASSERT_POINT2D_EQUAL(o, e, eps) do { \
+ CU_ASSERT_DOUBLE_EQUAL(o.x, e.x, eps); \
+ CU_ASSERT_DOUBLE_EQUAL(o.y, e.y, eps); \
+} while(0);
+
+#define ASSERT_POINT4D_EQUAL(o, e, eps) do { \
+ CU_ASSERT_DOUBLE_EQUAL(o.x, e.x, eps); \
+ CU_ASSERT_DOUBLE_EQUAL(o.y, e.y, eps); \
+ CU_ASSERT_DOUBLE_EQUAL(o.z, e.z, eps); \
+ CU_ASSERT_DOUBLE_EQUAL(o.m, e.m, eps); \
+} while(0);
+
/* Utility functions */
void do_fn_test(LWGEOM* (*transfn)(LWGEOM*), char *input_wkt, char *expected_wkt);
*/
extern int lwline_add_lwpoint(LWLINE *line, LWPOINT *point, int where);
+/**
+ * Interpolate one or more points along a line
+ */
+extern POINTARRAY* lwline_interpolate_points(const LWLINE *line, double length_fraction, char repeat);
+
/******************************************************************
* LWPOLY functions
******************************************************************/
extern char* lwpoint_to_latlon(const LWPOINT *p, const char *format);
extern int lwgeom_startpoint(const LWGEOM* lwgeom, POINT4D* pt);
-extern void interpolate_point4d(POINT4D *A, POINT4D *B, POINT4D *I, double F);
+extern void interpolate_point4d(const POINT4D *A, const POINT4D *B, POINT4D *I, double F);
/**
* Ensure the outer ring is clockwise oriented and all inner rings
* F=.2 : A-I-------B
*/
void
-interpolate_point4d(POINT4D *A, POINT4D *B, POINT4D *I, double F)
+interpolate_point4d(const POINT4D *A, const POINT4D *B, POINT4D *I, double F)
{
#if PARANOIA_LEVEL > 0
double absF=fabs(F);
}
+POINTARRAY* lwline_interpolate_points(const LWLINE *line, double length_fraction, 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;
+
+ /* Empty.InterpolatePoint == Point Empty */
+ if ( lwline_is_empty(line) )
+ {
+ return ptarray_construct_empty(has_z, has_m, 0);
+ }
+
+ /* 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 opa;
+ }
+
+ /* Interpolate points along the line */
+ length = ptarray_length_2d(ipa);
+ points_to_interpolate = repeat ? (uint32_t) floor(1 / length_fraction) : 1;
+ opa = ptarray_construct(has_z, has_m, points_to_interpolate);
+
+ const POINT2D* p1 = getPoint2d_cp(ipa, 0);
+ for ( i = 0; i < ipa->npoints - 1 && points_found < points_to_interpolate; i++ )
+ {
+ const POINT2D* p2 = getPoint2d_cp(ipa, i+1);
+ double segment_length_frac = distance2d_pt_pt(p1, p2) / 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 )
+ {
+ POINT4D p1_4d = getPoint4d(ipa, i);
+ POINT4D p2_4d = getPoint4d(ipa, i+1);
+
+ double segment_fraction = (length_fraction - length_fraction_consumed) / segment_length_frac;
+ interpolate_point4d(&p1_4d, &p2_4d, &pt, segment_fraction);
+ ptarray_set_point4d(opa, points_found++, &pt);
+ length_fraction += length_fraction_increment;
+ }
+
+ length_fraction_consumed += segment_length_frac;
+
+ p1 = p2;
+ }
+
+ /* 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);
+ }
+
+ return opa;
+}
+
#include "access/htup_details.h"
-/***********************************************************************
- * Simple Douglas-Peucker line simplification.
- * No checks are done to avoid introduction of self-intersections.
- * No topology relations are considered.
- *
- * --strk@kbt.io;
- ***********************************************************************/
-
-
/* Prototypes */
Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS);
Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS);
+Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS);
Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS);
Datum ST_MinimumBoundingRadius(PG_FUNCTION_ARGS);
Datum ST_MinimumBoundingCircle(PG_FUNCTION_ARGS);
static int point_in_ring(POINTARRAY *pts, const POINT2D *point);
static int point_in_ring_rtree(RTREE_NODE *root, const POINT2D *point);
+/***********************************************************************
+ * Simple Douglas-Peucker line simplification.
+ * No checks are done to avoid introduction of self-intersections.
+ * No topology relations are considered.
+ *
+ * --strk@kbt.io;
+ ***********************************************************************/
PG_FUNCTION_INFO_V1(LWGEOM_simplify2d);
Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS)
* Interpolate a point along a line, useful for Geocoding applications
* SELECT line_interpolate_point( 'LINESTRING( 0 0, 2 2'::geometry, .5 )
* returns POINT( 1 1 ).
- * Works in 2d space only.
- *
- * Initially written by: jsunday@rochgrp.com
- * Ported to LWGEOM by: strk@refractions.net
***********************************************************************/
-
-Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS);
-
PG_FUNCTION_INFO_V1(LWGEOM_line_interpolate_point);
Datum LWGEOM_line_interpolate_point(PG_FUNCTION_ARGS)
{
GSERIALIZED *gser = PG_GETARG_GSERIALIZED_P(0);
GSERIALIZED *result;
- double distance = PG_GETARG_FLOAT8(1);
- LWLINE *line;
- LWGEOM *geom;
- LWPOINT *point;
- POINTARRAY *ipa, *opa;
- POINT4D pt;
- int nsegs, i;
- double length, slength, tlength;
-
- if ( distance < 0 || distance > 1 )
+ double distance_fraction = PG_GETARG_FLOAT8(1);
+ int repeat = PG_NARGS() > 2 && PG_GETARG_BOOL(2);
+ int srid = gserialized_get_srid(gser);
+ LWLINE* lwline;
+ LWGEOM* lwresult;
+ POINTARRAY* opa;
+
+ if ( distance_fraction < 0 || distance_fraction > 1 )
{
elog(ERROR,"line_interpolate_point: 2nd arg isn't within [0,1]");
+ PG_FREE_IF_COPY(gser, 0);
PG_RETURN_NULL();
}
if ( gserialized_get_type(gser) != LINETYPE )
{
elog(ERROR,"line_interpolate_point: 1st arg isn't a line");
+ PG_FREE_IF_COPY(gser, 0);
PG_RETURN_NULL();
}
- /* Empty.InterpolatePoint == Point Empty */
- if ( gserialized_is_empty(gser) )
- {
- point = lwpoint_construct_empty(gserialized_get_srid(gser), gserialized_has_z(gser), gserialized_has_m(gser));
- result = geometry_serialize(lwpoint_as_lwgeom(point));
- lwpoint_free(point);
- PG_RETURN_POINTER(result);
- }
+ lwline = lwgeom_as_lwline(lwgeom_from_gserialized(gser));
+ opa = lwline_interpolate_points(lwline, distance_fraction, repeat);
- geom = lwgeom_from_gserialized(gser);
- line = lwgeom_as_lwline(geom);
- ipa = line->points;
+ lwgeom_free(lwline_as_lwgeom(lwline));
+ PG_FREE_IF_COPY(gser, 0);
- /* If distance is one of the two extremes, return the point on that
- * end rather than doing any expensive computations
- */
- if ( distance == 0.0 || distance == 1.0 )
+ if (opa->npoints <= 1)
{
- if ( distance == 0.0 )
- getPoint4d_p(ipa, 0, &pt);
- else
- getPoint4d_p(ipa, ipa->npoints-1, &pt);
-
- opa = ptarray_construct(lwgeom_has_z(geom), lwgeom_has_m(geom), 1);
- ptarray_set_point4d(opa, 0, &pt);
-
- point = lwpoint_construct(line->srid, NULL, opa);
- PG_RETURN_POINTER(geometry_serialize(lwpoint_as_lwgeom(point)));
+ lwresult = lwpoint_as_lwgeom(lwpoint_construct(srid, NULL, opa));
+ } else {
+ lwresult = lwmpoint_as_lwgeom(lwmpoint_construct(srid, opa));
}
- /* Interpolate a point on the line */
- nsegs = ipa->npoints - 1;
- length = ptarray_length_2d(ipa);
- tlength = 0;
- for ( i = 0; i < nsegs; i++ )
- {
- POINT4D p1, p2;
- POINT4D *p1ptr=&p1, *p2ptr=&p2; /* don't break
- * strict-aliasing rules
- */
+ result = geometry_serialize(lwresult);
+ lwgeom_free(lwresult);
- getPoint4d_p(ipa, i, &p1);
- getPoint4d_p(ipa, i+1, &p2);
-
- /* Find the relative length of this segment */
- slength = distance2d_pt_pt((POINT2D*)p1ptr, (POINT2D*)p2ptr)/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.
- */
- if ( distance < tlength + slength )
- {
- double dseg = (distance - tlength) / slength;
- interpolate_point4d(&p1, &p2, &pt, dseg);
- opa = ptarray_construct(lwgeom_has_z(geom), lwgeom_has_m(geom), 1);
- ptarray_set_point4d(opa, 0, &pt);
- point = lwpoint_construct(line->srid, NULL, opa);
- PG_RETURN_POINTER(geometry_serialize(lwpoint_as_lwgeom(point)));
- }
- tlength += slength;
- }
-
- /* Return the last point on the line. This shouldn't happen, but
- * could if there's some floating point rounding errors. */
- getPoint4d_p(ipa, ipa->npoints-1, &pt);
- opa = ptarray_construct(lwgeom_has_z(geom), lwgeom_has_m(geom), 1);
- ptarray_set_point4d(opa, 0, &pt);
- point = lwpoint_construct(line->srid, NULL, opa);
- PG_FREE_IF_COPY(gser, 0);
- PG_RETURN_POINTER(geometry_serialize(lwpoint_as_lwgeom(point)));
+ PG_RETURN_POINTER(result);
}
/***********************************************************************
* --jsunday@rochgrp.com;
AS 'MODULE_PATHNAME', 'LWGEOM_line_interpolate_point'
LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL;
+-- Availability: 2.5.0
+CREATE OR REPLACE FUNCTION ST_LineInterpolatePoints(geometry, float8, repeat boolean DEFAULT true)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME', 'LWGEOM_line_interpolate_point'
+ LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL;
+
-- Availability: 1.2.2
-- Deprecation in 2.1.0
CREATE OR REPLACE FUNCTION ST_line_interpolate_point(geometry, float8)
select 'line_interpolate_point', ST_AsText(ST_LineInterpolatePoint('LINESTRING(0 0, 1 1)', 0));
select 'line_interpolate_point', ST_AsText(ST_LineInterpolatePoint('LINESTRING(0 0 10, 1 1 5)', 0.5));
+--
+--- ST_LineInterpolatePoints
+--
+
+select 'line_interpolate_points', ST_AsText(ST_LineInterpolatePoints('LINESTRING(0 0, 1 1)', 0.7));
+select 'line_interpolate_points', ST_AsText(ST_LineInterpolatePoints('LINESTRING(0 0, 1 1)', 0.3));
--
-- ST_AddMeasure
line_substring_12|POINT Z (0.5 0.5 7.5)
line_interpolate_point|POINT(0 0)
line_interpolate_point|POINT Z (0.5 0.5 7.5)
+line_interpolate_points|POINT(0.7 0.7)
+line_interpolate_points|MULTIPOINT(0.3 0.3,0.6 0.6,0.9 0.9)
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
projects / postgis / commitdiff