I'm looking to generate elevation profiles from a number of GPX traces (which I also have GeoJSON conversions of). I have 20m resolution DEMs of the area, in GeoTIFF format. Ultimately I want to generate a static image that combines a number of these elevation profiles, but I'll be happy if I can get something like a CSV of elevation values sampled at regular intervals.

Summarised, how can I turn:

  • 1+ vector paths (.gpx or .geojson)
  • 1 DEM (.tif)

into:

  • 1+ .svg (raw plot, no axes or labels); or
  • 1+ .csv

using only command line tools. I have managed to achieve this manually using QGIS (and the Profile Tool plugin), but need to automate it. I'd prefer to avoid writing code - but Python is much preferable to R. A workflow that added elevations to the GeoJSON file would be better than nothing. I'm really not looking to implement a sampling algorithm from scratch.

OS X or Ubuntu are ok.

(Note to moderators: there seem to be lots of similar questions, but many of them are just "how do I get an elevation profile from this GPX file". The key differences are that I have the suitable terrain file already, and I need to automate this for dozens of files.)

  • I've since discovered the Mapbox Surface API. I believe the source is open but I don't know how to deploy it locally. – Steve Bennett Jul 12 '15 at 13:17
up vote 2 down vote accepted

If you load your LineString gpx/geojson

ogr2ogr -f "PostgreSQL" PG:"dbname=" linestring.geojson

and DEM tiff into PostGIS (replace XXXX with EPSG code, eg for WGS84 use 4326)

raster2pgsql -d -C -I -M -s XXXX -t auto dem.tiff dem

you can then use ST_Segmentize to add extra vertices at the resolution of your DEM (gdalinfo's Pixel Size), then ST_DumpPoints to convert the LineString to a bunch of points,

CREATE TABLE line_points AS select (ST_DumpPoints(ST_Segmentize(wkb_geometry, 5))).*, 0.0 AS ele from line;

then ST_Value to get the elevation at each point,

CREATE INDEX ON line_points USING gist (geom);
UPDATE line_points SET ele = (SELECT ST_Value(dem.rast, geom) AS ele FROM dem WHERE ST_Intersects(line_points.geom, dem.rast) ORDER BY ele LIMIT 1);

then update the geometry to include the elevation

UPDATE line_points SET geom = ST_MakePoint(ST_X(geom), ST_Y(geom), ele);

then ST_MakeLine to reconstruct the points with elevation back into a LineString.

CREATE AGGREGATE array_accum (anyelement)
(
    sfunc = array_append,
    stype = anyarray,
    initcond = '{}'
);
CREATE TABLE linez AS SELECT ST_MakeLine(array_accum(geom)) AS geom FROM (SELECT * FROM line_points ORDER BY path) AS line;

Export back to GeoJSON/CSV/etc:

ogr2ogr -f GeoJSON linez.geojson PG: linez

This could be optimised to generate multiple profiles in bulk, so you're only calling each SQL statement once.

  • Wow, nice answer. Thanks Andrew! – Steve Bennett Jul 26 '17 at 1:20

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.