Adding elevation to a line from a DEM in PostGIS and maintaining accurate measures

Posted by: Gavin Fleming | in PostGIS, QGIS | 2 months, 4 weeks ago | Comments

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 all about meeting two requirements that the team planning the marathon routes had:

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 identify tool function that shows the X, Y, Z and M coordinates of the nearest vertex in the 'derived' section, like this:

That point is 2.27km along the 5km fun run route and 13m above sea level. Note: In QGIS 3.4 this will also show the measure (distance) interpolated to exactly where you clicked.

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

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:

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.

Currently unrated

Comments

Template by Blacktie Mezzanine theme by CodingHouse