Top 5 Recent Articles

ARTICLES CATEGORIES

ST_Parallel for PostGIS

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.

 -- Drop function
 DROP FUNCTION IF EXISTS st_parallel(geometry, NUMERIC, INTEGER, INTEGER);
 -- Create Function
 CREATE OR REPLACE FUNCTION st_parallel(p_geometry     geometry,
                                          p_distance     NUMERIC,
                                          p_round_factor INTEGER,
                                          p_curved       INTEGER DEFAULT 0)
   RETURNS geometry AS
 $BODY$
 DECLARE
     v_rowcount        INTEGER;
     v_part_no         INTEGER;
     v_result          geometry;
     bAcute            BOOLEAN := FALSE;
     bCurved           BOOLEAN := ( 1 = p_curved );
     v_dims            INTEGER;
     v_delta           CoordType;
     v_prev_delta      CoordType;
     v_adj_Coord       CoordType;
     v_int_1           CoordType;
     v_int_coord       CoordType;
     v_int_2           CoordType;
     v_startCoord      CoordType;
     v_endCoord        CoordType;
     v_distance        NUMERIC;
     v_ratio           NUMERIC;
     v_az              NUMERIC;
     v_dir             INTEGER;
     v_return_geom     geometry;
     v_geom_string     text;
   c_parts CURSOR ( p_geom         geometry,
                    p_Geometrytype text )
   IS
     SELECT p_geom AS geom
      WHERE p_geometrytype IN ('ST_LineString' )
     UNION ALL
     SELECT ST_GeometryN(p_geom,generate_series(1,ST_NumGeometries(p_geom))) AS geo m
      WHERE p_GeometryType = 'ST_MultiLineString';
   c_vector CURSOR ( p_linestring geometry )
   IS
     SELECT startCoord, endCoord
       FROM (SELECT (ST_GetVector(p_linestring)).* ) o;
 BEGIN
   /* Problem with this algorithm is that any existing circular curved elements will not be honoured unless bCurved is 0 */
   IF ( p_round_factor IS NULL ) THEN
      raise exception 'p_round_factor may not be null';
   END IF;
   IF ( p_geometry IS NULL
        OR
        ST_GeometryType(p_geometry) NOT IN ('ST_MultiLineString','ST_LineString') ) THEN
      Raise Exception 'p_linestring is null or is not a linestring or multilinestring';
   END IF;
   IF ST_HasArc(p_geometry) THEN /* Should never happen as GeometryType captures anything with a curve */
      Raise Exception 'Compound linestrings are not supported.';
   END IF;
   /* Process all parts of a, potentially, multi-part linestring geometry */
    v_part_no := 0;
    FOR part IN c_parts(p_geometry,ST_GeometryType(p_geometry)) loop
     BEGIN
       v_part_no := v_part_no + 1;
       v_dims := ST_NDims(part.geom);
       IF ( v_part_no > 1 ) THEN
          v_geom_string := v_geom_string || CASE WHEN bCurved THEN ',' ELSE ',(' END;
       END IF;
       /* Process geometry one vector at a time */
       v_rowcount := 0;
       FOR rec IN c_vector(part.geom) LOOP
         v_rowcount := v_rowcount + 1;
         /* Compute base offset */
         v_startCoord := rec.startCoord;
         v_endCoord   := rec.endCoord;
         v_az := ST_Azimuth(ST_MakePoint(v_startCoord.X,v_startCoord.Y),
                            ST_MakePoint(v_endCoord.X,  v_endCoord.Y  ) );
         v_dir     := CASE WHEN v_az < pi() THEN -1 ELSE 1 END;
         v_delta.x := ABS(COS(v_az)) * p_distance * v_dir;
         v_delta.y := ABS(SIN(v_az)) * p_distance * v_dir;
         IF  NOT ( v_az > pi()/2 AND
                   v_az < pi() OR
                   v_az > 3 * pi()/2 ) THEN
           v_delta.x := -1 * v_delta.x;
         END IF;
         /* merge vectors at this point? */
         IF (v_ROWCOUNT > 1) THEN
            /* Get intersection of two lines parallel at distance p_distance from current ones */
            v_int_coord := rec.startCoord;
            SELECT *
              FROM findlineintersection(v_adj_coord.x,                   v_adj_coord.y,
                                        v_startCoord.x + v_prev_delta.x, v_startCoord.y + v_prev_delta.y,
                                        v_endCoord.x   + v_delta.x,      v_endCoord.y   + v_delta.y,
                                        v_startCoord.x + v_delta.x,      v_startCoord.y + v_delta.y
              INTO v_int_coord.x, v_int_coord.y,
                   v_int_1.x,     v_int_1.y,
                   v_int_2.x,     v_int_2.y);
            IF ( v_int_coord.x = 1E+308 ) THEN
              /* No intersection point as lines are parallel */
              bAcute := TRUE;
              /* int coord could be computed from start or end, doesn't matter. */
              v_int_coord.x := Round((v_startCoord.x + v_prev_delta.x)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
              v_int_coord.y := Round((v_startCoord.y + v_prev_delta.y)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
              v_int_coord.z := v_startCoord.z::DOUBLE PRECISION;
              v_int_coord.m := v_startCoord.m::DOUBLE PRECISION;
            ELSE
              bAcute := ( Round(v_int_coord.x::NUMERIC,p_round_factor::INTEGER) = Round(v_int_1.x::NUMERIC,p_round_factor::INTEGER) )
                    AND ( Round(v_int_coord.x::NUMERIC,p_round_factor::INTEGER) = Round(v_int_2.x::NUMERIC,p_round_factor::INTEGER) )
                    AND ( Round(v_int_coord.y::NUMERIC,p_round_factor::INTEGER) = Round(v_int_1.y::NUMERIC,p_round_factor::INTEGER) )
                    AND ( Round(v_int_coord.y::NUMERIC,p_round_factor::INTEGER) = Round(v_int_2.y::NUMERIC,p_round_factor::INTEGER) );
            END IF;
            IF ( bCurved AND NOT bAcute) THEN
              /* First point in "intersection" curve */
              v_int_1.x := Round((v_startCoord.x + v_prev_delta.x)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
              v_int_1.y := Round((v_startCoord.y + v_prev_delta.y)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
              v_int_1.z := v_startCoord.z;
              v_int_1.m := v_startCoord.m;
              /* Intersection point (top of circular arc) */
              /* Need to compute coordinates at mid-angle of the circular arc formed by last and current vector */
              v_distance := ST_Distance(ST_SetSRID(ST_Point( v_int_coord.x, v_int_coord.y),ST_SRID(p_geometry)),
                                        ST_SetSRID(ST_Point(v_startCoord.x,v_startCoord.y),ST_SRID(p_geometry)));
              v_ratio := ( p_distance / v_distance ) * SIGN(p_distance);
              v_adj_coord   := rec.startCoord;
              v_adj_coord.x := Round((v_startCoord.x + (( v_int_coord.x - v_startCoord.x ) * v_ratio ))::NUMERIC,p_round_factor::INTEGER);
              v_adj_coord.y := Round((v_startCoord.y + (( v_int_coord.y - v_startCoord.y ) * v_ratio ))::NUMERIC,p_round_factor::INTEGER);
              /* Last point in "intersection" curve */
              v_int_2.x := Round((v_startCoord.x + v_delta.x)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
              v_int_2.y := Round((v_startCoord.y + v_delta.y)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION;
              v_int_2.z := v_startCoord.z::DOUBLE PRECISION;
              v_int_2.m := v_startCoord.m::DOUBLE PRECISION;
            ELSE
              /* Intersection point */
              v_adj_coord   := v_int_coord;
              v_adj_coord.x := Round(v_int_coord.x::NUMERIC,p_round_factor::INTEGER);
              v_adj_coord.y := Round(v_int_coord.y::NUMERIC,p_round_factor::INTEGER);
            END IF;
         ELSE  /* c_vector%ROWCOUNT = 1 */
            v_geom_string := /*GeometryType(p_geometry) || */ '(';
           /* Translate start point with current vector */
           v_adj_coord.x := Round((v_startCoord.x + v_delta.x)::NUMERIC, p_round_factor::INTEGER)::DOUBLE PRECISION;
           v_adj_coord.y := Round((v_startCoord.y + v_delta.y)::NUMERIC, p_round_factor::INTEGER)::DOUBLE PRECISION;
           v_adj_coord.z := v_startCoord.z::DOUBLE PRECISION;
           v_adj_coord.m := v_startCoord.m::DOUBLE PRECISION;
         END IF;
         /* Now add computed coordinates to output geometry */
         IF (NOT bCurved) OR bAcute OR (v_ROWCOUNT = 1 ) THEN
       v_geom_string := v_geom_string ||
                        CoordAsText(v_adj_coord,v_dims)::text ||
                        ','::text;
         ElsIf (bCurved) THEN
           /* With any generated circular curves we need to add the appropriate type */
       v_geom_string := v_geom_string ||
                        CoordAsText(v_int_1,v_dims)::text ||
                        '),CIRCULARSTRING('::text || CoordAsText(v_int_1,    v_dims)::text || ','::text ||
                                                     CoordAsText(v_adj_coord,v_dims)::text || ','::text ||
                                                     CoordAsText(v_int_2,    v_dims)::text || ')'::text ||
                        ',(' || CoordAsText(v_int_2,v_dims)::text || ','::text;
         END IF;
         v_prev_delta := v_delta;
       END LOOP;
       /* Append Point at End of Linestring */
       v_geom_string := v_geom_string ||
              CoordAsText(create_coordtype(Round((v_endcoord.x + v_delta.x)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION,
                                           Round((v_endcoord.y + v_delta.y)::NUMERIC,p_round_factor::INTEGER)::DOUBLE PRECISION,
                                           v_endcoord.z::DOUBLE PRECISION,
                                           v_endcoord.m::DOUBLE PRECISION),v_dims)::text || ')';
     END;
   END Loop;
   IF (bCurved AND v_geom_string LIKE '%CIRCULARSTRING%' ) THEN
       v_geom_string := CASE WHEN ST_GeometryType(p_geometry) = 'ST_MultiLineString'
                             THEN 'MULTICURVE('
                             ELSE 'COMPOUNDCURVE('
                         END || v_geom_String || ')';
    ELSE
       v_geom_string := CASE WHEN ST_GeometryType(p_geometry) = 'ST_MultiLineString'
                             THEN GeometryType(p_geometry) || '(' || v_geom_String || ')'
                             ELSE GeometryType(p_geometry) || v_geom_String
                         END ;
    END IF;
   v_result := ST_GeometryFromText(v_geom_string,ST_SRID(p_geometry));
   RETURN v_Result;
 END;
 $BODY$
   LANGUAGE plpgsql IMMUTABLE STRICT
   COST 100;

Here is a test:

 WITH geometries AS (
           SELECT ST_GeomFromText('LINESTRING(1 1,1 10)') AS geom,
                   10.0 AS offset,2 AS roundFactor,0 AS curved
 UNION ALL SELECT ST_GeomFromText('LINESTRING(0 0,1 1,1 2)') AS geom,
                    0.5 AS offset,2 AS roundFactor,
                   generate_series(0,1,1) AS curved
 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,
                   10.0 AS offset,2 AS roundFactor,
                   generate_series(0,1,1) AS curved
 UNION ALL SELECT ST_GeomFromText('MULTILINESTRING((0 0,1 1,1 2),(2 3,3 2,5 4))') AS geom,
                    0.5 AS offsetRight,2 AS roundFactor,
                   generate_series(0,1,1) AS curved
 )
 SELECT ST_AsText(g.geom) AS origGeom, g.offset,g.curved,
         ST_AsText(ST_Parallel(g.geom,g.offset,g.roundFactor,g.curved)) AS geomWKTLeft,
         ST_AsText(ST_Parallel(g.geom,0.0-g.offset,g.roundFactor,g.curved)) AS geomWKTRight
   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.