Opened 4 months ago

Last modified 6 days ago

#5776 new enhancement

ST_LineSubstring 3D support

Reported by: vla123 Owned by: pramsey
Priority: medium Milestone: PostGIS 3.6.0
Component: postgis Version: 3.4.x
Keywords: Cc:

Description

It seems that ST_LineSubstring function return wrong results for 3D linestrings.

select 
	st_3dlength(geometry) as full_line ,
	st_3dlength(st_linesubstring(geometry, 0.0, 0.50)) as half_line
from (
		select st_makeline(ARRAY[
		st_setsrid(st_pointz(0,0,0), 3857), 
		st_setsrid(st_pointz(0,2,5), 3857), 
		st_setsrid(st_pointz(0,10,10), 3857)
		]) as geometry
	) line

Output :
  
"full_line"	"half_line"
14.819145939191106	8.92290773165573

This should be explicitly stated on documentation.

Change History (4)

comment:1 by pramsey, 4 months ago

Milestone: PostGIS 3.4.3PostGIS 3.4.4

comment:2 by vla123, 4 months ago

For those who want to split 3D LINESTRING, here is a plpgsql function :

CREATE or REPLACE FUNCTION ST_3DLineSubstring(
line_geom geometry, start_fraction float, end_fraction float
)
RETURNS setof geometry AS $$
DECLARE 
srid INTEGER = st_srid(line_geom);
n_pts INTEGER;
vertex geometry;
next_vertex geometry;
vertices geometry [];
i INTEGER;
distance float = 0.0;
d float;
fraction float;
start_length float = st_3dlength(line_geom) * start_fraction;
end_length float = st_3dlength(line_geom) * end_fraction;
start_vertex geometry;
end_vertex geometry;
new_x float;
new_y float;
new_z float;
start_added bool = false;
BEGIN

if line_geom is null then return next
NULL;
else
n_pts := ST_NPoints(line_geom) -1;
for i in 1..n_pts LOOP
	vertex := ST_PointN(line_geom, i);
	next_vertex := ST_PointN(line_geom, i + 1);
	d := abs(ST_3DDistance(vertex, next_vertex));

	-- vertex between start and end
	if distance < end_length and distance > start_length then
	vertices := array_append(vertices, vertex);
	end if;
	
	-- first vertex
	if distance + d > start_length and start_added is false then
	fraction := (start_length - distance)/d;
	new_x := st_x(vertex) + (st_x(next_vertex) - st_x(vertex)) * fraction;
	new_y := st_y(vertex) + (st_y(next_vertex) - st_y(vertex)) * fraction;
	new_z := st_z(vertex) + (st_z(next_vertex) - st_z(vertex)) * fraction;
	start_vertex := st_setsrid(st_makepoint(new_x, new_y, new_z), srid);
	start_added := true;
	end if;

	-- last vertex
	if distance + d >= end_length then
	fraction := (end_length - distance)/d;
	new_x := st_x(vertex) + (st_x(next_vertex) - st_x(vertex)) * fraction;
	new_y := st_y(vertex) + (st_y(next_vertex) - st_y(vertex)) * fraction;
	new_z := st_z(vertex) + (st_z(next_vertex) - st_z(vertex)) * fraction;
	end_vertex := st_setsrid(st_makepoint(new_x, new_y, new_z), srid);
	EXIT;
	end if;

	distance := distance + d;
	
END LOOP;

vertices := array_append(vertices, end_vertex);
vertices := array_prepend(start_vertex, vertices);

return QUERY
select(st_makeline(vertices));
end if;
END;

$$ LANGUAGE plpgsql;


comment:3 by pramsey, 4 months ago

Hm. Implication is the need for an ST_3DLineSubstring because silently processing the 3D coordinates in ST_LineSubstring will have the same problem in reverse.

comment:4 by robe, 6 days ago

Milestone: PostGIS 3.4.4PostGIS 3.6.0
Type: defectenhancement
Note: See TracTickets for help on using tickets.