[postgis-users] splitting 4D (XYZM) linestrings and interpolating/retaining M

[Phil Hurvitz]

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
phurvitz at u.washington.edu | 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.
>