Hillshading is the art and science of making a graphical representation of an illuminated landscape. It has a long history that includes names as Leonardo da Vinci (with the chiaroscuro technique) and many famous and not so famous Swiss and Japanese cartographers (have a look at some examples of mastery on manual hillshading by Eduard Imhof). A very complete account can be found in the excellent and very enjoyable paper “Hill shading and the reflectance map”, by Berthold Horn [1].

Shading has an intuitive appeal because it conveys information on shape in a straightforward manner. This is best envisaged as the difference in perception that a person experiences during a cloudy or a sunny day. In a cloudy day light is diffuse and does not produce shadows, giving the impression of flatness, whereas in a sunny day shadows allow us to make an immediate mental model of the shapes and enables better estimation of distances. This is very relevant in many practical situations, in particular [1]

… when the interpreter’s time is limited, as in aviation, for users that are not trained cartographers, and for small scale maps, where contours degenerate into messy tangles of lines.

Most importantly, the techniques involved in computing shaded maps are fundamental for applying topographic correction to satellite images. They can also become important when it is relevant to convey information about topography as context to something else, but when a exact quantitative specification is not necessary nor desirable, as in tourist maps or as background information in geoscientific studies, as shown below.

Every modern GIS package provides a convenient way to compile a shaded relief map, often including very involved corrections. We here focus on the fundamental aspects, leaving refinements aside, and will confine ourselves to the important and sufficiently general case of Lambertian surfaces, using as an example the DEM we compiled in previous posts of the region around Vesuvius volcano, in Italy.

The central concept here is to understand that the irradiance, the amount of light coming upon a surface element, depends on the position of the light source relative to the surface element. The radiance, the amount of light being reflected by the surface element, depends on the irradiance but also on the reflectance properties of the surface. As we mentioned, we here focus on Lambertian surfaces, reflecting diffusively and equally on all directions, in which case we can safely assume that radiance will be proportional to irradiance. For natural surfaces illuminated by the Sun the problem of finding the irradiance, or solar incidence, reduces to an application of the spherical law of cosines:

$cos \gamma_{i} = cos \theta_{p} cos \theta_{z} + sin \theta_{p} sin \theta_{z} cos (\phi_{a} - \phi_{o})$

where $\gamma_{i}$ is the solar incidence, $\theta_{p}$ is the slope of the surface, $\theta_{z}$ is the zenith angle, or the angle between the Sun position and the vertical, $\phi_{a}$ is the the Sun azimuth from the North, and $\phi_{o}$ is the surface aspect, that is, the direction in which the surface is sloping.

We begin by calculating the slope and aspect from the DEM:

1. slope:

$\theta_{p} = atan ( \sqrt{ (\frac{dz}{dx})^{2} + (\frac{dz}{dy})^{2} })$

we can easily apply this manipulation using grdmath, the GMT reverse polish notation calculator, on our DEM grid, here denoted dem.grd:

grdmath dem.grd DDX dem.grd DDY R2 SQRT ATAN R2D = slope.grd

Note that the function R2 is equal to the sum of the squares of the first two terms, and R2D is the Radian To Degree converter.

2. aspect:

$\phi_{o} = atan (\frac{dx}{dy})$

and similarly:

grdmath dem.grd DDX dem.grd DDY ATAN2 R2D 180 ADD = aspect.grd

Note that we used the function ATAN2, which specifically accounts for the signs of the numerator and denominator and thus gives the angles in the right quadrant. We add 180 degrees to get angles from 0 to 360, instead in the range -180 – +180.

Our partial results look like this:

We can now turn into computing cosine of the irradiance. As Sun azimuth and elevation we choose 315° and 45°, that is: illumination from the northwest at half sky, following the cartographic convention for shaded maps:

We do:

grdmath slope.grd D2R COS 45 D2R COS MUL slope.grd D2R SIN 45 D2R SIN MUL 315 D2R aspect.grd D2R SUB COS MUL ADD = cosgammai_315.grd

(this can also be accomplished using GMT’s grdgradient, which also has a function for specifying the properties of the reflective surface)

Our results look like this:

Presto! This is an image upon which we could drape a photograph or satellite image, or even draw the contours. This is also a fundamental piece of information needed for topographical correction of images, a topic of a previous post.

[1] Horn, B. K. (1981). Hill shading and the reflectance map. Proceedings of the IEEE, 69(1), 14-47. Find it in ResearchGate

[2] De Angelis, H. & Kleman, J. (2008). Palaeo‐ice‐stream onsets: examples from the north‐eastern Laurentide Ice Sheet. Earth Surface Processes and Landforms, 33(4), 560-572. Find it in ResearchGate.