From: Sandro Santilli Date: Thu, 9 Jun 2005 12:24:39 +0000 (+0000) Subject: NEW line_substring() function X-Git-Tag: pgis_1_1_0~347 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=744d8ecf04fa000af815ddabcc3f898cb55e3478;p=postgis NEW line_substring() function git-svn-id: http://svn.osgeo.org/postgis/trunk@1745 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/CHANGES b/CHANGES index 19726ca7f..cd6bc036b 100644 --- a/CHANGES +++ b/CHANGES @@ -9,6 +9,7 @@ PostGIS 1.1.0CVS - full autoconf-based configuration - added scale() and transscale() companion methods to translate() - initial implementation of postgis_proc_upgrade script + - NEW line_substring() function PostGIS 1.0.1 2005/05/24 diff --git a/TODO b/TODO index 8069080c6..e4141595b 100644 --- a/TODO +++ b/TODO @@ -41,12 +41,6 @@ Tue May 24 18:55:46 CEST 2005 Returns a float, representing the proportional distance the point is located along the line. -- line_interpolate_line(geometry, proportion, direction) - - Returns a line, representing the subset of the input line - that is the proportional distance along the line. Direction - determines what end to measure from. (??) - - MakeValidShape() to remove consecutive vertexes and any eventual other cleanup aimed at making pgsql2shp produce "valid" Shapefiles. diff --git a/doc/postgis.xml b/doc/postgis.xml index 106737340..1518a8e93 100644 --- a/doc/postgis.xml +++ b/doc/postgis.xml @@ -4415,22 +4415,28 @@ FROM geometry_table; - scale(geometry,float8,float8,float8) + Scale(geometry,float8,float8,float8) scales the geometry to a new size by multiplying the ordinates with the parameters. Ie: scale(geom, Xfactor, Yfactor, Zfactor). + + Availability: 1.1.0 + - transscale(geometry,float8,float8,float8,float8) + TransScale(geometry,float8,float8,float8,float8) First, translates the geometry using the first two floats, then scales it using the second two floats, working in 2D only. Using transscale(geom, X, Y, XFactor, YFactor) effectively is an efficient shortcut for scale(translate(geom,X,Y,0),Xfactor,YFactor,1). + + Availability: 1.1.0 + @@ -4496,6 +4502,41 @@ FROM geometry_table; + + Linear Referencing + + + + + line_interpolate_point(geometry, proportion) + + + Interpolates a point along a line. First argument must be a + LINESTRING. Second argument is a float between 0 and 1. Returns a + point. + + + + + + line_substring(geometry, start, end) + + + + Return a linestring being a substring of the input one starting + and ending at the given fractions of total length. Second + and third arguments are float8 values between 0 and 1. + + + Availability: 1.1.0 + + + + + + + + Misc @@ -4652,16 +4693,6 @@ FROM geometry_table; - - line_interpolate_point(geometry, proportion) - - - Interpolates a point along a line. First argument must be a - LINESTRING. Second argument is a float between 0 and 1. Returns a - point. - - - Accum(geometry set) @@ -4675,6 +4706,7 @@ FROM geometry_table; + diff --git a/lwgeom/liblwgeom.h b/lwgeom/liblwgeom.h index 4a3e7d13f..e38c01431 100644 --- a/lwgeom/liblwgeom.h +++ b/lwgeom/liblwgeom.h @@ -1005,6 +1005,7 @@ extern int ptarray_isclosed2d(const POINTARRAY *pa); extern int32 lwgeom_nrings_recursive(uchar *serialized); extern void ptarray_reverse(POINTARRAY *pa); +extern POINTARRAY *ptarray_substring(POINTARRAY *, double, double); // Ensure every segment is at most 'dist' long. // Returned LWGEOM might is unchanged if a POINT. diff --git a/lwgeom/lwgeom_functions_analytic.c b/lwgeom/lwgeom_functions_analytic.c index ee7e99404..1ac2b0abe 100644 --- a/lwgeom/lwgeom_functions_analytic.c +++ b/lwgeom/lwgeom_functions_analytic.c @@ -853,3 +853,43 @@ Datum LWGEOM_snaptogrid(PG_FUNCTION_ARGS) /*********************************************************************** * --strk@keybit.net ***********************************************************************/ + +Datum LWGEOM_line_substring(PG_FUNCTION_ARGS); + +PG_FUNCTION_INFO_V1(LWGEOM_line_substring); +Datum LWGEOM_line_substring(PG_FUNCTION_ARGS) +{ + PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + double from = PG_GETARG_FLOAT8(1); + double to = PG_GETARG_FLOAT8(2); + LWLINE *iline; + LWLINE *oline; + POINTARRAY *ipa, *opa; + PG_LWGEOM *ret; + + if( from < 0 || from > 1 ) { + elog(ERROR,"line_interpolate_point: 2nd arg isnt within [0,1]"); + PG_RETURN_NULL(); + } + + if( to < 0 || to > 1 ) { + elog(ERROR,"line_interpolate_point: 3rd arg isnt within [0,1]"); + PG_RETURN_NULL(); + } + + if( lwgeom_getType(geom->type) != LINETYPE ) { + elog(ERROR,"line_interpolate_point: 1st arg isnt a line"); + PG_RETURN_NULL(); + } + + iline = lwline_deserialize(SERIALIZED_FORM(geom)); + ipa = iline->points; + + opa = ptarray_substring(ipa, from, to); + + oline = lwline_construct(iline->SRID, 0, opa); + ret = pglwgeom_serialize((LWGEOM *)oline); + PG_FREE_IF_COPY(geom, 0); + lwgeom_release((LWGEOM *)oline); + PG_RETURN_POINTER(ret); +} diff --git a/lwgeom/lwpostgis.sql.in b/lwgeom/lwpostgis.sql.in index d0d8ac38b..ee186a266 100644 --- a/lwgeom/lwpostgis.sql.in +++ b/lwgeom/lwpostgis.sql.in @@ -3318,6 +3318,11 @@ CREATEFUNCTION line_interpolate_point(geometry, float8) AS '@MODULE_FILENAME@', 'LWGEOM_line_interpolate_point' LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable); +CREATEFUNCTION line_substring(geometry, float8, float8) + RETURNS geometry + AS '@MODULE_FILENAME@', 'LWGEOM_line_substring' + LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable); + CREATEFUNCTION segmentize(geometry, float8) RETURNS geometry AS '@MODULE_FILENAME@', 'LWGEOM_segmentize2d' diff --git a/lwgeom/ptarray.c b/lwgeom/ptarray.c index 4e4d8fc09..95bdf8930 100644 --- a/lwgeom/ptarray.c +++ b/lwgeom/ptarray.c @@ -479,3 +479,92 @@ ptarray_compute_box3d_p(const POINTARRAY *pa, BOX3D *result) return 1; } + +POINTARRAY * +ptarray_substring(POINTARRAY *ipa, double from, double to) +{ + POINTARRAY *opa; + int nopts=0; + uchar *opts; + uchar *optsptr; + POINT4D pt; + POINT2D p1, p2; + int nsegs, i; + int ptsize; + double length, slength, tlength; + int state = 0; // 0=before, 1=inside, 2=after + + /* Allocate memory for a full copy of input points */ + ptsize=TYPE_NDIMS(ipa->dims)*sizeof(double); + opts = lwalloc(ptsize*ipa->npoints); + optsptr = opts; + + /* Interpolate a point on the line */ + nsegs = ipa->npoints - 1; + length = lwgeom_pointarray_length2d(ipa); + + //lwnotice("Total length: %g", length); + + from = length*from; + to = length*to; + + //lwnotice("From/To: %g/%g", from, to); + + tlength = 0; + getPoint2d_p(ipa, 0, &p1); + for( i = 0; i < nsegs; i++ ) { + + getPoint2d_p(ipa, i+1, &p2); + + /* Find the relative length of this segment */ + slength = distance2d_pt_pt(&p1, &p2); + + /* + * 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 ( state == 0 ) + { + if ( from < tlength + slength ) { + double dseg = (from - tlength) / slength; + pt.x = (p1.x) + ((p2.x - p1.x) * dseg); + pt.y = (p1.y) + ((p2.y - p1.y) * dseg); + pt.z = 0; pt.m = 0; + memcpy(optsptr, &pt, ptsize); + optsptr+=ptsize; + nopts++; state++; + } + } + + if ( state == 1 ) + { + if ( to < tlength + slength ) { + double dseg = (to - tlength) / slength; + pt.x = (p1.x) + ((p2.x - p1.x) * dseg); + pt.y = (p1.y) + ((p2.y - p1.y) * dseg); + pt.z = 0; pt.m = 0; + memcpy(optsptr, &pt, ptsize); + optsptr+=ptsize; + nopts++; break; + } else { + memcpy(optsptr, &p2, ptsize); + optsptr+=ptsize; + nopts++; + } + } + + + tlength += slength; + memcpy(&p1, &p2, sizeof(POINT2D)); + } + + //lwnotice("Out of loop, constructing ptarray with %d points", nopts); + + opa = pointArray_construct(opts, + TYPE_HASZ(ipa->dims), + TYPE_HASM(ipa->dims), + nopts); + + return opa; +}