How to create an elevation profile

Q:

How to create an elevation profile from intersecting shapefiles

A:

You want to create an elevation profile for a walking track described in one file, when the only available elevation data is contained in another file describing the terrain contours for the same region as the track. It is assumed both files are available to GRASS in shapefile format. In order to ascertain the elevation along the length of the track, you must calculate the elevation at the points where the track intersects with the contour lines, whose elevation is known. GRASS can handle data in both vector and raster form. In this instance we will start with the two shapefiles in vector form, one representing the coordinates of the track, the other the coordinates and the elevation of the contour lines. We will then convert them to raster form in order to establish the elevation at the points where the track intersects the contours, and finally produce the track elevation profile data.

1. Importing the shapefiles into GRASS

v.in.ogr -o dsn=/home/george/GRASSDATA/Contour.shp \
   output=contours_vect min_area=0.0001 type=line snap=-1

v.in.ogr -o dsn=/home/george/GRASSDATA/Track_cl.shp \
   output=tracks_vect min_area=0.0001 type=line snap=-1

At this point you can optionally visualize the map display using the GRASS Display Manager which opens up automatically when you start GRASS. You can also check that the elevation values have been correctly imported using: v.db.select using elevation as the property to list.

2. Set the computational region resolution and convert vector to raster

In order to convert the vector files to raster we will use the v.to.rast program. However before you proceed you MUST ensure that the computational region is set to the correct resolution. The v.to.rast documentation has an innocent looking warning which reads as follows:

"v.to.rast will only affect data in areas lying inside the boundaries of the current geographic region. Before running v.to.rast, the user should therefore ensure that the current geographic region is correctly set and that the region resolution is at the desired level."

Do not ignore this warning as the most likely outcome will be a probably useless raster map consisting of large, chunky squares for each pixel. To change the resolution of the computational region use the g.region command with an appropriate resolution. Some experimentation might be required. For example use:

g.region -s -p res=10

which will create a resolution where every pixel on the map represents 10 metres on the ground. Then from the Map Display window icons bar click the "zoom to..." icon (a magnifying glass over a small map) and from the drop menu select "Set computational region extents to match display". That's it. You can now safely convert your vector files to raster, as follows:

v.to.rast input=contours_vect output=contours_rast use=attr \
   type=line layer=1 column=ELEVATION value=1 rows=4096 

v.to.rast input=tracks_vect output=tracks_rast use=val \
   type=line layer=1 value=1 rows=4096

3. Creating the elevation profile

r.surf.contour input=contours_rast output=elevmap_raster

Few seconds later it should have completed its job without errors. (See the Contour lines to DEM wiki page for details.)

Then go ahead with:

r.profile -i input=elevmap_rast output=- null=*

The -i options allows interactive following of the track on the display map by clicking with the left mouse button along the section of the track you are interested in, and a two column list of values will be printed on stdout, containing the distance between points and the elevation. Use any spreadsheet application to plot the values. Alternatively, you can omit the -i option and pipe the list of coordinates to r.profile as follows:

The *profile* parameter can be set to comma separated geographic coordinates for profile line endpoints. (The interactive flag (*-i*) overrides this option). Alternatively the coordinate pairs can be piped from stdin, one comma separated pair per line.

You can get the coordinate pairs from your track vector with v.out.ascii, and use the fs=, option to make the output file "comma-separated".

4. An alternative has been suggested (not tested yet):

Use the command v.drape to get profiles on raster maps along lines from vector maps. The command to be executed is something like.

v.drape input=perfil out=perfil_n_mg type=point rast=eigen \
   method=cubic

And to convert the profiles to ascii files, like tables, you can run:

v.out.ascii perfil_n_mg > perfil_n_mg.txt

An example using this method is described here: http://casoilresource.lawr.ucdavis.edu/drupal/node/375