Kartoza - Adding Elevation to a Line from a DEM in PostGIS and Maintaining Accurate Measures

This is the second in a three part series on the behind-the-scenes GIS work that can go into planning a complex event, in this case the Cape Town Marathon.

Following on from creating point markers along a line in PostGIS this post delves into 4D coordinates, PostGIS functions and rasters in PostGIS.

What I show here is about meeting two requirements of the marathon route planning team:

• Automatically update the length of a route after each edit with the most accurate calculation possible
• Enable a QGIS user to query the distance along a line at any point by clicking on it

For the first one, we want to calculate the length along the spheroid, which is easy enough, but we can make it even more accurate by taking into account elevation as the route goes up and down hills.

For the second one, we want to take advantage of the QGIS 3.4 identify tool function that shows the X, Y, Z and M coordinates of the nearest vertex as well as the interpolated M and Z values in the 'derived' section, like this:

That point is 15.74km along the marathon route and 7.5m above sea level.

We are going to update the Z (elevation) and M (measure) values of every vertex every time the geometry is edited and saved.

For elevation we're going to use a raster DEM (digital elevation model) within the database. The City of Cape Town has an accessible, open data portal and I found a 10m DEM there. So I downloaded it and then loaded it into my PostGIS database:

raster2pgsql -C -I -M -F -Y -s 40019 -t 100x100 10m_BA.tif dem | psql -d mydb

The marathon and long trail routes in QGIS 3D

Then I had to make sure that each route's geometry field catered for 4 dimensions:

ALTER TABLE marathon ALTER COLUMN geom TYPE geometry(LinestringZM,4326) USING ST_Force4D(geom);

Next I created this function:

--you need a DEM called 'dem' in SRID WGS19 loaded in the database for this to work
--adapted from http://blog.mathieu-leplatre.info/drape-lines-on-a-dem-with-postgis.html
CREATE OR REPLACE FUNCTION update_ZM()
  RETURNS trigger AS
$BODY$
BEGIN
 EXECUTE format('WITH line AS
    (SELECT ST_Transform($1.geom,40019) as geom),
  points2d AS
    (SELECT (ST_DumpPoints(geom)).geom AS geom FROM line),
  cells AS
    (SELECT p.geom AS geom, ST_Value(d.rast, 1, p.geom) AS val
     FROM %I.dem d, points2d p
     WHERE ST_Intersects(d.rast, p.geom)),
  points3d AS
    (SELECT ST_SetSRID(ST_MakePoint(ST_X(geom), ST_Y(geom), val), 40019) AS geom FROM cells),
  line3d AS
    (SELECT ST_MakeLine(ST_Transform(geom,4326)) AS geom FROM points3d),
  line4d AS
    (SELECT ST_AddMeasure(geom,0,st_length(geom::geography)/1000) AS geom FROM line3d)
  UPDATE %I.%I a
     SET geom = line4d.geom
     FROM line4d
     WHERE a.id = $1.id', TG_TABLE_SCHEMA, TG_TABLE_SCHEMA, TG_TABLE_NAME) USING NEW;                    
 RETURN NEW;
END;
$BODY$
LANGUAGE plpgsql;

The update_ZM()  function does the following:

• projects the line to the local CRS
• extracts all the vertices as points
• intersects the points with the DEM
• adds the elevation from the DEM to the Z value of the point coordinates
• builds the points back into a line
• adds accurate 3D distances to the M value of the coordinates as described in more detail in the previous article
• replaces the previous version of the line with the one generated by the above process

The function is a dynamic trigger function, which means it will run when a trigger is fired and work on any table that has a line geometry.

The last step is to add the trigger to any line table you want to run this function on:

DROP TRIGGER update_zm ON marathon;
CREATE TRIGGER update_zm
  AFTER INSERT OR UPDATE OF geom
  ON marathon
  FOR EACH ROW
  WHEN (pg_trigger_depth() = 0)
  EXECUTE PROCEDURE update_ZM();

This trigger runs the update_ZM() function on each record that you add or edit. If you are editing in QGIS, this means as soon as you hit Save.