Go to content Go to navigation and search

Home

Current PostGIS Blog Articles


Search

Browse

RSS / Atom

Email me

textpattern

Creative Commons License
All Blog Articles, Data Models and Free Source Code by Simon Greener, The SpatialDB Advisor is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.

ST_Parallel for PostGIS

Thursday May 10 2012 at 15:23

KeywordsST_OffsetCurve parallel postgis
Summary

I present a pgPlSql based function that, for a (m)linestring creates an offset linestring left or right of the original. This function is predominantly for versions of PostGis before 2.0 as 2.0 contains ST_OffsetCurve.

Many people have asked me to convert my Oracle Parallel function to PostGIS. I did it a year or so ago but the results were not as I wished. I have now corrected much that was wrong with the function so I hope it is of more value to those who have asked for it in the past.

Why is it that, while many wanted the function, no one ever offered to help getting it to work? Isn’t collaboration part and parcel of how “open source” works?

Anyway, here is the function with some tests. I have performed reasonable checks by throwing different geometries and parameter settings at it but it may still break as the construction of the output WKT is not that systematic.

Some of the required support functions are not in this article but are in the referenced source code file at the bottom of the article.

Finally, this function, and the Oracle original, do not handle the situation where two vectors, in the acute case, at certain parallel distances, decompose to nothing.

  1. -- Drop function
  2. DROP FUNCTION IF EXISTS st_parallel(geometry, NUMERIC, INTEGER, INTEGER);
  3. -- Create Function
  4. CREATE OR REPLACE FUNCTION st_parallel(p_geometry     geometry,
  5.                                        p_distance     NUMERIC,
  6.                                        p_round_factor INTEGER,
  7.                                        p_curved       INTEGER DEFAULT 0)
  8.   RETURNS geometry AS
  9. $BODY$
  10. DECLARE
  11.     v_rowcount        INTEGER;
  12.     v_part_no         INTEGER;
  13.     v_result          geometry;
  14.     bAcute            BOOLEAN := FALSE;
  15.     bCurved           BOOLEAN := ( 1 = p_curved );
  16.     v_dims            INTEGER;
  17.     v_delta           CoordType;
  18.     v_prev_delta      CoordType;
  19.     v_adj_Coord       CoordType;
  20.     v_int_1           CoordType;
  21.     v_int_coord       CoordType;
  22.     v_int_2           CoordType;
  23.     v_startCoord      CoordType;
  24.     v_endCoord        CoordType;
  25.     v_distance        NUMERIC;
  26.     v_ratio           NUMERIC;
  27.     v_az              NUMERIC;
  28.     v_dir             INTEGER;
  29.     v_return_geom     geometry;
  30.     v_geom_string     text;
  31.   c_parts CURSOR ( p_geom         geometry,
  32.                    p_Geometrytype text )
  33.   IS
  34.     SELECT p_geom AS geom
  35.      WHERE p_geometrytype IN ('ST_LineString' )
  36.     UNION ALL
  37.     SELECT ST_GeometryN(p_geom,generate_series(1,ST_NumGeometries(p_geom))) AS geom
  38.      WHERE p_GeometryType = 'ST_MultiLineString';
  39.   c_vector CURSOR ( p_linestring geometry )
  40.   IS
  41.     SELECT startCoord, endCoord
  42.       FROM (SELECT (ST_GetVector(p_linestring)).* ) o;
  43. BEGIN
  44.   /* Problem with this algorithm is that any existing circular curved elements will not be honoured unless bCurved is 0 */
  45.   IF ( p_round_factor IS NULL ) THEN
  46.      raise exception 'p_round_factor may not be null';
  47.   END IF;
  48.   IF ( p_geometry IS NULL
  49.        OR
  50.        ST_GeometryType(p_geometry) NOT IN ('ST_MultiLineString','ST_LineString') ) THEN
  51.      Raise Exception 'p_linestring is null or is not a linestring or multilinestring';
  52.   END IF;
  53.   IF ST_HasArc(p_geometry) THEN /* Should never happen as GeometryType captures anything with a curve */
  54.      Raise Exception 'Compound linestrings are not supported.';
  55.   END IF;
  56.   /* Process all parts of a, potentially, multi-part linestring geometry */
  57.    v_part_no := 0;
  58.    FOR part IN c_parts(p_geometry,ST_GeometryType(p_geometry)) loop
  59.     BEGIN
  60.       v_part_no := v_part_no + 1;
  61.       v_dims := ST_NDims(part.geom);
  62.       IF ( v_part_no > 1 ) THEN
  63.          v_geom_string := v_geom_string || CASE WHEN bCurved THEN ',' ELSE ',(' END;
  64.       END IF;      
  65.       /* Process geometry one vector at a time */
  66.       v_rowcount := 0;
  67.       FOR rec IN c_vector(part.geom) LOOP
  68.         v_rowcount := v_rowcount + 1;
  69.         /* Compute base offset */
  70.         v_startCoord := rec.startCoord;
  71.         v_endCoord   := rec.endCoord;
  72.         v_az := ST_Azimuth(ST_MakePoint(v_startCoord.X,v_startCoord.Y),
  73.                            ST_MakePoint(v_endCoord.X,  v_endCoord.Y  ) );
  74.         v_dir     := CASE WHEN v_az < pi() THEN -1 ELSE 1 END;
  75.         v_delta.x := ABS(COS(v_az)) * p_distance * v_dir;
  76.         v_delta.y := ABS(SIN(v_az)) * p_distance * v_dir;
  77.         IF  NOT ( v_az > pi()/2 AND
  78.                   v_az < pi() OR
  79.                   v_az > 3 * pi()/2 ) THEN
  80.           v_delta.x := -1 * v_delta.x;
  81.         END IF;
  82.         /* merge vectors at this point? */
  83.         IF (v_ROWCOUNT > 1) THEN
  84.            /* Get intersection of two lines parallel at distance p_distance from current ones */
  85.            v_int_coord := rec.startCoord;
  86.            SELECT *
  87.              FROM findlineintersection(v_adj_coord.x,                   v_adj_coord.y,
  88.                                        v_startCoord.x + v_prev_delta.x, v_startCoord.y + v_prev_delta.y,
  89.                                        v_endCoord.x   + v_delta.x,      v_endCoord.y   + v_delta.y,
  90.                                        v_startCoord.x + v_delta.x,      v_startCoord.y + v_delta.y
  91.              INTO v_int_coord.x, v_int_coord.y,
  92.                   v_int_1.x,     v_int_1.y,
  93.                   v_int_2.x,     v_int_2.y);
  94.            IF ( v_int_coord.x = 1E+308 ) THEN
  95.              /* No intersection point as lines are parallel */
  96.              bAcute := TRUE;
  97.              /* int coord could be computed from start or end, doesn't matter. */
  98.              v_int_coord.x := Round((v_startCoord.x + v_prev_delta.x)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
  99.              v_int_coord.y := Round((v_startCoord.y + v_prev_delta.y)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
  100.              v_int_coord.z := v_startCoord.z::DOUBLE PRECISION;
  101.              v_int_coord.m := v_startCoord.m::DOUBLE PRECISION;
  102.            ELSE
  103.              bAcute := ( Round(v_int_coord.x::NUMERIC,p_round_factor::INTEGER) = Round(v_int_1.x::NUMERIC,p_round_factor::INTEGER) )
  104.                    AND ( Round(v_int_coord.x::NUMERIC,p_round_factor::INTEGER) = Round(v_int_2.x::NUMERIC,p_round_factor::INTEGER) )
  105.                    AND ( Round(v_int_coord.y::NUMERIC,p_round_factor::INTEGER) = Round(v_int_1.y::NUMERIC,p_round_factor::INTEGER) )
  106.                    AND ( Round(v_int_coord.y::NUMERIC,p_round_factor::INTEGER) = Round(v_int_2.y::NUMERIC,p_round_factor::INTEGER) );
  107.            END IF;
  108.            IF ( bCurved AND NOT bAcute) THEN
  109.              /* First point in "intersection" curve */
  110.              v_int_1.x := Round((v_startCoord.x + v_prev_delta.x)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
  111.              v_int_1.y := Round((v_startCoord.y + v_prev_delta.y)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
  112.              v_int_1.z := v_startCoord.z;
  113.              v_int_1.m := v_startCoord.m;
  114.              /* Intersection point (top of circular arc) */
  115.              /* Need to compute coordinates at mid-angle of the circular arc formed by last and current vector */
  116.              v_distance := ST_Distance(ST_SetSRID(ST_Point( v_int_coord.x, v_int_coord.y),ST_SRID(p_geometry)),
  117.                                        ST_SetSRID(ST_Point(v_startCoord.x,v_startCoord.y),ST_SRID(p_geometry)));
  118.              v_ratio := ( p_distance / v_distance ) * SIGN(p_distance);
  119.              v_adj_coord   := rec.startCoord;
  120.              v_adj_coord.x := Round((v_startCoord.x + (( v_int_coord.x - v_startCoord.x ) * v_ratio ))::NUMERIC,p_round_factor::INTEGER);
  121.              v_adj_coord.y := Round((v_startCoord.y + (( v_int_coord.y - v_startCoord.y ) * v_ratio ))::NUMERIC,p_round_factor::INTEGER);
  122.              /* Last point in "intersection" curve */
  123.              v_int_2.x := Round((v_startCoord.x + v_delta.x)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
  124.              v_int_2.y := Round((v_startCoord.y + v_delta.y)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
  125.              v_int_2.z := v_startCoord.z::DOUBLE PRECISION;
  126.              v_int_2.m := v_startCoord.m::DOUBLE PRECISION;
  127.            ELSE
  128.              /* Intersection point */
  129.              v_adj_coord   := v_int_coord;
  130.              v_adj_coord.x := Round(v_int_coord.x::NUMERIC,p_round_factor::INTEGER);
  131.              v_adj_coord.y := Round(v_int_coord.y::NUMERIC,p_round_factor::INTEGER);
  132.            END IF;
  133.         ELSE  /* c_vector%ROWCOUNT = 1 */
  134.            v_geom_string := /*GeometryType(p_geometry) || */ '(';
  135.           /* Translate start point with current vector */
  136.           v_adj_coord.x := Round((v_startCoord.x + v_delta.x)::NUMERIC, p_round_factor::INTEGER)::DOUBLE PRECISION;
  137.           v_adj_coord.y := Round((v_startCoord.y + v_delta.y)::NUMERIC, p_round_factor::INTEGER)::DOUBLE PRECISION;
  138.           v_adj_coord.z := v_startCoord.z::DOUBLE PRECISION;
  139.           v_adj_coord.m := v_startCoord.m::DOUBLE PRECISION;
  140.         END IF;
  141.         /* Now add computed coordinates to output geometry */
  142.         IF (NOT bCurved) OR bAcute OR (v_ROWCOUNT = 1 ) THEN
  143.       v_geom_string := v_geom_string ||
  144.                        CoordAsText(v_adj_coord,v_dims)::text ||
  145.                        ','::text;
  146.         ElsIf (bCurved) THEN
  147.           /* With any generated circular curves we need to add the appropriate type */
  148.       v_geom_string := v_geom_string ||
  149.                        CoordAsText(v_int_1,v_dims)::text ||
  150.                        '),CIRCULARSTRING('::text || CoordAsText(v_int_1,    v_dims)::text || ','::text ||
  151.                                                     CoordAsText(v_adj_coord,v_dims)::text || ','::text ||
  152.                                                     CoordAsText(v_int_2,    v_dims)::text || ')'::text ||
  153.                        ',(' || CoordAsText(v_int_2,v_dims)::text || ','::text;
  154.         END IF;
  155.         v_prev_delta := v_delta;
  156.       END LOOP;
  157.       /* Append Point at End of Linestring */
  158.       v_geom_string := v_geom_string ||
  159.              CoordAsText(create_coordtype(Round((v_endcoord.x + v_delta.x)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION,
  160.                                           Round((v_endcoord.y + v_delta.y)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION,
  161.                                           v_endcoord.z::DOUBLE PRECISION,
  162.                                           v_endcoord.m::DOUBLE PRECISION),v_dims)::text || ')';
  163.     END;
  164.   END Loop;
  165.   IF (bCurved AND v_geom_string LIKE '%CIRCULARSTRING%' ) THEN
  166.       v_geom_string := CASE WHEN ST_GeometryType(p_geometry) = 'ST_MultiLineString'
  167.                             THEN 'MULTICURVE('
  168.                             ELSE 'COMPOUNDCURVE('
  169.                         END || v_geom_String || ')';
  170.    ELSE
  171.       v_geom_string := CASE WHEN ST_GeometryType(p_geometry) = 'ST_MultiLineString'
  172.                             THEN GeometryType(p_geometry) || '(' || v_geom_String || ')'
  173.                             ELSE GeometryType(p_geometry) || v_geom_String
  174.                         END ;
  175.    END IF;                              
  176.   v_result := ST_GeometryFromText(v_geom_string,ST_SRID(p_geometry));
  177.   RETURN v_Result;
  178. END;
  179. $BODY$
  180.   LANGUAGE plpgsql IMMUTABLE STRICT
  181.   COST 100;

Here is a test:

  1. WITH geometries AS (
  2.           SELECT ST_GeomFromText('LINESTRING(1 1,1 10)') AS geom,
  3.                  10.0 AS offset,2 AS roundFactor,0 AS curved
  4. UNION ALL SELECT ST_GeomFromText('LINESTRING(0 0,1 1,1 2)') AS geom,
  5.                   0.5 AS offset,2 AS roundFactor,
  6.                  generate_series(0,1,1) AS curved
  7. UNION ALL SELECT ST_GeomFromText('LINESTRING(0.0 0.0, 45.0 45.0, 90.0 0.0, 135.0 45.0, 180.0 0.0, 180.0 -45.0, 45.0 -45.0, 0.0 0.0)') AS geom,
  8.                  10.0 AS offset,2 AS roundFactor,
  9.                  generate_series(0,1,1) AS curved
  10. UNION ALL SELECT ST_GeomFromText('MULTILINESTRING((0 0,1 1,1 2),(2 3,3 2,5 4))') AS geom,
  11.                   0.5 AS offsetRight,2 AS roundFactor,
  12.                  generate_series(0,1,1) AS curved
  13. )
  14. SELECT ST_AsText(g.geom) AS origGeom, g.offset,g.curved,
  15.        ST_AsText(ST_Parallel(g.geom,g.offset,g.roundFactor,g.curved)) AS geomWKTLeft,
  16.        ST_AsText(ST_Parallel(g.geom,0.0-g.offset,g.roundFactor,g.curved)) AS geomWKTRight
  17.   FROM geometries AS g;

Results ….

OrigGeom Offset Curved GeomWKTLeft geomWKTRight
LINESTRING (1 1,1 10) 10.0 0 LINESTRING (11 1,11 10) LINESTRING (-9 1,-9 10)
LINESTRING (0 0,1 1,1 2) 0.5 0 LINESTRING (0.35 -0.35,1.5 0.79,1.5 2) LINESTRING (-0.35 0.35,0.5 1.21,0.5 2)
LINESTRING (0 0,1 1,1 2) 0.5 1 COMPOUNDCURVE ( (0.35 -0.35,1.35 0.65),CIRCULARSTRING (1.35 0.65,1.46 0.81,1.5 1), (1.5 1,1.5 2)) LINESTRING (-0.35 0.35,0.5 1.21,0.5 2)
LINESTRING (0 0,45 45,90 0,135 45,180 0,180 -45,45 -45,0 0) 10.0 0 LINESTRING (7.07 -7.07,45 30.86,90 -14.14,135 30.86,170 -4.14,170 -35,49.14 -35,7.07 7.07) LINESTRING (-7.07 7.07,45 59.14,90 14.14,135 59.14,190 4.14,190 -55,40.86 -55,-7.07 -7.07)
LINESTRING (0 0,45 45,90 0,135 45,180 0,180 -45,45 -45,0 0) 10.0 1 COMPOUNDCURVE ( (7.07 -7.07,45 30.86,82.93 -7.07),CIRCULARSTRING (82.93 -7.07,90 -10,97.07 -7.07), (97.07 -7.07,134.71 31.15,170 -4.14,170 -35,49.14 -35,7.07 7.07)) COMPOUNDCURVE ( (-7.07 7.07,37.93 52.07),CIRCULARSTRING (37.93 52.07,45 55,52.07 52.07), (52.07 52.07,89.71 13.85,127.93 52.07),CIRCULARSTRING (127.93 52.07,135 55,142.07 52.07), (142.07 52.07,187.07 7.07),CIRCULARSTRING (187.07 7.07,189.16 4.01,190 0), (190 0,190 -45),CIRCULARSTRING (190 -45,187.13 -52.01,180 -55), (180 -55,45 -55),CIRCULARSTRING (45 -55,41.27 -54.28,37.93 -52.07), (37.93 -52.07,-7.07 -7.07))
MULTILINESTRING ( (0 0,1 1,1 2), (2 3,3 2,5 4)) 0.5 0 MULTILINESTRING ( (1.65 2.65,3 1.29,5.35 3.65)) MULTILINESTRING ( (2.35 3.35,3 2.71,4.65 4.35))
MULTILINESTRING ( (0 0,1 1,1 2), (2 3,3 2,5 4)) 0.5 1 MULTICURVE ( (1.65 2.65,2.65 1.65),CIRCULARSTRING (2.65 1.65,3 1.5,3.35 1.65), (3.35 1.65,5.35 3.65)) MULTILINESTRING ( (2.35 3.35,3 2.71,4.65 4.35))

This function depends on some other structures and functions which are included in the following file PostGIS_Parallel.sql

I hope this helps someone.

Creative Commons License

post this at del.icio.uspost this at Diggpost this at Technoratipost this at Redditpost this at Farkpost this at Yahoo! my webpost this at Windows Livepost this at Google Bookmarkspost this to Twitter

Comment [2]

Hi Simon,

Thanks for sharing this. I’m curious about the workings of this functions but can’t get it to work.
While adding the functions postgres tells me ‘ERROR: type “coordtype” does not exist’
Could that be a type of your own?

Cheers,
Tom

— Tom · 21 July 2011, 21:13 · #

Tom,

I think as Simon noted in the article he doesn’t have the helper functions included in this article, but will in a later article (I do recall the coord type he has mentioned in ST_DumpPoints article or something of that sort)

Simon,

Just as a note to your PostGIS readers. If you are using upcoming PostGIS 2.0 alpha, there is a new function called ST_OffsetCurve that does something similar. Check it out.

http://www.postgis.org/documentation/manual-svn/ST_OffsetCurve.html
(pictures will be up soon if they aren’t by the time you read this)
For windows users, we have experimental binaries available for PostGIS 2.0 – PostgreSQL 8.4-9.1beta3

http://www.postgis.org/download/windows/experimental.php#PostGIS_2_0_0

PostGIS 2.0 is still a work in progress so expect it to change before release (Release hopefully in like 2 months from now).

— Regina · 25 July 2011, 22:28 · #