cumulatively calculate the distance along a line

Question

I am using Postgis to:

• calculate the cumulative distance in metres along the line
• store the cumulative distance in "m"

I have a table with points on this line and each row holds lat and lon coordinates. The "m" value is 0 for all rows.

The following code nicely gets me the distance between the first point and the next.

SELECT a.geom
FLOOR(
ST_Distance(ST_Transform(a.geom, 3857), ST_Transform(b.geom, 3857)) * cosd(42.3521))
AS dist
FROM line a
LEFT JOIN line b
ON a.id = b.id + 1

however, it does not add up to a total nor does it update the "m" value yet.

I tried a set of permutations of the following, but none of them store a cumulative sum of the previous row's value for "m" plus the calculation of the distance between this point and the previous.

WITH next AS (
    SELECT
    ST_Distance(ST_Transform(a.geom, 3857), ST_Transform(b.geom, 3857)) * cosd(42.3521)
    AS dist,
    a.id AS rowid,
    FROM line a
    JOIN line b
    ON a.id = b.id + 1
)
UPDATE line a
SET m = FLOOR(next.m + next.dist)
FROM next
WHERE a.id = next.rowid
RETURNING a.m, next.dist;

or

update line a
    set m = FLOOR(a.m + prev.dist) 
    from (select
            l.*,
            ST_Distance(ST_Transform(
                lag(geom) over (order by l.id asc)
            , 3857), ST_Transform(geom, 3857)) * cosd(42.3521) as dist
            from line l
          ) AS prev
    where prev.id = a.id - 1
    returning prev.m, prev.dist, a.m

Any ideas?

Answer 1

Fiddling around with cosd() is silly, especially when you do not know the actual latitude. ST_Distance can compute the desired value directly if you give it geography values; if you have lon/lat geometries, you can simply cast them.

To add up the lengths step by step, use a recursive CTE:

WITH RECURSIVE cum(id, geom, dist) AS (
  SELECT id,
         geom,
         0::double precision
  FROM line
  WHERE id = 1

  UNION ALL

  SELECT line.id,
         line.geom,
         cum.dist + ST_Distance(cum.geom::geography, line.geom::geography)
  FROM cum
  JOIN line ON cum.id + 1 = line.id
)
UPDATE line
SET m = dist
FROM cum
WHERE line.id = cum.id;