Hi all, I figured this out, following this StackExchange post:
https://gis.stackexchange.com/questions/40622/how-to-add-vertices-to-existing-linestrings
The method is to find the location of the intersection using
ST_LineLocatePoint, which gives a value between 0 and 1, indicating the
relative position along the line where the split should occur, and then
to use ST_LineSubstring twice, to get the portions of the linestring
before and after the split.
with
--make a simple xyzm line
l as (select st_makeline(st_makepoint(0,0,0,0), st_makepoint(0,1,10,20))
as geom_l)
--an xy line for splitting
, c as (select st_makeline(st_makepoint(-1, 0.25), st_makepoint(1,
0.25)) as geom_c)
--get the point of intersection
, i as (select st_intersection(l.geom_l, c.geom_c) as geom_i from l, c
where st_intersects(l.geom_l, c.geom_c))
--locate the point of intersection along the line
, lp as (select st_linelocatepoint(l.geom_l, i.geom_i) as p from l, i)
--split the line in two based on the location
, s as (
select st_linesubstring(l.geom_l, 0, p) as geom from l, lp
union all
select st_linesubstring(l.geom_l, p, 1) as geom from l, lp
)
--text now shows interpolation in Z and M
select st_astext(geom) from s;
st_astext
----------------------------------------
LINESTRING ZM (0 0 0 0,0 0.25 2.5 5)
LINESTRING ZM (0 0.25 2.5 5,0 1 10 20)
-P.
**************************************************************
Philip M. Hurvitz, PhD | Research Assistant Professor | UW-CBE
Urban Form Lab | 1107 NE 45th Street, Suite 535 | Box 354802
University of Washington, Seattle, Washington 98195-4802, USA
[email protected] | http://gis.washington.edu/phurvitz
"What is essential is invisible to the eye." -de Saint-Exupéry
**************************************************************
On 2018-08-14 11:41, Phil Hurvitz wrote:
Hi all, I'm trying to split 4D linestrings and interpolate the M value.
Here is a sample data set:
with --make a simple xyzm line l as (select
st_makeline(st_makepoint(0,0,0,0), st_makepoint(0,1,0,20)) as g0)
--show as text, clearly shows M value select st_asewkt(g0) from l;
st_asewkt
------------------------------
LINESTRING(0 0 0 0,0 1 0 20)
Splitting with a point as the blade interpolates the M value. If I
split the line at 25% along its length, the M value is interpolated
appropriately.
--split by point
with
--make a simple xyzm line
l as (select st_makeline(st_makepoint(0,0,0,0),
st_makepoint(0,1,0,20)) as g0)
--an xy point
, p as (select st_makepoint(0, 0.25) as g1)
--split
, s as (select st_split(l.g0, p.g1) as g2 from l, p)
--show as text, clearly shows first element ends and second element
starts at M = 5, which is 25% along
select st_asewkt((st_dump(g2)).geom) from s;
st_asewkt
---------------------------------
LINESTRING(0 0 0 0,0 0.25 0 5)
LINESTRING(0 0.25 0 5,0 1 0 20)
However, splitting with a line seems to cause the M value to disappear
(set to zero):
--split by line
with
--make a simple xyzm line
l as (select st_makeline(st_makepoint(0,0,0,0),
st_makepoint(0,1,0,20)) as g0)
--an xy line for splitting
, c as (select st_makeline(st_makepoint(-1, 0.25), st_makepoint(1,
0.25)) as g1)
--split
, s as (select st_split(l.g0, c.g1) as g2 from l, c)
--show as text, M is now zero
select st_asewkt((st_dump(g2)).geom) from s;
st_asewkt
----------------------------
LINESTRING(0 0 0,0 0.25 0)
LINESTRING(0 0.25 0,0 1 0)
Also the same happens splitting with a polygon, in this case, a buffer:
--split with buffer
with
--make a simple xyzm line
l as (select st_makeline(st_makepoint(0,0,0,0),
st_makepoint(0,2,0,20)) as g0)
--a buffer for splitting
, b as (select st_force4d(st_buffer(st_makepoint(0,1,0,5), 0.5, 10))
as g1)
--split with buffer
, x as (select st_split(l.g0, b.g1) as g2 from l, b)
--M goes to zero
select st_asewkt((st_dump(g2)).geom) from x;
st_asewkt
-----------------------------
LINESTRING(0 0 0,0 0.5 0)
LINESTRING(0 0.5 0,0 1.5 0)
LINESTRING(0 1.5 0,0 2 0)
Can anyone provide advice on how to split the 4D line with a non-point
geometry and calculate interpolated M values?
Thanks for any guidance.
_______________________________________________
postgis-users mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/postgis-users
