Creating perpendicular line transects in PostGIS

Question

I have a river network (Line with arrow as shown in image below) and I need to create perpendicular lines to the river at every 15 meters as you walk along it (see image below).

Is there a function in PostGIS capable of doing this?

Answer 1

No clear existing function to do this, but working back from the great suggestions on this post: Making perpendicular lines through each polygon edge using PostGIS?

Assuming you are good with PostGIS syntax, this might help you:

1. Split the line into 15m segments, like this: https://www.northrivergeographic.com/qgis-split-a-line-at-an-interval
2. Loop every segment in the newly segmented line and for each segment:
3. Create parallel lines with ST_OffsetCurve with a desired distance D (can be as longer than the lake/sea as shown above, we will clip them later), on both sides of the line
4. For every 15m segment, make a new line from the centroid of that segment to the centroid of the corresponding parallel lines with ST_MakeLine
5. You should have perpendicular lines, spaced 15m extending beyond the river basin/lake
6. Clip the lines with the sea region

Answer 2

Unfortunately, I can't see your figure and as a result, I suggest to create the following pgSQL geo-function for implementation, as I perceived the content of the question.

There are many simple and complex geo-tools that can solve your complex problem and this is one of them.

1. The first function can solve the problem for rectilinear objects.

Create the geo-function:

CREATE OR REPLACE FUNCTION ST_TransectsOnStraightLine(
    geom GEOMETRY,
    sl numeric)
   RETURNS GEOMETRY AS 
$BODY$
    WITH
        geodata AS (SELECT row_number() over() AS id, (ST_Dump($1)).geom geom),
        linecut AS (SELECT id, ST_LineSubstring(d.geom, substart, CASE WHEN subend > 1 THEN 1 ELSE subend END) geom FROM (SELECT id, geom, ST_Length(((geom)::geometry)) len, $2 sublen FROM geodata) AS d CROSS JOIN LATERAL (SELECT i,  (sublen * i)/len AS substart, (sublen * (i+1))/len AS subend
    FROM generate_series(0, floor(d.len/sublen)::integer) AS t(i) WHERE (sublen * i)/len <> 1.0) AS d2),
        rotate AS (SELECT id, ST_Rotate(ST_Collect(geom), -pi()/2, ST_LineInterpolatePoint(geom, 0.5)) geom FROM linecut GROUP BY id, geom)
        SELECT ST_Union(geom) geom FROM rotate;
$BODY$
LANGUAGE SQL

Run for your geodata (the function action is tested for geodata located in a flat rectangular coordinate system).

SELECT ST_TransectsOnStraightLine(geom, 15) geom FROM <name_line_table>

As a result you should get a geo image with perpendicular transects of 15 meters long, as shown in the figure below.

Picture

Note: If necessary, crop them with the ST_Buffer function with the required size.

2. The second function is capable of solving the problem with slight simplifications for complex convoluted linear objects.

CREATE OR REPLACE FUNCTION ST_TransectsOnCurveLine(
    geom GEOMETRY,
    sl numeric)
   RETURNS GEOMETRY AS 
$BODY$
WITH
    geodata AS (SELECT row_number() over() AS id, (ST_Dump($1)).geom geom),
    linecut AS (SELECT row_number() over() AS id, ST_LineSubstring(d.geom, substart, CASE WHEN subend > 1 THEN 1 ELSE subend END) geom FROM (SELECT id, geom, ST_Length(((geom)::geometry)) len, $2 sublen FROM geodata) AS d CROSS JOIN LATERAL (SELECT i,  (sublen * i)/len AS substart, (sublen * (i+1))/len AS subend
    FROM generate_series(0, floor(d.len/sublen)::integer) AS t(i) WHERE (sublen * i)/len <> 1.0) AS d2),
    rotate1 AS (SELECT id, ST_Rotate(ST_Collect(geom), -pi()/2, ST_LineInterpolatePoint(geom, 0.5)) geom FROM linecut GROUP BY id, geom),
    tr1 AS  (SELECT id,  ST_MakeLine(ST_StartPoint(geom), ST_EndPoint(geom)) geom FROM (SELECT id, (ST_Dump(geom)).geom geom FROM rotate1) foo),
    rotate2 AS (SELECT id, ST_Rotate(ST_Collect(geom), pi()/2, ST_LineInterpolatePoint(geom, 0.5)) geom FROM linecut GROUP BY id, geom),
    tr2 AS (SELECT id, ST_MakeLine(ST_StartPoint(geom), ST_EndPoint(geom)) geom FROM (SELECT id, (ST_Dump(geom)).geom geom FROM rotate2) foo),
    line_incoherence AS (SELECT a.id, ST_MakeLine(ST_StartPoint(a.geom), ST_EndPoint(b.geom)) geom FROM tr1 a JOIN tr2 b ON true AND a.id=b.id), 
    int_pnt AS (SELECT a.id, ST_ShortestLine(ST_LineInterpolatePoint(a.geom, 0.5), b.geom) geom FROM line_incoherence a JOIN linecut b ON a.id=b.id GROUP BY a.id, a.geom, b.geom ORDER BY id), 
    transectline AS (SELECT id, ST_Rotate(ST_Collect(geom), -pi(), ST_EndPoint(geom)) geom FROM int_pnt GROUP BY id, geom)
    SELECT ST_Union(geom) geom FROM (SELECT geom FROM int_pnt UNION SELECT geom FROM transectline) foo;
$BODY$
LANGUAGE SQL

Run

SELECT ST_TransectsOnCurveLine(geom, 15) geom FROM <name_line_table>

As a result you should get a geo image with perpendicular transects of 15 meters long, as shown in the figure below.

Picture

(-: CM-FOGS :-)

Answer 3

Here's another function to create transects, providing control over the transect length as well as the transect separation length

It operates as follows:

• split the input line into sections of length secLen
• for each section compute the transect angle and centre point of the section (the angle is Pi - azimuth(startPt, endPt))
• compute the transect lines from two points offset from the centre point

CREATE OR REPLACE FUNCTION ST_Transects(
    lineGeom geometry,
    secLen float,
    transectLen float)
    RETURNS geometry 
    LANGUAGE sql AS 
$BODY$
WITH
geodata AS (SELECT lineGeom, ST_Length(lineGeom) AS lineLen),
sections AS (SELECT ST_LineSubstring(lineGeom, secStart, CASE WHEN secEnd > 1 THEN 1 ELSE secEnd END) geom
    FROM geodata AS t
    CROSS JOIN LATERAL 
        (SELECT i * secLen/lineLen AS secStart, (i+1) * secLen/lineLen AS secEnd
            FROM generate_series(0, floor(lineLen / secLen)::integer) AS t(i) 
            WHERE (secLen * i)/lineLen <> 1.0) AS t2 ),
sectAnglePt AS (SELECT pi() - ST_Azimuth(ST_StartPoint(geom), ST_EndPoint(geom)) AS ang,
                  ST_LineInterpolatePoint(geom, 0.5) AS centre
                  FROM sections)
SELECT ST_Collect(ST_MakeLine(
                    ST_Point( ST_X(centre) - transectLen * cos(ang), 
                              ST_Y(centre) - transectLen * sin(ang)),
                    ST_Point( ST_X(centre) + transectLen * cos(ang), 
                              ST_Y(centre) + transectLen * sin(ang)) )) AS geom 
  FROM sectAnglePt;
$BODY$;

Call it like:

SELECT ST_Transects( line, 0.01, 0.005);