Another stab at creating a Tanaka-style contour map

In a recent post Anita Graser (aka underdark) showed how to create illuminated or Tanaka contours in QGIS using various functions available in the toolbox and some custom functions.

Here I want to explore a slightly different way to achieve the same, using GRASS GIS to compute the azimuth, brightness and line width. I’ll use the command line, but you can do the same using the menu in GRASS, or the corresponding GRASS functions in the QGIS processing toolbox.

The map itself is still created in QGIS, utilizing the powerful expression based symbology and layer blending mode. For this example I used the same DEM as used in the blog post Fancy shaded and scaled contour lines in QGIS.

Create contour lines and break them in segments

First define the angle of the sun. In the example below I use a variable (SA), so I can later easily change this values without having to go through the script. Next create the contour lines with r.contour and split them with v.split.

SA=320
r.contour -t input=dem output=tmp1 step=100 cut=2 --overwrite
v.split input=tmp1 output=tmp2 vertices=2 --overwrite

Create attribute table

Next we need to assign an unique category value to each line segment. To do this, use v.category to drop the old categories and then to create new ones. Next, create an attribute table, with columns for the azimuth, colour and line width (we will populate these with values in the next step).

v.category input=tmp2 output=tmp3 option=del cat=-1 --overwrite
v.category input=tmp3 output=contours option=add
v.db.addtable map=contours columns="azimuth double precision,color integer,width double precision

Compute the azimuth of the line segments

Compute the azimuth of each line segment, using the v.to.db function. Shift the azimuth values to set the angle of the sun at zero (first two db.execute statements) and rescale the range of values 180-365 to 180-0 (third db.execute statement).

v.to.db map=contours option=azimuth columns=azimuth units=degrees
db.execute sql="UPDATE contours SET azimuth = azimuth + $SA WHERE azimuth < (365-$SA)"
db.execute sql="UPDATE contours SET azimuth = azimuth - 365 + $SA WHERE azimuth >= (365-$SA)"
db.execute sql="UPDATE contours SET azimuth = 180 - (azimuth - 180) WHERE azimuth > 180"

Compute the brightness and line width values

Next compute the brightness and width values per line segment. For the brightness I need to transform the azimuth values to values between 0 (black for lines in the shadow) to 100 (white for lines fully exposed to the sun). Line widths range from 1 for lines in the sun or shadow) to 0.2 for lines in orthogonal direction.

The values will be written respectively to the columns color and width in the attribute table.

db.execute sql="UPDATE contours SET colour = cast(round(azimuth / 180 *100,0) as integer)"
db.execute sql="UPDATE contours SET width = (azimuth-90)/90*-0.8+0.2  WHERE azimuth <=90"
db.execute sql="UPDATE contours SET width = (90-azimuth)/90*-0.8+0.2  WHERE azimuth >90 AND azimuth <=180"

Export the vector layer as shapefile

As last step, export the vector layer as shapefile (and for good order remove the temporary vector layers).

v.out.ogr input=contours output=/home/paulo format=ESRI_Shapefile  --overwrite
g.remove type=vector pattern=tmp* -f

Create the map in QGIS

Open the DEM raster layer and the shapefile you just exported from GRASS in QGIS. Use the expression based symbology options to define the brightness based on the values in the column colour of the attribute table.

To dynamically define the colour of the line segments, based on their values in the “colour” column of the contour attribute table use an expression as in the screenshot above.

Likewise, use the values in the column width of the attribute table to define the width of the contour lines.

The line width is simply defined by the values in the column “width” of the contour attribute table. Or multiply or divide that value for an overall increase or decrease of the line width.

Finally, use the Overlay layer blending mode to blend contours and DEM colour.

The result… a beautiful altitude map. I am sure things can be done more efficiently, but it’ll gives you an idea of the possibilities.

Enjoy 🙂