Tuesday, 29 October 2019

Varying size buffer along line with PostGIS


I want to create a buffer along a line where the buffer width corresponds to a fraction of the length of the line, something like that:


tg


To do that my idea is to compute for each node of the line (except for the first and last node) the bisector angle between the previous and the next segment, example for the node, bisector angle in red #3:


example1


At the moment I'm just trying to create one side (a line) of the buffer:


example2


Here's my attempt:


WITH POINT_LINE AS
(

WITH SUB AS
(
-- Extract all the segments and calculate for each segment the cumulative length
WITH segments AS (
SELECT gid, (pt).path[1]-1 as path, ST_MakeLine(lag((pt).geom, 1, NULL) OVER (PARTITION BY gid ORDER BY gid, (pt).path), (pt).geom) AS geom
FROM (SELECT gid, ST_DumpPoints(geom) AS pt FROM MYLINE) as dumps
)
SELECT
gid,
sum(st_length(geom)) over (order by path asc rows between unbounded preceding and current row)-st_length(geom) as l_start,

sum(st_length(geom)) over (order by path asc rows between unbounded preceding and current row) as l_end,
ST_StartPoint(geom) as g_start,
ST_EndPoint(geom) as g_end,
path,
geom
FROM segments WHERE geom IS NOT NULL
)
SELECT
a.gid,
a.path,

a.l_end,
b.g_start,
-- HERE I try to calculate the bisector angle:
atan((ST_Y(a.g_start)-ST_Y(a.g_end))/(ST_X(a.g_start)-ST_X(a.g_end))) + asin((ST_X(a.g_end)-ST_X(a.g_start))/SQRT((ST_X(a.g_end)-ST_X(a.g_start))^2+(ST_Y(a.g_end)-ST_Y(a.g_start))^2)) as phi,
a.geom
FROM SUB a, SUB b
WHERE a.path = b.path-1
AND a.gid = b.gid
)
-- Intermediate node

SELECT gid, path+1 as path, ST_SetSRID(ST_MakePoint(ST_X(g_start)+cos(phi)*l_end,ST_Y(g_start)+sin(phi)*l_end),2056) as geom
FROM POINT_LINE
UNION
-- First node of the line
SELECT gid, path as path, ST_StartPoint(geom) as geom
FROM POINT_LINE
WHERE path = 1

Unfortunately, the node are not always in the right position, so I guess that my angle is sometimes wrong.


Do you have an idea on how to solve this problem?




Answer



I've improved the answer given by @Cyril, it is way faster, here is the solution:


WITH 
step1 AS
(SELECT gid,
ST_DumpPoints(geom) AS dump,
ST_Length(geom) AS len,
geom
FROM mylines),
step2 AS

(SELECT gid,
(dump).path[1],
ST_Buffer((dump).geom, ST_LineLocatePoint(geom, (dump).geom)*len/10 + 0.01) AS geom
FROM step1),
step3 AS
(SELECT gid,
ST_ConvexHull(ST_Union(geom, LEAD(geom) OVER(PARTITION BY gid ORDER BY gid, path))) AS geom
FROM step2)

SELECT gid, ST_Union(geom) AS geom FROM step3 GROUP BY gid


What this request does ?


Step 1:


SELECT gid, ST_DumpPoints(geom) AS dump, ST_Length(geom) AS len, geom FROM mylines

We extract every points of each line. I keep the GID so I will be able to perform the same operation on several line at the same time.


Step 2:


SELECT gid, (dump).path[1], st_buffer((dump).geom, ST_LineLocatePoint(geom, (dump).geom)*len/10 + 0.01) AS geom FROM step1

We apply a buffer to each point, the buffer size correspond to the corresponding length of the line at this point divided by ten (ten is arbitrary). I need to add a small constant to the buffer size (here 0.01) so the first buffer size is not 0. We keep the path (order) of each point.



Step 3:


SELECT gid, ST_ConvexHull(ST_Union(geom, LEAD(geom) OVER(PARTITION BY gid ORDER BY gid, path))) AS geom FROM step2

We use ST_ConvexHull for each consecutive pair of circular buffer.


SQL FUNCTION:


CREATE OR REPLACE FUNCTION ST_VariableBufferFromLine(
geom GEOMETRY,
Length_BufferSize_Ratio NUMERIC
)


RETURNS GEOMETRY AS

$BODY$

WITH
step1 AS
(SELECT ST_DumpPoints(geom) AS dump,
ST_Length(geom) AS len,
geom),
step2 AS

(SELECT (dump).path[1],
ST_Buffer((dump).geom, GREATEST(ST_LineLocatePoint(geom, (dump).geom)*len/Length_BufferSize_Ratio,0.001)) AS geom
FROM step1),
step3 AS
(SELECT
ST_ConvexHull(ST_Union(geom, LEAD(geom) OVER(ORDER BY path))) AS geom
FROM step2)
SELECT ST_Union(geom) AS geom FROM step3

$BODY$


LANGUAGE SQL;

So now we can simply use:


SELECT gid, ST_VariableBufferFromLine(geom,10.0) AS geom FROM mylines

If you need a fixed buffer size (for example the end of the line should have of buffer of 100m) we can simply replace this part in the function:


GREATEST(ST_LineLocatePoint(geom, (dump).geom)*len/Length_BufferSize_Ratio,0.001))

by this one:



GREATEST(ST_LineLocatePoint(geom, (dump).geom)*end_width,0.001))

With end_width = the buffer size at the end of the line. Don't forget to adapt the variable name.


Example result with 2 lines:


enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...