Combining multiple shaded reliefs with hypsometric colouring and digital orthophotos using free and open-source software

: In this paper we combined layers created from several terrain rendering techniques, namely a shaded relief rendered in the free and open-source 3D computer graphics software Blender, a hillshade created in the free and open-source Geographic Information System (GIS) software QGIS, a hypsometric coloured Digital Elevation Model (DEM) and a draped digital orthophoto. Following a recent trend in the cartographic community towards using Blender, we tried to improve the standard relief visualization in common GIS software by blending it with a shaded relief rendered in Blender. Using different QGIS blending modes and opacity values we achieved different graphic visualizations. To compare and evaluate the suitability of different rendering techniques we chose national park Risnjak located in Croatia because of its specific and diverse terrain landforms. After comparing different input layers and parameter sets, we selected the blending combination which is best suited for visualizing terrain characteristics of all Croatian national parks. The result is a shaded relief created for every national park which is combined from a shaded relief rendered in Blender, a hillshade created in QGIS, a hypsometric coloured DEM and a draped digital orthophoto.


Introduction
Hillshading, or creating shaded relief, is an effective method for visualizing the realistic three-dimensional shape of the earth on a two-dimensional map by using illumination from one or more light sources. When illuminating the surface of an object a shadow effect is produced which can be used in maps to visually enhance different topographic forms of Earth's surface. The computer-based process of generating shaded relief from a digital elevation model (DEM) is called analytical relief shading (URL1). The main design principles for creating effective shaded relief date back to the time when it was manually drawn (Imhof, 1982). Today shaded relief is almost exclusively computer-generated using either GIS software such as ArcGIS, QGIS or GRASS GISS, or graphic and 3D editors such as Adobe Photoshop or Blender. Although software for creating shaded relief from DEMs is around for several decades, not until recently blending multiple layers to render a shaded relief has not been possible using GIS software. Also, there is a recent trend in the cartographic community towards using Blender (URL 2) for creating renders of shaded relief. Several famous cartography blogs like Daniel Huffman's somethingaboutmaps (URL 3), Morgan Hite's The Wandering Cartographer (URL 4) and David Garcia's Mapmaker (URL 5) introduced relief shading tutorials using Blender. They have shown that shaded reliefs rendered in Blender create a more natural, softer, and attractive shaded relief with clearer structures thanks to the improved modelling of lighting in a threedimensional workspace (URL 3). Besides this recent trend, our research was also inspired by a poster of Britain's National Parks created by Joe Harrison from the GeoDataViz team of the Ordnance Survey (URL 6). The poster is a beautiful showcase of a compiled hillshade created in Blender with a natural colour palette terrain style used in traditional maps. In this paper we compared different visualizations of shaded relief created by combining multiple shaded reliefs with hypsometric colouring and digital orthophotos. The main goal was to visually enhance analytical relief shading in QGIS by blending it together with a shaded relief created in Blender. Additionally, we have used hypsometric tinting, also called hypsometric colouring, elevation colouring or layer "painting", which has proven to be a very effective method for portraying the landscape (Jenny and Hurni, 2006). By adding a colour palette to the blended relief, a visually improved depiction of different elevation levels within every single national park was achieved. Furthermore, we applied the same colour palette with the same values for pixel minimum and maximum to achieve visually comparable elevation levels between each national park. Finally, we have draped a digital orthophoto over as well as and under the blended relief. The basic idea was to investigate if the texture of the orthophoto will give a more natural three-dimensional map (Nighbert, 2000). The software used in this relief modelling process is free and open source and thus available to everyone. Data used in this research is open data and hence publicly available. The result is a shaded relief created for every national park which is combined from a shaded relief rendered in Blender, a hillshade created in QGIS, a hypsometric coloured DEM and a draped digital orthophoto by using different blending modes and opacity values.

Study area
For the study area we have decided to use Croatian national parks similar to the Ordnance Survey poster. The reason for this is the variety of terrain types available in all the 8 national parks. Figure 1 shows an overview map of the location of the parks with their according number labelled. Three of them are located on islands. National park (NP) Mljet (1) is part of the eponymous island, whereas the other two are island archipelagos. Brijuni (7) NP consists of 14 islands and Kornati (2) of around 140 islands. NP Risnjak (8) is located in a mountainous and heavily forested region. Paklenica (4) is a karts and steep river canyon stretching from the see to the southern slopes of the Velebit mountain, whereas North Velebit (6) covers the northern part of the same mountain. Plitvice (5) is the largest NP and like Krka (3) stretches in a karts hilly area with several large cascading waterfalls.

Data
The main data source for this research is a DEM. For our study area we used the Digital Elevation Model over Europe (EU-DEM) which is a publicly available (open) data source of the Copernicus programme, managed by the European Environment Agency (URL 7). The EU-DEM is in fact a digital surface model (DSM) created by fusing the global Space Shuttle Radar Topography Mission (SRTM) DEM and ASTER Global DEM (ASTER GDEM). It has a 25 m spatial resolution and a stated ± 7 m vertical accuracy. EU-DEM data covers the European territory and comes in the ETRS89 Lambert Azimuthal Equal Area coordinate reference system (EPSG code: 3035). For our study area we only had to download one EU-DEM tile, namely the E40N20. The digital orthophoto (DO) image was downloaded from the Web Mapping Service (WMS) of the State Geodetic Administration (SGA) of the Republic of Croatia. For every national park, a GeoTIFF raster was exported with a spatial resolution of 5 m per pixel together with a corresponding Tiff World File (data of the SGA, collected on November 12 th 2020, DOF5 from 2017 and 2018).
For the boundaries of the NPs, we used publicly available data from the Web Feature Service (WFS) of the Bioportal (URL 8), which is part of the Nature Protection Information System managed by the Institute for Environmental Protection and Nature of the Ministry of Economy and Sustainable Development of the Republic of Croatia. More specifically, for the NPs in the inland we used the layer of boundaries of protected areas and for the NPs located on islands the layer of habitat mapping which includes the seashore line.

Data processing
First, we had to extract the exact boundaries of our study areas i.e., the NPs. Polygons for inland NPs were exported easily using QGIS from the WFS layer of boundaries of protected areas. For NPs located on islands this layer does not contain coastal borders vectorized at the seashore line, rather a polygon surrounding all islands including the see. Therefore, to create the exact boundary we used the WFS layer containing the seashore line. After deleting all unnecessary lines, the polygons for NPs Brijuni and Kornati were created by converting lines to polygons and exporting them into individual files. For NP Mljet we had to check the boundary of the NP going through the inland of the island and after we made sure it was closed, we converted it to a polygon and exported. After creating the polygons, we used them to clip the DEM and DO to the exact boundaries of each NP. Finally, the DEM had to be rescaled i.e., the elevation range was stretched to cover values between 0 and 65535 and saved as a 16-bit Unsigned Integer TIFF since Blender reads only elevation data rounded to integers.

Rendering shaded relief in Blender
To import the processed DEM files to Blender we used the BlenderGIS addon (URL 9). Out of our DEMs a plane mesh was build using 'DEM raw data build (slow)' mode. After import, a grey 3D terrain render is floating in Blenders coordinate space ( Figure 2). The next few settings were set in Blender only for the first render and were subsequently used to create other renders for all NPs. First, we adjusted the vertical exaggeration of the plane which is in Blender called 'Scale Z'. We set it to 0.5 because our plane had stretched elevations and therefore, we had to lower the vertical exaggeration. Then we changed the Blender rendering engine from Eevee to Cycles. Both rendering engines have advantages and disadvantages, but we used the later because it is commonly used in all relief shading tutorials in Blender (see Introduction). In the next step we set up a new material for our plane. Materials are used in 3D models to simulate variations in surface colour, surface texture, or roughness. By reflecting, absorbing, and scattering light Blender manages to give hillshades a more realistic look compared to GIS software (URL 10). After testing different surface options for material, we used the default Principled BSDF (Bidirectional Scattering Distribution Function). In the material settings we changed the base colour of the surface by setting the R, G and B values to 0.615 to make it darker. For NPs Krka and Plitvice, we set the R, G and B values to 0.532 to make the surface even darker. Reason for this is to create a more lightabsorbing material i.e., a more distinctive difference between shadows and highlights. Another setting common to all our renders is true displacement of the surface which was changed from 'Bump only' to 'Displacement only' to make the rendered image more realistic. Creating a new camera to position the orthographic view to capture the whole plane was easily accomplished using the Georef camera, a functionality built in BlenderGIS addon. A new camera had to be created for every render because every DEM has different width and height. Beside the material of the surface, the possibility to add multiple light sources and different light colours is what distinguishes relief shading in Blender from GIS software. Not only multiple light sources but the strength, colour, type, and position of the light source affect the scene's lighting, and thus the final render. When creating a shaded relief in GIS software usually there are three key parameters: vertical exaggeration, elevation of the sun and direction (azimuth) of the sun. Although Blender has different types of light sources (point, sun, spot, area) we want to simulate the Sun as the only light source for relief shading. For every render of a NP, we adjusted the sun's strength depending on to the terrain steepness and height differences (Table 1). By changing the strength, we decide how bright the light source is. A higher number means a brighter sunlight hence a lighter relief. The elevation angle of the sun was adjusted by changing the Y rotation of the light source. Important to know is that the rotation angle is the complement of the elevation angle, so if we want the sun to be 30° above the horizon, we had to set the Y rotation to 60°. To change the azimuth of the sun we adjusted the Z rotation of the light source. The Z rotation begins at 0° for East and increases counterclockwise which means if we want the sun coming from Northwest, we had to set the Z rotation to 135°. Table 1 shows the different settings we chose to create our final renders. Finally, we had to set the X and Y resolution of every NP render individually to match the pixel dimensions of our input DEMs. In this way we could georeference every single hillshade by applying the world file of the corresponding input DEM. Blender outputs were exported as TIFF files with 300 rendering samples for each pixel.

Blending layers together
Graphic artist and cartographers are widely using blending (compositing) techniques to fine tune how layers are blended together as well as how they interact with layers below. Usually in GIS software layers are overlayed by controlling layer transparency. Blending modes were historically used only in graphic design software but fortunately QGIS also offers the possibility to blend various layers (Tzvetkov, 2018). In our research we have created 8 visualizations ( Figure  3) of NP Risnjak by blending together different raster layers created from the clipped EU-DEM, from the clipped DO and from the shaded relief created in Blender.
Blending was achieved only using blending modes and by adjusting global opacity of layers. Three layers were created from the EU-DEM: a hypsometrically coloured DEM, a one-directional hillshade, and a multi-directional hillshade. After comparing all visualizations and finding the best blending combination we created the final hillshades for all NPs.

Hypsometric colouring
In QGIS hypsometric colouring is achieved by changing the layer style of the DEM to singleband pseudocolor. Two options are available for colour tinting: discrete (a specific colour assigned to a specific elevation interval) and continuous. The latter creates a colour ramp with gradually varying colours. We used it to represent the full range of elevation values. After trying different colour palettes available in QGIS we concluded that they were not satisfactory for our visualization purpose. Therefore, we imported several colour gradients from an archive called 'cpt-city' (URL 11). Although some colour pallets of this archive are already available in QGIS we used a palette manually imported called 'tpgl hwf'. To show height differences within a NP but also between individual NPs, it is necessary to carefully select the elevation range that the colour ramp will cover. After trying different values, we selected a range from 0 to 2000 m.

Hillshading in QGIS
Both hillshades of Risnjak NP, the one-directional and the multi-directional, were created using QGIS built-in raster analysis tool. As the input layer we used the clipped EU-DEM. QGIS calculates DEM shading according to the azimuth (horizontal angle; clockwise direction) and the vertical angle (altitude in degrees) of the sun (URL 12). Vertical exaggeration can be adjusted by changing the Z factor. Our parameters for both hillshades were as follows: 45° altitude, 315° azimuth (NW direction), 1.5 Z factor.

Shaded reliefs of NPs
In order to create the final result i.e., a shaded relief for every national park, we had to analyse the visualizations of Risnjak NP shown in the previous chapter and choose the one that suits all the NPs terrain landforms best. Choosing the best relief visualization is clearly subjective because we are forced to visually assess the graphical results of our research. However, when comparing all outputs, we can see that visualizations without a DO (Figures 3.a, 3.c, 3.h) do not represent the detailed land cover features and elevation differences in a satisfactory manner. Setting the DO above hillshades (Figure 3.g) darkens the image and therefore cancels the shading effect. Compared to other reliefs, the relief in Figure 3.b (without the QGIS hillshade) least emphasizes the various landforms and looks pale. In our opinion the best relief visualizations are d), e) and f). After a detailed zoomed in examination of visualization in Figure 3.f which contains the multidirectional hillshade, we concluded that it creates lighter shadows and does not visualize the steep river canyon in the northeastern part of the NP as well as the other two. Also, we tested additional layer style options to enhance the relief visualization i.e., we adjusted brightness and contrast of the DO to 15% and -15%, respectively. Finally, we created our shaded reliefs for every NP shown in Figures 4-10 using following layers and parameters: coloured DEM layer (Normal, 100 %), Blender hillshade layer (Normal, 60 %), QGIS onedirectional hillshade layer (Multiply, 75 %), DO (Multiply, 75 %).

Conclusions
Although creating a visually appealing relief representation can be challenging, it is achievable using only free and open-source software. Blender is a powerful computer graphic software which offers different settings not available in conventional hillshading algorithms used in GIS software. In this research we combined layers created with different terrain rendering techniques to create more natural and softer shaded reliefs. The layer blending combination of multiple hillshades, hypsometric DEM colouring and DOs can be achieved using blending modes in QGIS. This method significantly enhances terrain visualization by improving the visual details of land cover features. Also, by applying this method to all eight Croatian NPs, we showed that it is applicable to a series of different terrain landforms and sizes. The final shaded reliefs of all NPs are applicable for multiple purposes and will be used in future to create a large format poster.