How to create illuminated contours, Tanaka-style

In the category “last night on Twitter”, a challenge I couldn’t resist: creating illuminated contours (aka Tanaka contours) in QGIS. Daniel P. Huffman started the thread by posting this great example:

CFnWnA5UkAAuFm9

This was quickly picked up by Hannes Kröger who blogged about his first attempt at reproducing the effect using QGIS and GIMP. Obviously, that left the challenge of finding a QGIS-only solution.

Everything that’s needed to create this effect is a DEM. As Hannes describes in his post, the DEM can then be used to compute the contour lines, e.g. with Raster | Extraction | Contour:

gdal_contour -a ELEV -i 100.0 C:\Users\anita\Geodata\misc\mt-st-helens\10.2.1.1043901.dem C:/Users/anita/Geodata/misc/mt-st-helens/countours

Screenshot 2015-05-24 11.17.49

contours

In order to be able to compute the brightness of the illuminated contours, we need to compute the orientation of every subsection of the contours. Therefore, we need to split the contour lines at each node. One way to do this is using v.split from the Processing toolbox:

Screenshot 2015-05-24 11.23.11

When we split the contours and visualize the result using arrows, we can see that they all wrap around the mountain in clockwise direction (light DEM cells equal higher elevation):

split_contours

After the split, we can compute the orientation of the contour subsections using, for example, a user-defined function:

Screenshot 2015-05-24 19.09.12

"""
Define new functions using @qgsfunction. feature and parent must always be the
last args. Use args=-1 to pass a list of values as arguments
"""

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def azimuth(x1,y1,x2,y2, feature, parent):
	p1 = QgsPoint(x1,y1)
	p2 = QgsPoint(x2,y2)
	a = p1.azimuth(p2)
	if a < 0:

		a += 360
	return a

This function can then be used in a Field calculator expression:

Screenshot 2015-05-24 19.11.53

Based on the orientation, we can then write an expression to control the contour line color. For example, if we want the sun to appear in the north west (-45°) we can use:

color_hsl( 0,0,
  scale_linear( abs(
    ( CASE WHEN &quot;azimuth&quot;-45 &lt; 0
      THEN &quot;azimuth&quot;-45+360
      ELSE &quot;azimuth&quot;-45
    END )
  -180), 0, 180, 0, 100)
  )

This will color the lines which are directly exposed to the sun white hsl(0,0,100) while the ones in the shadows will be black hsl(0,0,0).

Screenshot 2015-05-24 11.55.50

Use the Overlay layer blending mode to blend contours and DEM color:

illuminated_contours

The final step, to get as close to the original design as possible, is to create the effect of discrete elevation classes instead of a smooth color gradient. This can easily be achieved by changing the color interpolation mode of the DEM from Linear to Discrete:

Screenshot 2015-05-24 12.11.01

This leaves us with the following gorgeous effect:

tanaka_contours

As Hannes pointed out, another important aspect of Tanaka’s method is to also alter the contour line width. Lines in the sun or shadow should be wider (1 in this example) than those in orthogonal direction (0.2 in this example):

scale_linear(
abs( abs(
  ( CASE WHEN &quot;azimuth&quot;-45 &lt; 0
    THEN  &quot;azimuth&quot;-45+360
    ELSE  &quot;azimuth&quot;-45
  END )
-180) -90),
0, 90, 0.2, 1)

datadefined_line_width

Enjoy!

Advertisements
29 comments
  1. Terry said:

    I had a problem calculating azimuth. I kept getting a spread of values from 180 – 360 with a bunch of null values for anything below 180. I fixed it by changing a couple of lines:

    from qgis.core import*
    from qgis.gui import*

    @qgsfunction(args=”auto”, group=’Custom’)
    def azimuth(x1,y1,x2,y2,feature,parent):
    p1 = QgsPoint(x1,y1)
    p2 = QgsPoint(x2,y2)
    a = p1.azimuth(p2)
    if a < 181:
    a += 180
    return a

    Then, the shading didn't work properly, I assume due to the changes mentioned above. This was easily resolved by changing 45 degrees to 225 degrees:

    color_hsl( 0,0,
    scale_linear( abs(
    ( CASE WHEN "azimuth"-225 < 0
    THEN "azimuth"-225+360
    ELSE "azimuth"-225
    END )
    -180), 0, 180, 0, 100)
    )

    In the end, though, I got a beautiful map. Thanks so much for this.

    • Thanks a lot for your feedback! Would be interesting to find out why QgsPoint.azimuth() seems to work differently for you …

      • Terry said:

        I’m using 2.8.1 instead of the dev build. Could that have something to do with it?

      • Maybe though it shouldn’t. For me, 2.8.2 and dev produce the same results. Maybe this uses geos (or something similar) and we have different versions.

  2. that looks awesome. would be cool to create such a theme in openlayers although the on the fly processing for the style might be to process intensive.

    • True, I don’t think trying to render this in OpenLayers would end well. But one could prepare some nice tiles.

  3. Nice tutorial.Thanks.

  4. Ross said:

    Very cool. It reminds me of the default OpenSUSE wallpaper (who knows why they choose to default to that).

  5. Miroslav Umlauf said:

    On the first look it is beautiful but… It is dependent on the contour starting and ending point and if you have more mountains and up and downs it ends up with making lighter contours on downhills and dark on uphills. All is needed is the azimuth actually needs to be calculated from slope direction, not from contour direction.

    • If you care to share details about the issue you discovered, a picture would be useful. Of course the effect depends on the contour start and end point. Therefore it heavily depends on the consistency of contour line orientation generated by the contour tool.

      • Miroslav Umlauf said:

        Thank you. I found out that yesterday when trying to figure out why all contours in your images are consistent and mine not. I have used matplotlib.pyplot.contour() to generate mine. Well if I will need style dependent on direction now I know I need to use gdal_contour. Thank you again.

  6. volaya said:

    Quick suggestion:

    This can all be automated and put into a graphical model quite easily, and i guess it would make a good exercise.

    You can use the “Advanced python calculator” for the azimuth calculation, and the “Set style for vector layer” algorithm for setting the style (which you should create before and save to a known location)

    For systems were GRASS is not installed, the “Explode lines” algorithm im Propcessing does the same thing that you are doing with v.split.vert

    Great work!

    • Good point! Might be a good follow-up post.

  7. The Tanaka contours reference article that you cite, http://wiki.gis.com/wiki/index.php/Tanaka_contours, says that the contour lines perpendicular to the light source are drawn thinner, but the lines in your image are all the same size. Also the lines do not “change from white in the northwest to black in the southeast”. Is this correct, and if so, is this image really ‘Tanaka’?

    Trying to figure out what the image means, it seems to me that:
    – colour from blue to white on the image (but not the contour lines) represents elevation
    – colour on the contour lines represents both slope AND aspect. For example lines on the right side of the cone are white because they are perpendicular to the sun and because their aspect is between 270 and 360 degrees, that is, they slope TOWARD the sun; lines in the lower right side of the image that are also perpendicular to the sun are dark because their aspect is between 90 and 180 degrees, that is because they slope AWAY FROM the sun. Is this correct?

    Great posting, thanks for everything.

    • The last paragraph covers the varying line width topic.

      The line color simply shows contour line orientation. That way it represents the aspect, i.e. “the compass direction that a slope faces”. The slope is encoded in the distance between contour lines I’d say.

  8. Ben C said:

    Who is Tanaka? Can we give Tanaka credit for the original idea?

  9. I have problem in the field calculator step that didn’t crop up originally. It says I have the wrong number of parenthesis when I type in: azimuth(xat(0),yat(0),xat(-1),yat(-1)). I’ve tried different data, restarting the program and I’m just about reinstall qgis. Any suggestions for other troubleshooting because I’m running out of ideas. or possible another route to getting the azimuth field.

    • Which QGIS version do you use? I don’t think reinstalling will do any good.

      • stephen engley said:

        woops, sorry, its early. yeah, pisa and then lyon. Thanks for the tutorial, while it worked it was brilliant.

      • Some function names changed in Lyon I think. Will check later.

  10. stephen engley said:

    I was using pisa and then lyon

    • Sorry for the late reply. Everything seems to work fine here in QGIS 2.13. Did you create the azimuth function in the Function Editor before trying to use it in the Field Calculator?

      • stephen engley said:

        It was a coding error, no matter how many times I types in the code I typed in the same error. I put it down for a couple weeks and then I saw the error. Thanks for checking back in.

%d bloggers like this: