Terrain representation using orientation

ABSTRACT A terrain data model using orientation rather than elevation permits more efficient analysis and stores its data in a multi-band raster. Representation techniques from the computer graphics industry are readily adopted with this data model. A common data model for terrain surfaces–the raster digital elevation model (DEM)–carries with it limitations by emphasizing height. Derived products such as relief shading require additional processing to determine orientation, even though they are used more frequently than those relying on elevation (e.g. hypsometric tinting). We show some of the benefits of encoding and analyzing terrain based on surface orientation, an approach that uses normal vectors stored as multi-band raster, the data storage convention in the computer graphics industry. A change in the data model and the conceptualization of the surface’s defining characteristics allows relief shading methods to run faster than conventional tools. Processing efficiencies are especially useful for more advanced shading models that may employ hundreds of relief shading calculations. In addition, an orientation-focused approach to terrain permits cartographic techniques to parallel common computer graphics methods. This project explores one such method, normal-mapping, an effect that adds texture to conventional relief shading by perturbing surface normal vectors.


Introduction
Computer graphics software is able to render highquality, detailed displays with modest computing resources. The 3D computer graphics software packages are no longer the exclusive province of specialized animation studios. Cartographers are beginning to effectively use open-source software such as Blender to augment their terrain representation workflows. GIScience does not yet have similar capabilities for modeling terrain. One reason for this is that the default data model used for terrains is different than that used for 3D rendering in the computer graphics industry. The digital elevation model (DEM) is the de-facto standard model used to represent terrain. This height field dataset is the common input data for virtually all terrain analysis workflows. The model has an impact on how we think of the terrain, which in turn affects analysis and display. For example, shaded relief analyzes and displays surface orientation based on a complex approximation of the surface normal vector, which is derived from the height field internally.
This research explores surface normal vectors as a means of modeling a terrain surface and shows a pair of case studies illustrating the advantages. Terrain rendering workflows can be made to run faster.
Additionally, terrain orientation allows for mature applications to be adapted from the computer graphics industry. We demonstrate how textures can be properly embossed onto a shaded relief map using a bumpmapping technique using surface normal vectors.

Background
Terrain utilizes a number of data models, but commonly, it is modeled and distributed using a gridded raster data format-the DEM (Maune & Nayegandhi, 2018). The raster DEM holds the average height or elevation associated with a cell. Such raster models are well-suited and widely used for the representation of continuous phenomena such as elevation (Bolstad, 2016). Most terrain analyses in GIS operate on this data model and it is the de-facto standard format used for many surface analyses (Wilson & Gallant, 2000).
Surface analysis in analytic cartography (Bolstad, 2016;Moellering, 2000) can be categorized by whether the surface position or the surface orientation is being analyzed. Position refers to where in three-dimensional space a given point can be found, including its elevation or height. Orientation refers to the direction in which the surface is facing at that point, regardless of height. Hodgson (1998) argues that orientation metrics, specifically slope and aspect, are more commonly used than elevation in terrain analysis. Esri indicated that in 2020, surface tools that rely on orientation (Slope, Aspect, Curvature, and Hillshade) are 2.5 times more frequently used than tools that rely on elevation or height (Contour, Viewshed, Cutfill) based on geoprocessing tool usage statistics provided through their Esri User Experience Improvement Program.
Both types of analysis currently use the height field of the DEM as input. Elevation data directly from the DEM can be readily used for operations such as defining contours. Orientation analyses such as slope and aspect derive their outputs by first interpreting orientation via a focal raster function operating on the DEM (Clarke & Lee, 2007;Jones, 1998;Skidmore, 1989). Shaded relief also relies on orientation, which captures variations in both slope and aspect.
Shaded relief is an important geovisualization beyond slope and aspect as it enhances reader comprehension of terrain shape, giving a realistic impression of how light falls on the modeled surface. Shaded relief's effectiveness has been the subject of many usability studies (Biland & Çöltekin, 2017;Farmakis-Serebryakova & Hurni, 2020;Raposo & Brewer, 2013), confirming its importance in conveying shape to the reader. It is an important complement to other symbols such as contour lines (Collier et al., 2003), and of benefit to nonspecialist readers who may lack experience with abstract symbology (Patterson, 2002). Additional processing to add atmospheric realism , further augments the effectiveness of shaded relief.
Shaded relief created by hand is an expensive and time-consuming process (Marston & Jenny, 2015). Analytic methods have made the automated generation of shaded relief layers much more available, and fast to produce. Such techniques for relief shading rely on mapping surface reflectance (Horn, 1981). Most typically, reflectance maps are based on Lambert's cosine emission law. The canonical Lambertian shading method considers a single light source to compute brightness values on the surface. This point light source is set infinitely far away, rendering all light rays effectively parallel (Ding & Densham, 1994).
More complex shaded relief mapping relies on multiple light sources, which also increases the level of computational effort. The Multi-Directional Oblique Weighted method (Mark, 1992) and other techniques derived from it (Loisios et al., 2007) compute four conventional shaded layers, then sums them according to a weighting scheme. The uniform sky model (P. J. Kennelly & Stewart, 2006) and subsequent general sky model (P. J. Kennelly & Stewart, 2014) aggregates many more single-light shaded relief layers, potentially many hundreds of layers, in a weighted sum as a means of approximating real-world illumination. Sky models allow for any number of lights, arranging them in space according to various distribution models (clear sky, direct light, turbid, overcast, etc.).
Such complex and repetitive computations can benefit from an orientation model using surface normal vectors-the unit-length vector emanating from the locally tangent plane at a point on the surface. In developing a terrain orientation model, we chose to use vectors and their notation derived from surface gradient (Horn & Sjoberg, 1979) to align with the computer graphics industry. This differs from the typical approach in GIS of deriving surface normal vectors from slope and aspect. The resulting surface normal vectors can then be utilized for computations with Lambert's cosine emission law (Wiechel, 1878). Lambert predicts surface reflectance to be proportional to the angle formed between a illumination vector and the surface normal vector (Lambert & DiLaura, 2001). Illumination direction is specified as a unit-length vector as opposed to azimuth and inclination above the horizon.
The industry standard of 3D computer graphics (animations, video games, etc.), is to operate directly on surface normal and illumination vectors (F. F. Dunn & Parberry, 2011;Foley, 1995). In this model, the shader-that part of the workflow responsible for computing shade values-is a lightweight piece of the rendering calculations, with a limited role. Relief shading, by contrast, operates on the height field as its primary input and only arrives at orientation obliquely using calculations of slope and aspect. Orientation must be inferred internally (Horn, 1981;Ritter, 1987;Sobel, 2014), encumbering the shading algorithm with orientation calculations.
While the Lambert cosine emission law relies on computing the cosine of the angle theta described in Figure 1, it is not necessary to know the angle itself. The cosine of the angle between two vectors can be more readily computed via the dot product (Horn & Sjoberg, 1979;Zhou, 1992) of those vectors.
where n is the surface normal vector and i is the illumination vector. Equation 1 is a simplified form of the cosine emission law; it depends on both of these vectors having previously been normalized to a length of one unit. Relief shading calculations of a Lambertian shader (Esri, n.d.b;Grass Development Team, 2021) are algebraically equivalent to this vector-centered approach but rely on trigonometric functions after slope and aspect are computed separately. Surface orientation as a distinct data layer from DEMs allows for a more modular approach to relief shading. For example, a different model of surface reflectance can be readily implemented by changing the shader element, leaving other rendering calculations unchanged. The surface normal data can also be manipulated separately from shading to achieve other effects, such as applying a texture to the surface Texture can serve as a design element in terrain representation (Bertin, 2010;P. J. Kennelly, 2015). Texture can be added either before or after relief shading. Nighbert Nighbert's technique (Nighbert, 2002) can be broadly classified as a type of bump-mapping. By contrast, (Patterson, 2002; Figure 2(b)) embossed rock outcrops and trees after relief shading, using an image editing tool.
Bump-mapping describes computer graphics methods that manipulate either the surface geometry or its surface normal vectors to apply a texture. Offsetting the surface, as implemented by Nighbert, is more specifically referred to as displacement mapping. A much more commonly used method in computer graphics is normal-mapping (Blinn, 1978;Mikkelsen, 2008) which Figure 2. Two methods to add texture to shaded relief. In a) texture is created by offsetting z-values of the DEM [Image from Nighbert (2003)]. In b) texture is added after relief shading [Image from Patterson (2002)]. adjusts surface normal vectors prior to shading. A terrain orientation model based on surface normal vectors is well aligned with normal-mapping, as detailed in our methodology.

Surface normal vectors
Surface orientation is derived by first calculating the rate of change in z along the x and y directions (Horn & Sjoberg, 1979;Zhou, 1992). The "rate-of-change" phrasing describes the partial derivatives @z @x and @z @y associated with gradient along the axes. These gradients may be used to compose the surface normal vector without requiring slope or aspect to be calculated first (Horn & Sjoberg, 1979;Ritter, 1987;Zhou, 1992).
Gradient may be estimated in a number of ways. (See, Jones (1998), Clarke and Lee (2007), and Skidmore (1989) for comparisons of several approaches). Gradients are often calculated via Equation (2) (Horn, 1981;Sobel, 2014) using the kernel shown in Figure 3.
The cell grid width and height are Δx and Δy, respectively. Esri's ArcMap and ArcGIS Pro use this equation in their analysis tools (M. M. Dunn & Hickey, 1998;Esri, n.d.b;Hodgson, 1998) as does the shaded relief tool within GRASS (Grass Development Team, 2021). These partial derivative values lead to the surface normal vector using a method described in Horn and Sjoberg (1979) and Zhou (1992), summarized below in Equation (3). Changes in z along the x and y axes (computed using Equation (2)) are used to compose the basis vectors of the locally tangent plane (p and q). The surface normal vector extends orthogonally from this plane, so may be obtained via cross product of the p and q vectors.
It can be shown that this is equivalent to Which is then normalized to unit length: The conventional method within the computer graphics industry for storing the three components of the surface normal vector is to encode them separately in three bands of an RGB color image (Preppernau, 2020). This convention stores the x, y, and z vector components in the following bands: N x in red, N y in green, and N z in blue. This multi-band storage standard has the additional advantage that bands of data within the same raster share spatial reference, projection, cell size, and alignment (Esri n.d.a). The three bands representing the components of the surface normal vector for an octagonal pyramid rising from a plane are included in Figure 4. The pyramid sides of this sample terrain slope at 45 degrees from horizontal. The planimetric nature of the terrain surface results in similar histograms for N x and N y , but a different histogram pattern for N z . The x component, N x (in red) shows histogram peaks at 0 for flat areas and those areas facing north and south; � 0:7 for those 45 slope sections facing east and west; and � 0:5 for those panels that face the intercardinal directions. N x and N y have similar histograms as each can occupy the full range of possible values between -1 and 1. For terrains, N z can only be positive (pointing upwards) and approaches zero as the terrain becomes more vertical. Flat areas have a z component of exactly 1; sloped areas in this scene are +0.7, regardless of aspect.
When these bands are rendered with a standard image viewer, the surface normal vectors display as a RGB color image ( Figure 5), blending the bands. This color representation of the surface normal vector field is a useful visualization of the surface orientation common in computer graphics (Preppernau, 2020). This symbology, however, may not be intuitive to map users, as cartographic displays typically use grayscale relief shading to illustrate changes in surface orientation Figure 5(b). The RGB display of surface normal vectors does not require illumination information to reveal orientation. The blue cast to Figure 5(a) is the result of a large positive z value throughout the scene. Facets to the SE and E appear magenta due to mixing of red and blue. Facets to the NE and N appear cyan due to mixing of green with blue. The same scene rendered as shaded relief using light from the northwest at 45 degrees above the horizon Figure 5(b) shows a three-dimensional effect in shaded relief but does not uniquely describe surface orientation. For example, the NE and SW facets have the same brightness value.

Manipulating surface normal vectors
The multi-band raster containing the three components of the surface normal vector can now be easily manipulated using conventional computer graphics approaches. The aforementioned normal-map involves the steps shown in  profile view in Figure 6. Figure 6(a) shows the desired effect of triangular bumps. Figure 6(b) shows the surface normal vectors that would be associated with Figure 6(a). Finally, Figure 6(c) shows the normal-map, a planar surface. If this normal-map were applied to a flat terrain surface, the surface normal vectors for relief shading would be as shown.
For each point on sloping terrain surfaces, the planar normal-map is held tangent, and the normal-map's reoriented surface normal vector replaces the surface normal vector of the terrain. (F. F. Dunn & Parberry, 2011). The resulting shading gives the illusion that the textural pattern associated with the normal-map is embossed on the terrain surface.
The example in Figure 7 illustrates the orientationspecific effect of perturbing a terrain surface's normal vectors using a normal-map. The underlying surface (Figure 7(a)) is Wizard Island, in Crater Lake, Oregon, USA. The sample normal-map (Figure 7(b)) is a "diamond plate" pattern with regular features. Applying this texture to Wizard Island shows the resulting effect (Figure 7(c)). The bumps are shaded as if oriented to the underlying surface rather than if a consistent displacement were added to all of the terrain.

Data sources & techniques
Our study area is a larger extent of Crater Lake, Oregon, USA. Source elevation data was obtained as a DEM raster from the collection of reference models described in P. J. Kennelly et al. (2020) and downloaded from https://shadedrelief.com/SampleElevationModels/. Land cover from the National Land Cover Dataset (NLCD; Homer & Fry, 2012; obtained from https://www.mrlc. gov/data) is employed to define areas for varying textures with normal-maps.
Methods presented here have been implemented in a custom toolbox coded in Python for use within ArcGIS Pro. This toolbox is available as a download from the Zenodo repository at Trantham and P. Kennelly (2021) "Orientation Toolbox for ArcGIS Pro" https://doi.org/ 10.5281/zenodo.5801684. All of the tools related to methods discussed in this work are summarized in Appendix A (see Supplementary Materials).
We begin by converting the DEM into a multi-band raster dataset that uniquely identifies the orientation of each grid cell. The associated surface normal vectors are saved as their own dataset, which is consumed by all subsequent processes in lieu of a DEM. These processes include a relief shading from a single illumination  source, a more efficient sky model methodology, and a normal-mapping technique that applies a properly illuminated normal-map texture to relief shading.
Our relief shading tool using one illumination direction mimics the output of conventional relief shading but uses the dot product of vectors instead of trigonometric functions based on slope and aspect derived from the DEM. The output from our Lambertian shader produces brightness values as floating-point numbers in the range ½À 1::1� to agree with the convention in computer graphics, meaning the values are unclamped. The long-standing convention within GIScience is to clamp brightness values to the range ½0::255� as an unsigned integer (Esri, n.d.b;Horn, 1981). Clamping can discard some values and may also re-scale the remaining output to fill the defined range, which may result in loss of detail as described in Preppernau (2020).
Typical shaded relief brightness values below zero are "self-shaded" areas and are displayed as black by resetting negative numbers to zero. This occurs when the angle between the illumination and surface normal vectors is greater than 90 degrees (Ding & Densham, 1994). The remaining values are then scaled to the range ½0::255� and cast as integers. This mechanism is useful to facilitate the storage of BVs as 8-bit unsigned integers in the output raster. This is the default output range and behavior for relief shading calculated by Esri, as well as open-source tools GRASS (Grass Development Team, 2021) and GDAL (GDAL, n.d.) (upon which QGIS relies for much of its geoprocessing library). Our method honors all values between negative one and one, and uses soft shading techniques as described by Preppernau (2020) to preserve detail.
Our sky model tool allows us to measure performance differences between our methodology and traditional relief shading by executing relief shading hundreds of times (P. J. Kennelly & Stewart, 2014). This method mimics the appearance of real-world lighting, in that ambient light is modeled by the inclusion of many lights of low intensity (Figure 8(b)). Our toolset is packaged with several predefined sky model configurations, but still permits the user to supply their own custom lighting configuration via a text file. The pre-defined configuration we use is lit from the northnorthwest (azimuth 337.5) as discussed in Biland and Çöltekin (2017).
In addition to more efficient relief shading, an orientation approach to terrain surfaces allows us to parallel computer graphics methods such as normalmapping. The land cover dataset (Figure 9) is used as a mask to delineate the areas receiving various textures, much as Nighbert (2003) described. We apply a normal-map wave pattern to cells in the"Open Water" category and a forest canopy pattern to cells in the "Evergreen Forest" category. These textures provide an additional visual variable that overprints relief shading in a consistently illuminated manner and complements hue to differentiate the land covers.
The two normal-map tiles used in our example are shown as shaded relief in Figure 10. These were obtained from on-line libraries at https://polyhaven.com/textures and https://ambientcg.com/. Creating custom textures is also possible using standard image manipulation software such as Adobe Photoshop or GIMP. These normalmap images carry no spatial reference and are not intended for geospatial applications. In our method, normal-map images inherit pixel size and projection information from the underlying terrain.
The forest normal-map we use includes isolated areas where the effective slope of the bump is as much as 19 degrees. Assuming a light elevation of 45 degrees above the horizon, applying such a normal-map to a terrain sloping as little as 26 degrees away from the illumination source will result in an incident angle greater than 90 degrees and cause self-shading. This underscores the importance of using unclamped brightness values to provide the highest level of detail possible in relief shading with normal-mapping.

Relief shading comparisons
Before exploring efficiencies of our approach, we compare results based on BV of our relief shading tool using one illumination direction with the traditional tool packaged with ArcGIS Pro. To do so, we clamp brightness values to the same range of [0.255]. To measure differences in outputs, a raster comparison subtracts the conventional shaded relief output from our output. The outputs, while not identical, are comparable. Of the 27,040,000 raster cells in the 5200 × 5200 Crater Lake study area, 99.773% show no difference between the two methods' output. The remaining 61,396 cells show a difference of � 1 BV, with a similar number of grid cells in each category.
Multi-directional relief shading such as sky models involve hundreds of shaded relief layers. Performance gains per layer will aggregate to achieve significant savings on the method as a whole. The implementation of the sky model method we compare against is provided as a component of the TerrainTools suite (Field & Beale, 2016). TerrainTools is a toolbox meant to be accessed from within Esri desktop software. The sky model tool is also available from P. J. Kennelly and Stewart (2014).  Performance improvements are observed with our sky model method when compared to Terrain Tools. These improvements aggregate and are proportional to the number of relief shaded layers. Figure 11 summarizes each tool's performance as a function of number of illumination directions. Our test configurations vary the number of light sources (10-250) while holding constant other configuration parameters (location of dominant light source, sky illumination model). Data for runs with and without shadows is included to illustrate performance differences only. Shadows are included in the output of each tool.
The performance differences between TerrainTools and our sky model method cannot be ascribed exclusively to the change in the data model. The TerrainTools implementation requires that the source DEM be read from disk multiple times, for example. The elimination of this I/O overhead is a substantial factor in the difference in execution times seen in Figure 11. By reducing run time, our sky model method is more accessible to cartographers with modest computing resources. It also improves usability by simplifying the configuration interface and providing pre-defined illumination models (e.g. clear sky, overcast sky, etc.).

Normal-mapping
Textures and colors are applied to open water and evergreen forest land cover types as defined by NLCD data (Figure 9). The normal-maps are designed as seamless repeatable tiles to permit their application over a scene of arbitrary size. Figure 12 compares a conventional relief shading with one in which the normal-mapped pattern of trees and water waves are applied. The texture from normal-mapping is properly shaded with respect to the orientation of the terrain and helps to disguise the  otherwise apparent pixelation in the color of the forested areas. Also, drainages in the south-central portion of the map display are still visible due to the proper normal-map texture shading of the evergreen forest.

Conclusion
The common conceptualization of terrain as a raster DEM data model is limiting for symbolizing and analyzing terrain. The DEM's focus on z-value requires that additional processing be undertaken for common methods such as relief shading. An orientation data model adapted from the computer graphics industry has advantages for terrain representation and surface analysis. Our surface normal vector model offers performance benefits as well as access to techniques from the computer graphics industry such as normalmapping.
A multi-band raster data set holding surface normal vectors is a pre-processed data structure ideal for inferring orientation. This adds consistency and efficiency to the analysis because the surface normal calculation is made once and the result stored. The persistent dataset of surface normal vectors becomes the authoritative source for orientation in all subsequent analyses.
Relief shading based explicitly on an orientation dataset benefit by a single-direction shading operation executing faster. This performance gain accrues with the number of light sources in the method. Multidirectional methods such as sky models realize significant performance gains.
New opportunities for representing terrain arise when the surface data model changes from elevationto orientation-based. Normal-mapping, as demonstrated in this research, manipulates surface normal vectors to apply a rendered texture to relief shading. A texture normal-mapped into the terrain effectively becomes part of the terrain; it responds to lighting conditions just as the terrain does. The normalmapped texture is shaded using the same lighting as that which shades the terrain. This approach complements other techniques such as color overlays by imprinting textures into the relief shading as a visual variable.
A modular workflow in which the shader is separated from other calculations can add flexibility to shaded relief. It may be an advantage or a preference, depending on the terrain, to apply a non-Lambertian model of surface reflectance. A Lommel-Seeliger shader, for example, is well suited to modeling a "dusty" surface effect, and has a history of being used in cartography (Batson et al., 1975). Other shading models, such as Oren-Nayar (Oren & Nayar, 1995) adds a parameter to control roughness, which lends a matte finish to the surface. Both of these shaders and many others are readily adapted to the data model we propose. Specular reflectors for terrain can also be derived from surface normal vectors, illumination vectors, and an additional vector indicating the observer's viewpoint.
Surface normal vector data may also be manipulated at any stage between calculating the multi-band raster and shading. Normal-mapping is one such method, as illustrated above. As another example, terrain orientation can be generalized, for example, using low-pass filters on the components of the normal vector rather than on the DEM (P. Kennelly, 2021). Any kernel filter can be applied all three components, or to individual bands in the multi-band raster.
Surface normal vectors have broad application outside of the context of a raster representation of the terrain. Triangular Irregular Networks (TINs) are more similar than DEMs to the polygon meshes used in computer graphics, so techniques translate more directly. Each polygon within the TIN may be normalmapped to provide the otherwise flat surface with properly shaded texture. Other techniques such as Phong shading (Phong, 1975) and related techniques (Gouraud, 1971) make use of surface normal data to remove sharp transitions between facets of a polygon mesh.
With the enabling technology of surface normal vectors as a standard, the techniques that have matured and been optimized within a computer graphics context can find ready adoption in cartography. The capability to render objects very quickly, with fine control over the interplay of light, surface, and observer orientations may in turn open other possibilities for cartographic displays. In addition to the surface manipulations used to achieve specialized visual effects, the speed-up in shading would allow more interactive relief shading and enable exploratory interaction with the terrain. This interaction would lend itself to providing the cartographer with effective previews of configuration choices during relief shading.

Data Availability Statement
Python source code for the ArcGIS Pro toolbox developed as part of this project, along with associated data and configuration files, is available in the Zenodo repository at Trantham and P. Kennelly (2021) "Orientation Toolbox for ArcGIS Pro" https://doi.org/10.5281/zenodo.5801684

Disclosure statement
No potential conflict of interest was reported by the author(s).