Creating Aspect-slope Maps with Python

Introduction

Aspect-slope maps are an interesting way to visualize reliefs. In essence, these maps are created by first discretizing both metrics, slope and aspect, and secondly combining the resulting values into a single one. In a third step, a custom color code is assigned to each computed number. This blog post shows how to create such a map using Python. Figure 1 shows the elevation data (left) used in this blog posts can be converted to the respective aspect-slope map (right).

To a large degree, the discussed approach has been taken from here and to a smaller degree from here. Only minor modifications have been made to achieve more desirable effects.

Creating aspect and slope

Before the map can be created, both, aspect and slope need to be derived from a digital elevation model (DEM). For this blog post, we are going to use a high-quality DEM from the SwissALTI3D dataset, which can be found here. More specifically, we are going to use a tile that covers the Jungfrau mountain located in Switzerland, which can be obtained via this link. Next, we need to derive the aspect and slope datasets and store them in separate files. The probably easiest way to achieve this is to use the respective QGIS functionality, which can be found in “Raster” → “Analysis” → “Aspect” and “Slope”, respectively. The default settings are just fine. It is important, however, to generate the values as degrees and not radians.

Discretizing slope and aspect

The first step is to load the generated raster datasets using gdal and keep it as numpy arrays as follows:

# load slope dataset
slope_path = "slope.tif"

slope_ds = gdal.Open(slope_path)
slope_arr = slope_ds.ReadAsArray()

# load aspect dataset
aspect_path = "aspect.tif"

aspect_ds = gdal.Open(aspect_path)
aspect_arr = aspect_ds.ReadAsArray()

Once loaded, we are assigning discrete numbers to different ranges of slope and aspect as follows, overwriting the original values:

# discretize slope values
slope_arr[(slope_arr >= 0) & (slope_arr <= 4.999)] = 0
slope_arr[(slope_arr >= 5) & (slope_arr <= 19.999)] = 2
slope_arr[(slope_arr >= 20) & (slope_arr <= 39.999)] = 4
slope_arr[(slope_arr >= 40)] = 6

# discretize aspect values
aspect_arr[(aspect_arr >= 0) & (aspect_arr <= 22.499)] = 10
aspect_arr[(aspect_arr >= 22.5) & (aspect_arr <= 67.499)] = 20
aspect_arr[(aspect_arr >= 67.5) & (aspect_arr <= 112.499)] = 30
aspect_arr[(aspect_arr >= 112.5) & (aspect_arr <= 157.499)] = 40
aspect_arr[(aspect_arr >= 157.5) & (aspect_arr <= 202.499)] = 50
aspect_arr[(aspect_arr >= 202.5) & (aspect_arr <= 247.499)] = 60
aspect_arr[(aspect_arr >= 247.5) & (aspect_arr <= 292.499)] = 70
aspect_arr[(aspect_arr >= 292.5) & (aspect_arr <= 337.499)] = 80
aspect_arr[(aspect_arr >= 337.5) & (aspect_arr <= 360.5)] = 10

Note that for slope, the resulting discrete values have single digits, while the results for aspect have two digits and always end with zero. Hence, the combination of these slope and aspect values via addition will always result in unique values. Combining these numbers can be achieved as follows:

compound_classification = slope_arr + aspect_arr

Colorization

Now that the compound values have been generated, the next step is to assign a unique color to each of them. We are using the color palette proposed here, with only minor modifications. The color palette is stored as a Python map:

color_map_rgb = {10 : (181 , 181 , 181),
12 : (159 , 192 , 133),
14 : (157 , 221 , 94),
16 : (154 , 251 , 12),
18 : (155 , 255 ,  0),

20 : (181 , 181 , 181),
22 : (114 , 168 , 144),
24 : (61 , 171 , 113),
26 : (0 , 173 , 67),
28 : (128 , 128 , 128),

30 : (181 , 181 , 181),
32 : (124 , 142 , 173),
34 : (80 , 120 , 182),
36 : (0 , 104 , 192),
38 : (128 , 128 , 128),

40 : (181 , 181 , 181),
42 : (140 , 117 , 160),
44 : (119 , 71 , 157),
46 : (108 , 0 , 163),
48 : (128 , 128 , 128),

50 : (181 , 181 , 181),
52 : (180 , 123 , 161),
54 : (192 , 77 , 156),
56 : (202 , 0 , 156),
58 : (128 , 128 , 128),

60 : (181 , 181 , 181),
62 : (203 , 139 , 143),
64 : (231 , 111 , 122),
66 : (255 , 85 , 104),
68 : (128 , 128 , 128),

70 : (181 , 181 , 181),
72 : (197 , 165 , 138),
74 : (226 , 166 , 108),
76 : (255 , 171 , 71),
78 : (128 , 128 , 128),

80 : (181 , 181 , 181),
82 : (189 , 191 , 137),
84 : (214 , 219 , 94),
86 : (244 , 250 , 0),
88 : (128 , 128 , 128)}

Finally, we create a new numpy array compound_rgb using the dimensions of the array compound_classification and initialize it with the value 128 to receive a grey color by default. Then we iterate over all the items of the color palette. In each iteration, we determine those cells of compound_classification that equal the key within the current iteration and assign its corresponding value (the color as an rgb tuple) to all these cells of compound_rgb. The code is as follows:

# create the empty array
compound_rgb = np.zeros(shape=(compound_classification.shape[0], compound_classification.shape[1], 3))
compound_rgb[:,:,:] = np.asarray([128])

# iterate over color palette and assign values
for v, rgb in color_map_rgb.items():
    where_v = (compound_classification == v)
    compound_rgb[:,:,0][where_v] = rgb[0]
    compound_rgb[:,:,1][where_v] = rgb[1]
    compound_rgb[:,:,2][where_v] = rgb[2]

Finally, the result can be written to disk:

# write the slope aspect map
combined_path = "slope_aspect_combined_rgb.tif"

outdataset = driver.Create(combined_path, slope_arr.shape[1], slope_arr.shape[0], 3, gdal.GDT_Byte, ['COMPRESS=LZW'])
outdataset.SetGeoTransform(transform)
outdataset.SetProjection(projection)

outdataset.GetRasterBand(1).WriteArray(compound_rgb[:,:,0])
outdataset.GetRasterBand(2).WriteArray(compound_rgb[:,:,1])
outdataset.GetRasterBand(3).WriteArray(compound_rgb[:,:,2])

outdataset.FlushCache()
outdataset = None

Conclusion

In this post, we have shown how Python can be used to create aspect-slope maps. This comprised downloading an existing digital elevation model (DEM) and the derivation of its slope and aspect information as separate raster files. The values of the the slope and aspect raster files are then discretized and combined to a single file. Distinct color values are then assigned to each cell based on its combined value. The resulting slope-aspect map is then written to disk.