Author: Paulo van Breugel
Updated on: 2015-05-25
In this tutorial, I want to explore how to create illuminated or Tanaka contours, 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.
Using the newly created layers, I'll use QGIS to create the final map, 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.
Likewise, use the values in the column width of the attribute table to define the width of the contour lines.
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.
Software: @grassgis @qgis
Tools: @r_contour @r_split @v_out_ogr