Introduction
Relief shading has a long tradition in Cartography as a means to present the terrain for a map in an intuitive manner. As such, relief shadings are also an ongoing topic of research with recent advances going so far as to employ Neural Networks to mimic hand-drawn relief shadings. However, simpler analytical shading algorithms can commonly be found in Geographic Information Systems (GIS), such as QGIS and ArcGIS. Good explanations on this procedure can, for example, be found here and here. In this post, we will show how analytical relief shadings can be efficiently implemented in Python by making use of the equations described in these two sources. The focus will lie on the implementation with Numpy broadcasting, a technique that allows to efficiently process large amounts of data without explicitly writing loops.
Data
Like with our previous post on label generation for isolines, we are going to use a digital elevation model (DEM) originating from the Shuttle Radar Topography Mission (SRTM) that can be downloaded from CGIAR. In fact, we are going to use the exact same extent as in the previous post. As such, it covers Switzerland and has been projected to the Swiss coordinate reference system (CRS) CH1903/+LV95 (EPSG:2056). The DEM is shown in Figure 1.
Reading and Preparing the DEM
The de-factor standard way to access raster data is by making use of GDAL, the Geospatial Data Abstraction Library. GDAL can commonly be installed via PIP or via a management system such as Anaconda. Using GDAL, we can load a DEM as follows:
import numpy as np
import gdal
dem_path = "data/dem.tif"
dem_dataset = gdal.Open(dem_path)
dem = dem_dataset.ReadAsArray().astype(np.int64)
print(dem.shape)
geotransform = dem_dataset.GetGeoTransform()
cellsize = geotransform[1]After importing NumPy and GDAL (lines 1 and 2), we can open the file with gdal.Open (line 5) and can convert the raster values into a numpy array named dem. Numpy is the de-factor standard for operation on matrices in Python. Note that we are converting the data type of the DEM from uint8 (i.e., a single “byte” per raster cell) into int64. The larger data type prevents integer overflows that might occur in some of the calculations later on. The dimensions of the loaded matrix, i.e., the number of pixels along each axis, can be shown by printing the shape attribute of our dem variable (line 7). In this case, it shows 222 pixels for the x-axis and 332 pixels for the y-axis. Lastly, we require the resolution (i.e., cell size) of each raster pixel. This can be obtained by accessing the second entry of the geotransform. More information about the geotransform can be found in here.
Ultimately, relief shading produces a single brightness value in the range [0, 255] for each pixel of a DEM. To do so, algorithms take into account the neighborhood of each pixel. Typical names of cells can be found in Table 1 with the cell named “e” being in the center.
| a | b | c |
| d | e | f |
| g | h | i |
We will make use of the broadcasting functionality of NumPy to efficiently calculate the different steps of relief shading. Broadcasting is a method that allows common to perform arithmetics on high-dimensional arrays without looping over each value individually. The main requirement for this method to work is that the dimensions of the arrays to process are compatible. Usually, this means that their dimensions need to match. Details on the different rules for broadcasting to apply can be found here.
For the first step we need to take for applying braodcasting, thus is to make sure that the dimensions of the arrays match. We can use array slicing to obtain the distinct arrays each holding the values for each of the cell names shown in Table 1:
a = dem[:-2,:-2]
b = dem[:-2,1:-1]
c = dem[:-2,2:]
d = dem[1:-1,:-2]
e = dem[1:-1,1:-1]
f = dem[1:-1,2:]
g = dem[2:,:-2]
h = dem[2:,1:-1]
i = dem[2:,2:]For example, to obtain the cells that correspond to the name h, we need to crop the first two rows (indicated by “2:”) of the array as well as the first column and the last column (as indicated by “1:-1”). Note that negative numbers count the axis from the end of that dimension. For more information about how array slicing works, the reader is referred to this tutorial.
To provide a numerical example, suppose we are using the following artificially generated dumm_dem. Applying the slices as described before will yield the desired subsets:
dummy_dem = np.asarray([[ 1, 2, 3, 4],
[ 5, 6, 7, 8],
[ 9,10,11,12],
[13,14,15,16]])
a = dummy_dem[:-2,:-2] # a
# result:
# array([[1, 2],
# [5, 6]])
b = dummy_dem[:-2,1:-1] # b
# result:
# array([[2, 3],
# [6, 7]])
c = dummy_dem[:-2,2:] # c
# result:
# array([[3, 4],
# [7, 8]])
d = dummy_dem[1:-1,:-2] # d
# result:
# array([[ 5, 6],
# [ 9, 10]])
e = dummy_dem[1:-1,1:-1] # e
# result:
# array([[ 6, 7],
# [10, 11]])
f = dummy_dem[1:-1,2:] # f
# result:
# array([[ 7, 8],
# [11, 12]])
g = dummy_dem[2:,:-2] # g
# result:
# array([[ 9, 10],
# [13, 14]])
h = dummy_dem[2:,1:-1] # h
# result:
# array([[10, 11],
# [14, 15]])
i = dummy_dem[2:,2:] # i
# result:
# array([[11, 12],
# [15, 16]])Note that all the results have the same dimensions, exactly what we need. Now, we can apply broadcasting to do the first calculations.
Shading Calculation
The step for carrying out analytical relief shading is to calculate the differentials in the x- and y-directions of the DEM. The equations we are using for this purpose are as follows:
# calculate derivatives
dzdx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize)
dzdy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize)
Because the variables are numpy arrays, broadcasting is automatically applied and the result is, again, a numpy array with the same dimensions as the inputs. Note that we are making use of the cellsize variable that we obtained form the input raster. Likewise, we can calculate slope and aspect as follows:
exeggeration = 10
slope = np.arctan(exeggeration * np.sqrt(dzdx ** 2 + dzdy ** 2))
aspect = np.arctan2(dzdy, -dzdx)Note that we can apply an exaggeration parameter that adjusts the height of the relief, thereby increasing the effect if a value greater than 1 is used or lowering its effect if the value is smaller than 1. A value of exactly 1 uses the original height.
According to this explanation aspect needs to be further adapted. In particular, the following modifications need to be made:
- if dzdx is not zero and aspect is smaller than zero, then add 2 * pi to aspect
- if dzdx is zero and dzdy is greater than zero, set aspect to pi / 2
- if dzdx is zero and dzdy is smaller than zero, set aspect to 2 * math.pi – math.pi / 2
To apply these modifications in an efficient way, we can make use of boolean indexing. Boolean indexing can be used to only select those cells in an array for which a condition is true. More details can be found in this blog post. This means, that we first need to define the boolean masks to use for indexing. Once these are defined, they can be combined and used fo index and modify aspect:
# define masks
dzdx_zero = dzdx == 0
dzdx_nonzero = ~dzdx_zero
aspect_smaller_zero = aspect < 0
dzdy_greater_zero = dzdy > 0
dzdy_smaller_zero = dzdy < 0
# use masks for indexing
aspect[dzdx_nonzero & aspect_smaller_zero] = 2 * math.pi + aspect[dzdx_nonzero & aspect_smaller_zero]
aspect[dzdx_zero & dzdy_greater_zero] = math.pi / 2
aspect[dzdx_zero & dzdy_smaller_zero] = 2 * math.pi - math.pi / 2Note how we are combining different conditions with the & operator and make sure that the dimensions of the arrays match (line 11). Next, we need to define the altitude and the azimuth of the “imaginary” location where the light source is placed. This is defined in terms of its altitude and its azimuth. The values of 45° and 315° are typically used conventions. This makes the light shine from the “upper left” corner of the DEM:
# define altitude
altitude = 45
zenith_deg = 90 - altitude
zenith_rad = math.radians(zenith_deg)
# define azimuth
azimuth = 315
azimuth_math = 360 - azimuth + 90
if azimuth_math >= 360:
azimuth_math = azimuth_math - 360
azimuth_rad = math.radians(azimuth_math)Finally, we can put everything together to actually calculate the shaded value for each cell. Note that we are making use of the NumPy routines to calculate the trigonometric operations as these allow for broadcasting rather than using the equivalents of Python’s math library:
# calculate reflectance
shaded = 255 * (np.cos(zenith_rad) * np.cos(slope) + (np.sin(zenith_rad) * np.sin(slope) * np.cos(azimuth_rad - aspect)))The multiplication with 255 is necessary to obtain the grayscale value representing the shading at this location. Finally, we need to determine what happens with the border pixels, i.e., those locations that do not have enough neighbors to calculate the shading. In this case, we simply copy the values of the next neighboring pixels for which a shading value could be calculated. Again, we can use broadcasting for this purpose:
# create an array that is one pixel larger than the origional shading
shaded_padded = np.zeros((shaded.shape[0] + 2, shaded.shape[1] + 2))
# fill the inner part with the created shaded relief
shaded_padded[1:-1,1:-1] = shaded
# fill the borders
shaded_padded[0,1:-1] = shaded[0,:]
shaded_padded[-1,1:-1] = shaded[-1,:]
shaded_padded[1:-1,0] = shaded[:,0]
shaded_padded[1:-1,-1] = shaded[:,-1]
# fill the corners
shaded_padded[0,0] = shaded[0,0]
shaded_padded[0,-1] = shaded[0,-1]
shaded_padded[-1,0] = shaded[-1,0]
shaded_padded[-1,-1] = shaded[-1,-1]Writing the Shaded Relief
Again, we can make use of GDAL to write our result to disk. As our original DEM was georeferenced, we can make use of those information to likewise georeference the shaded relief:
# get projection information
transform = dem_dataset.GetGeoTransform()
projection = dem_dataset.GetProjection()
# create output dataset
driver = gdal.GetDriverByName("GTiff")
shaded_path = "data/shaded_relief.tif"
shaded_dataset = driver.Create(shaded_path, dem.shape[1], dem.shape[0], 1, gdal.GDT_Byte, ['COMPRESS=LZW'])
# copy transform and projection
shaded_dataset.SetGeoTransform(transform)
shaded_dataset.SetProjection(projection)
# write the shaded relief to disk
shaded_dataset.GetRasterBand(1).WriteArray(shaded_padded)
# make sure that all the data is written
shaded_dataset.FlushCache()
# close the dataset
shaded_dataset = NoneNote that we are using the transform and projection information from the original DEM so that they are placed exactly at the same location (lines 2, 3 and lines 14, 15). The format will be GeoTiff (line 6) and there will only be a single band of type byte (uint8) that is compressed via the LZW method (line 9). By flushing the cache (line 19) and setting the created dataset to None, we make sure that all data is written to the file and the dataset is closed (line 22). Finally, we can visualize the result:
The whole script can be found here:
# author: Magnus Heitzler
# website: heitzler-geoinformatik.de
import numpy as np
import gdal, math
dem_path = "data/dem.tif"
dem_dataset = gdal.Open(dem_path)
dem = dem_dataset.ReadAsArray().astype(np.int64)
print(dem.shape)
geotransform = dem_dataset.GetGeoTransform()
cellsize = geotransform[1]
# get the subsets
a = dem[:-2,:-2]
b = dem[:-2,1:-1]
c = dem[:-2,2:]
d = dem[1:-1,:-2]
e = dem[1:-1,1:-1]
f = dem[1:-1,2:]
g = dem[2:,:-2]
h = dem[2:,1:-1]
i = dem[2:,2:]
# calculate the derivatives
dzdx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * cellsize)
dzdy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * cellsize)
# calculate the slope
exeggeration = 10
slope = np.arctan(exeggeration * np.sqrt(dzdx ** 2 + dzdy ** 2))
# calculate the aspect
aspect = np.arctan2(dzdy, -dzdx)
# define masks
dzdx_zero = dzdx == 0
dzdx_nonzero = ~dzdx_zero
aspect_smaller_zero = aspect < 0
dzdy_greater_zero = dzdy > 0
dzdy_smaller_zero = dzdy < 0
# adapt aspect
aspect[dzdx_nonzero & aspect_smaller_zero] = 2 * math.pi + aspect[dzdx_nonzero & aspect_smaller_zero]
aspect[dzdx_zero & dzdy_greater_zero] = math.pi / 2
aspect[dzdx_zero & dzdy_smaller_zero] = 2 * math.pi - math.pi / 2
# calculate reflectance
altitude = 45
zenith_deg = 90 - altitude
zenith_rad = math.radians(zenith_deg)
azimuth = 315
azimuth_math = 360 - azimuth + 90
if azimuth_math >= 360:
azimuth_math = azimuth_math - 360
azimuth_rad = math.radians(azimuth_math)
shaded = 255 * (np.cos(zenith_rad) * np.cos(slope) + (np.sin(zenith_rad) * np.sin(slope) * np.cos(azimuth_rad - aspect)))
# add padding
shaded_padded = np.zeros((shaded.shape[0] + 2, shaded.shape[1] + 2))
# fill the inner part with the created shaded relief
shaded_padded[1:-1,1:-1] = shaded
# fill the borders
shaded_padded[0,1:-1] = shaded[0,:]
shaded_padded[-1,1:-1] = shaded[-1,:]
shaded_padded[1:-1,0] = shaded[:,0]
shaded_padded[1:-1,-1] = shaded[:,-1]
# fill the corners
shaded_padded[0,0] = shaded[0,0]
shaded_padded[0,-1] = shaded[0,-1]
shaded_padded[-1,0] = shaded[-1,0]
shaded_padded[-1,-1] = shaded[-1,-1]
# get projection information
projection = dem_dataset.GetProjection()
# create output dataset
driver = gdal.GetDriverByName("GTiff")
shaded_path = "data/shaded_relief.tif"
shaded_dataset = driver.Create(shaded_path, dem.shape[1], dem.shape[0], 1, gdal.GDT_Byte, ['COMPRESS=LZW'])
# copy transform and projection
shaded_dataset.SetGeoTransform(transform)
shaded_dataset.SetProjection(projection)
# write the shaded relief to disk
shaded_dataset.GetRasterBand(1).WriteArray(shaded_padded)
# make sure that all the data is written
shaded_dataset.FlushCache()
# close the dataset
shaded_dataset = NoneConclusion
Matrix operations can be efficiently executed using Python when employing the Numpy library. Yet, the syntax certainly requires some training until it can be effectively used. For GIS, this implies that map algebra of all sorts can be easily implemented using boiler-plate Python libraries. Relief shading as demonstrated in this blog post is only one example for this. It should also be noted that the implementations of relief shadings slightly differ from source to source and smaller deviations can commonly be found. In addition, nowadays other shading methods, like those based on Neural Networks, might be better suited, depending on the use case.