Original Research

The first sub-meter resolution digital elevation model of the Kruger National Park, South Africa

Kai Heckel, Marcel Urban, Jean-Sébastien Bouffard, Jussi Baade, Peter Boucher, Andrew Davies, Evan G. Hockridge, Wolfgang Lück, Jonas Ziemer, Izak Smit, Bernhard Jacobs, Mark Norris-Rogers, Christiane Schmullius

Received: 20 Mar. 2021; Accepted: 01 Oct. 2021; Published: 17 Dec. 2021

Copyright: © 2021. The Author(s). Licensee: AOSIS.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The use of digital elevation models has proven to be crucial in numerous studies related to savanna ecosystem research. However, the insufficient spatial resolution of the chosen input data is often considered to be a limiting factor when conducting local to regional scale ecosystem analysis. The elevation models and orthorectified imagery created in this study represent the first wall-to-wall digital elevation data sets produced for the Kruger National Park (KNP), South Africa, at very high spatial resolution. Using colour-infrared (CIR) aerial imagery from the archives of the Chief Directorate: National Geo-spatial Information (CDNGI), Department of Agriculture, Land Reform and Rural Development (DALRRD) aerial acquisition programme, we created digital surface models (DSMs), digital terrain models (DTMs) and CIR orthomosaics covering the entire KNP with a nominal ground sampling distance of 0.25 m. Elevation information was derived using state-of-the-art stereo matching algorithms that utilised semi-global matching (SGM) as a cost aggregation function throughout the image pairing, using the Enterprise software from CATALYST. The final products were validated against reference products, and showed excellent agreement with R² values of 0.99. Further, the validation of the DTM and DSM revealed median absolute vertical height error (LE90) across all sites of 1.02 m and 2.58 m, respectively. The orthomosaics were validated with in situ ground control points (GCPs) exhibiting a horizontal Circular Probable Error (CPE) of 1.37 m. The data resulting from this work will be distributed freely with the aim of fostering more scientific studies in the African science community and beyond.

Conservation implications: Accurate information about terrain and surface height are crucial inputs to a variety of scientific analysis, which are essential in protected areas, such as flood prediction or fire hazard estimation. Elevation data sets and orthomosaics in very high resolution can therefore serve as a crucial tool to improve park management and foster positive implications on conservation efforts.

Keywords: Kruger National Park; elevation; aerial imagery; semi-global matching; stereoscopy; 3D; protected areas.

Introduction

Since their introduction in 1958, digital elevation models (DEMs) have been used widely in fields such as ecology, hydrology, other Earth-related disciplines, modelling and remote sensing (Miller & Laflamme 1958). A DEM is a digital representation of a topographic surface that is defined through a two-dimensional discrete function of a morphometric variable and consists of XYZ coordinates (Florinsky 2012). Their accurate derivation is relevant not only from a scientific but also a socio-economic point of view. Originally, DEMs were created through classic field surveys. In recent decades, height structure analyses were mostly carried out using elevation information derived using remote sensing techniques such as photogrammetry based on stereoscopic imagery (Pieczonka, Bolch & Buchroithner 2011; Pulighe & Fava 2013), Light Detection and Ranging (LiDAR) (Lim et al. 2003; Lindberg & Holmgren 2017) or air-/spaceborne Interferometric Synthetic Aperture Radar (InSAR) (Hirt, Filmer & Featherstone 2010; Nelson, Reuter & Gessler 2009). On local to regional scales, aerial photos and LiDAR collected by unmanned aerial vehicles (UAVs) or airborne instruments are, because of their comparably lower degree of complexity, more commonly used than InSAR procedures from satellite-derived imagery. Since LiDAR data sets represent a direct physical measurement of height information, they are known to be more accurate than DEMs derived from aerial images (Baltsavias 1999). This arises mainly from the ability of full-waveform laser scanning systems to detect entire vertical structures of objects situated aboveground and to penetrate tree canopies (Reitberger, Krzystek & Stilla 2008). Major disadvantages of LiDAR campaigns are the costs of data acquisition, which are significantly higher than for flight missions acquiring aerial imagery (White et al. 2013), as well as the missing physical relation between plant structures and the radiation in comparison to true colour (RGB) or colour-infrared (CIR) data.

Photogrammetric methods have been used for more than a century to retrieve height information from aerial imagery. While they are now processed digitally, these methods still use the same fundamental principles of parallax measurement, taking into account different viewing angles of overlapping image pairs (Straub et al. 2013). Because of the lower flight altitude compared to satellite imagery, airborne flight acquisitions are able to reach higher spatial resolutions with adjustable overlap, and thus are less impacted by vertical as well as horizontal uncertainties. The high degree of overlap presents multiple measurements of the same point at different parallaxes with multi-angular observations, leading to more precise elevation measurements, and a more rigorous triangulation via a block bundle adjustment. Furthermore, the development of innovative techniques has increased the accuracy of elevation models that were calculated from airborne imagery in recent years (Reitberger et al. 2008; Stepper, Straub & Pretzsch 2015). These include innovations in computer vision technologies, such as the improvement of dense image matching algorithms, fostering new research using stereoscopic imagery from UAVs and airborne sensors for forest applications (Haala et al. 2010; Straub et al. 2013; Tian et al. 2017).

Global-scale data sets from DEMs have existed for decades, and many are now freely available. However, these global data sets provide only medium to coarse resolution products (≥ 30 m cell size), and fine-resolution DEMs are required to study local topographic phenomena at the regional level. High-resolution products are usually not freely available in most areas of the world. In protected areas, such as the Kruger National Park (KNP), reference data are usually scarce, yet the need for high precision and (very) high resolution elevation models is great. Several studies have indicated the need for such data products to minimise biases in vegetation modelling (O’Loughlin et al. 2016) and flood management (Wu et al. 2017).

This study aims to deliver the first publicly available very high resolution (VHR) elevation models and orthomosaics of the KNP derived from photogrammetric algorithms using freely available aerial imagery. We introduce and discuss three products: (1) a digital surface model (DSM), (2) a digital terrain model (DTM), and (3) an orthomosaic data set. A DSM represents the Earth’s surface including all aboveground features (such as vegetation or buildings), while a DTM characterises the elevation of bare ground underneath these features. The orthomosaic data sets represent a mosaic of aerial images that were corrected for geometric effects, such as perspective distortions and topographic influences. The final products were produced with pixel postings between 0.25 m and 5 m to provide the basis for a variety of ecosystem analysis applications. In addition, to allow users to utilise data products at a lower resolution, we resampled the data to 1 m and 5 m, as lower resolution products may suffice for specific types of applications (e.g. modelling vs. habitat mapping).

Data description

Study area

The KNP comprises protected land areas to the extent of approximately 19 500 km², making it South Africa’s largest protected area. It is located in the northeastern part of South Africa, bordering Mozambique in the west and Zimbabwe in the north. Elevations vary between 260 m and 839 m above sea level. The KNP has a mean annual precipitation of around 600 mm, and exhibits a gradient of rainfall that increases from north to south (Gertenbach 1980). Precipitation occurs heterogeneously over time and space. The climate alternates between contrasting seasons, with a wet and hot season (usually October to March) as well as a dry and mild season (April to September). Geologically, the park can be divided into two parts: (1) the basaltic formations (nutrient-rich) in the east, and (2) the granitoid bedrock (nutrient-poor) in the west (Venter 1990). The characteristic flora of the KNP is dominated by woody vegetation and herbaceous plants, with canopy cover ranging between 10% and 80% (Venter, Scholes & Eckhardt 2003). Vegetation is largely influenced by geology and rainfall, as well as disturbances like fire and herbivory.

Aerial imagery

To cover the entire KNP, an area of approximately 19 500 km², image data from 12 image acquisition tiles, also referred to as ‘entities’, were collected through flight campaigns accomplished at the end of the dry season in 2018. Entities comprise between 500-700 aerial images, each covering an area of roughly 20 km² adding up to an extent of approximately 60 by 60 km per tile, which can differ for some areas. The CIR raster data sets were stored in GeoTIFF file format and included three channels: near-infrared (NIR) (830 nm, red (652 nm) and green (544 nm). An overview of the aerial campaign entities is provided in Figure 1.

FIGURE 1: Leica digital mapping camera coverage of the Kruger National Park, South Africa represented as individual flight strip centroids for each acquisition tile.

The aerial images were acquired with the Leica Digital Mapping Camera (DMCIII) instrument which was utilised within the flight campaigns of the South African national mapping agency Chief Directorate: National Geo-spatial Information (CDNGI), Department of Agriculture, Land Reform and Rural Development (DALRRD) aerial acquisition programme. The DMCIII consists of eight individual monochromatic cameras, four pan cameras, and one for each spectral band. Each camera features its individual interior orientation (IO). During the calibration processing (Level 1 processing), images were projected from their IO to that of an idealised and common IO (or camera model). Pan bands were mosaicked, and multispectral data were stacked. Pan and multispectral data were then fused to create a pan-sharpened product. The Level 1 processing software uses a Brovery Image Fusion approach, allowing that only three spectral bands can be fused at a time and multispectral bands poorly correlated with the pan band get radiometrically distorted. The imagery gets captured in 12 bit and delivered in either 8 or 16 bit (16 bit imagery was used in this study). Data sets were taken from an altitude of approximately 5500 m – 6000 m above ground. The flight campaigns were carried out by GeoSpace International (Pty) Ltd. Latest data for the whole country of South Africa have been acquired at 0.5 m and 0.25 m pixel posting for the periods 2008–2016 and 2017–present, respectively. The imagery was acquired in September/October 2018 using a CIR channel combination using the NIR, red and green bands. Analysis was then carried out at a pixel posting of 0.25 m.

Validation data

To validate the products that were derived in this study, various independent data sources were used. For the elevation models, we utilised different types of height measurements as reference. In order to validate the height at ground level (DTM), we used differential Global Navigation Satellite System (GNSS)-based ground survey points that were collected during field campaigns between February 2014 and March 2016 using Leica GS10 and GS15 GNSS receivers (Baade & Schmullius 2016). The roughly 10 000 ground measurements were stratified over various land cover types following a transect sampling strategy and filtered to prevent the usage of points that were situated inside water bodies, which were dried out during the time of acquisition. For the validation of the DSM, we utilised VHR LiDAR UAV data. Data sets were collected at 100 m altitude above ground over numerous sites across the KNP, primarily focusing on experimental burn and herbivore exclosure plots, during January and February 2020 using a Riegl VUX-1 LR instrument attached to a DJI M600 Pro UAV. These sites covered a diverse set of vegetation types and structures, ranging from open to dense and from short to tall, thus providing an ideal data set for validation. The average size of each site was 110 ha with a flight overlap of 60% and a pulse repetition rate of 820 kHz, resulting in an average point density of ~150 points per m². Terrain information from the Shuttle Radar Topography (SRTM) mission was used in flight planning to homogenise LiDAR point density and beam divergence across the study area. In post-processing, UAV trajectories were refined with base station data, denoised, classified (Axelsson 2000), and matched together. Rasterised DSMs were created by calculating the maximum elevation of all 1st return points within each 0.25 m cell. After processing, the accuracy of the point clouds, and thus the DSM elevations, was expected to be ~0.10 m or less. The point clouds were pre-processed to calculate DSMs at a cell resolution of 0.25 m. The locations of all sites that were used to validate the DMC surface model are shown in Figure 2.

FIGURE 2: Distribution of reference data sets used to validate the elevation models for Kruger National Park, South Africa.

In order to validate the horizontal precision of the orthomosaics, we relied on the geolocation of triangulation stations, or ‘trig’ beacons. Shapefiles were provided by CDNGI with the exact easting and northing position, as well as the appearance and height of the trig object (e.g. pole). In total, more than 50 locations were used to evaluate the accuracy of the orthorectified aerial imagery (see Figure 2).

Methods

Image processing with CATALYST Enterprise

To extract height information from the stereoscopic DMC data, we utilised CATALYST’s (formerly known as PCI Geomatics) commercial Enterprise software environment, which is a geospatial image-processing software designed to leverage complex computing hardware and to be scalable, using the available computational resources efficiently. A major advantage is its underlying Python programming framework, enabling the development of new tools and the optimised management of resources (CATALYST 2020). Using CATALYST Professional’s archive of algorithms, a wide variety of data manipulation methods can be used for elevation extraction and orthomosaic creation. Thus, we created a workflow to implement different processing and quality assessment steps. These are automated within Enterprise, and thus could be easily applied to other study sites and data sets.

To derive height information from stereoscopic aerial imagery, we followed the workflow outlined in Figure 3. The processing steps made use of state-of-the-art routines that are applied to carry out height extraction algorithms on stereoscopic images from UAVs, airplanes, or satellites.

FIGURE 3: Workflow to extract elevation information from digital mapping camera aerial imagery using CATALYST Enterprise.

Data preparation

As visualised in Figure 3, the ingestion of aerial image data into the Enterprise production system could only be accomplished once metadata had been prepared. The metadata provide crucial information derived from the physical characteristics of the acquisition platform, in this case inertial measurement unit (IMU) and differential Global Positioning System (GPS) of an airplane. They include image centre position (easting & northing), flight altitude, as well as image orientation (Ѡ), phi (ϕ) and kappa (К), which determine the orientation of images through three different axes. These variables are essential inputs for all aerial image processing. For this study, the exterior orientations (EO) calculated by CDNGI were used. Once all available input information has been properly prepared, the data can be ingested into the Enterprise production system.

Bundle adjustments and tie point collection

The EO model of each image was refined for the purpose of improving the accuracy of the extracted DSMs and the ortho images. This was accomplished by collecting tie points and performing block bundle adjustments. The selection of tie points was carried out using a gridded sampling approach with 150 samples per overlapping area using the NIR channel as the matching variable between image pairs. This bandwidth was chosen because it is less sensitive to illumination differences, and thus more robust in terms of atmospheric effects (Zhu et al. 2019). A Fast Fourier Transform Phase matching (FFTP) approach (Brigham 1988) was utilised to automatically identify tie points between aerial images. This approach identifies points based on the relative shift between images (the phase difference in the frequency domain) as the criteria for point matching. If the correlation coefficient between two matching points was below a given threshold (0.85) within a search radius of five pixels, these points were rejected for the adjustment process, which further increased the quality of the matching process between adjacent images.

Extraction and filtering of a digital elevation model

To generate DEMs from airborne imagery, images needed to be projected from sensor to an epipolar geometry, with Y coordinates of overlapping pairs perfectly aligned and the remaining parallax represented in the X coordinates of the resulting data set, from which a measured parallax is then translated into elevation. The threshold for the minimum overlap between images was chosen to be at least 40% and the NIR channel was selected as the input for the epipolar generation. This bandwidth was chosen because the NIR band preserves contrast, and thus continues to exhibit bright reflectances over both bare surfaces and vegetation. The NIR channel is also less affected by haze and atmospheric effects, although its effective resolution is poorer than that of the shorter wavelength because of diffraction. Any remaining misalignment identified in the epipolar geometry was modelled and quantified to improve the subsequent DSM extraction procedure. Generally, height extraction from airborne data was carried out through several steps including matching cost computation, cost aggregation, disparity optimisation and refinement procedures. In recent years, semi-global matching (SGM) approaches have proven to be amongst the most popular and successful algorithms in the fields of stereo vision and photogrammetry (Klette et al. 2011; Michael et al. 2013). To extract the height information from the aerial imagery, we used this matching approach, which utilises intensity differences, mutual information (as the cost function) and an approximation of the global energy function that is being optimised path-wise (16 paths in this study) from all directions over the image. The cost function is significantly influenced by the use of penalty values, which were chosen based on performance tests and represent varying magnitudes of disparity changes. These variables have a strong impact on the matching performance and the robustness that is related to this processing step. The term ‘semi-global’ arises from the combination of both global and local methods in a way that the complexity of the process is lowered, and the quality of the matching is drastically improved. While the computation time for global methods is often considerably higher, the overall performance increases compared to local matching algorithms. Further, pixel-wise calculated matching cost, contrary to the calculation along image paths, poses negative effects of insufficient correspondences related to low texture and ambiguity (Hirschmüller 2007).

To minimise undesirable features in the elevation model, filtering is a necessary step in the extraction of height information. To minimise noise originating from correlator accuracy and to improve the accuracy of the extracted DSM, a series of spatial filters were applied. An intermediate level of filtering strength was used, this includes: noise filter, hybrid median filter, and a bi-level smoothing filter. The noise filter performs two operations: (1) a neighbourhood analysis to identify blunders by examining statistics, and (2) analysis to identify blunders by pinpointing pixels surrounded by failed pixel matches. The hybrid median filter is used to smooth noise after the interpolation and rasterisation of the correlation point cloud. This version of the median filter has been modified to preserve natural edge features and corners (e.g. trees, buildings), while still reducing noise. The bi-level smoothing is applied to further reduce noise, while preserving abrupt and steep elevation changes throughout the DSM. In addition to these spatial filters, a pattern suppression technique was applied to minimise artefacts caused by strong patterns in the landscape (e.g. rows in agricultural fields). This method examines characteristics of the accumulated costs to infer a blunder. After filtering, the DEM extraction for each tile was completed.

The impact of these filtering algorithms on the final spatial resolution of the DSM is believed to be minimal. This assumption is based on the evidence that the noise filter has no impact, and the hybrid median and bi-level smoothing may have only a minimal impact because of their small neighbourhood size and feature preserving attributes.

Digital terrain model conversion

The derivation of bare earth elevation data requires the removal of all elements that are above the ground surface (such as buildings or vegetation) of the DEM. An adaptive slope-based filtering approach was applied to remove surface features from the raster data set. To allow maximum flexibility in the elimination of aboveground elements, the filtering process was customised by defining object sizes (100 pixels) to remove surface features, as well as a small gradient threshold to preserve hills and natural slopes (slopes above five degrees of gradient were filtered). Additionally, filtering steps were introduced to remove unwanted bumps and pits (filter window of 15 pixels and a slope of five degrees) in the data that may persist from prior object removals. Thresholds for these filters and surface object extents were fixed for the entire KNP to retrieve height information with consistent parameters. The application of filters to smoothen out remaining objects that are situated above ground (above DTM height) led to a certain degree of blending of the height statistics of some of the pixels within their neighbourhood in the magnitude of the described window size. During this filtering step, it can be assumed that the height of the terrain surface was not altered as much as that of the bumps and pits, which were the main objective of filtering in this processing step.

Orthorectification of aerial photos

To create orthorectified aerial images with exact geolocations, the previously created height models were used to correct the flight strips for effects of terrain. In order to use the CIR aerial imagery for land cover and classification analysis, the aerial imagery needed to be orthorectified and prepared for image mosaicking. Because of its high spatial detail, it was essential to geolocate the data as precisely as possible. The main goal of orthorectification is the correction of view angle and relief-induced effects, as well as the georeferencing of the aerial scenes. Thus, the elevation models with a pixel size of 0.25 m served as elevation input to remove topographical effects from the imagery. The conical rays in unrectified aerial images were converted to parallel rays, which made them orthogonal to the earth’s surface (Jensen 2007).

Mosaicking

For each entity presented in Figure 1, all available flight paths were merged to create a single seamless mosaic of elevation models and orthomosaics. The so-called ‘image stitching’ of the elevation models was based on two main components, cutline computation and colour balancing, which minimised edge effects that could have occurred in the combined data set between the borders of neighbouring tiles. Firstly, cutlines were created between all adjacent images using a criterion of minimum relative difference in grey value levels in the overlapping area of neighbouring images. These boundaries define at which location adjacent images are combined to form seamless mosaics. Next, the radiometric differences within the image were reduced by applying colour balancing algorithms, which balanced histograms of the areas of overlap between the flight strips. The colour balancing process is important to remove possible patchwork effects from the mosaic. Thus, contrasts are levelled out and seams were removed, which positively impacted the image statistics.

Artefact elimination

To remove artefacts from the elevation models that remained after height extraction, filter techniques were applied to the data. Larger water bodies such as dams, lakes, or other standing water bodies led to faulty values in the elevation models, which were located in the greater negative feature space. This effect occurs because the algorithm is not able to detect tie points to connect adjacent aerial images over water surfaces, resulting in incorrect elevation values in these areas. Because of the changing appearance of the water surface and the very low contrast, elevation estimates over water are unreliable and highly variable. To account for these effects and remove them from the data sets, we used an edge filter to fill pits within water bodies using the elevation value of the surrounding bank areas. Besides the removal of water bodies, which could not be captured by the SGM algorithm, no post-processing filtering techniques were applied to the data.

Validation

Elevation models

A comparison between the highly precise in situ samples (from the GNSS ground points and UAV-LiDAR data) with the modelled elevation (DTM and DSM) from the DMC data sets was carried out using the absolute pixel-wise height difference Δh, which quantifies the deviation between the reference height information (h_r) and the respective target elevation model (h_DMC). Thus, a vertical error is measured relative to the known height. Consequently, the height differences Δh_DTM and Δh_DSM for DTM and DSM, respectively, are given as:

Since the analysis of the height deviation underlies a number of sources of uncertainty, further height measures were calculated to provide a substantial representation of the differences between data sets. Besides the pixel-specific height deviation, we calculated the median Δh_p50 of the height deviation, the standard deviation (SD), the 95th percentile, and the median of the 90th percentile of the absolute linear vertical error (LE90) of all individual reference sites. Additionally, we calculated the confidence intervals Δh_CI (95%) for Δh.

Orthomosaics

To quantify the horizontal accuracy of the geolocation of the orthorectified aerial imagery we used trigonometrical beacons as reference data. These are relatively evenly distributed across the KNP, allowing a fairly good representation of the area. In total, 89 trig stations were utilised, which vary in their appearance from beacons installed on hilltops to well-known man-made structures, such as outlooks or buildings. They are part of the physical framework for South Africa’s geodetic reference, and will be replaced with 50 GPS base stations in the next few years. However, these trig stations represented the most evenly distributed and precise source for geolocation accuracy assessments available for the KNP available for this study.

In order to quantify the horizontal accuracy of the ortho imagery, we computed several statistical measures. First, the single ortho images were merged by creating a virtual raster (VRT) data set using Geospatial Data Abstraction Library (GDAL) for the KNP. The VRT file format is a GDAL driver, allowing the computation of fully functional (statistically) raster stacks without the need to write and export the extensive data set to common raster file formats such as ‘.tif’ or ‘.pix’. This way, valuable disk space can be saved, while also speeding up the validation process (GDAL/OGR Contributors 2021). Second, the VRT’s estimated location of the trig stations was acquired. Finally, several statistics were calculated to quantify the horizontal error in X and Y direction, using the geographic coordinates. The equations for these statistics are provided in Table 1 and visualised in Figure 4.

FIGURE 4: Overview of circular error measures; error measures given: three-five sigma error, circular map accuracy standard, mean square positional error, circular probable error, mean circular standard error, directional standard error of X and Y.

Results

Products

Digital terrain model, DSM and orthomosaics were created covering the entire KNP. Cell sizes of the products are 0.25 m, 1 m and 5 m. All products are projected to UTM Zone 36, WGS84 using the WGS84 ellipsoid also as the vertical datum (if vertical information is stored). Digital elevation models and orthomosaics were partitioned according to the CDNGI 1:50 000 map sheet units for more convenient download. Each unit covers an area of maximum 640 km² or 25 km in X and Y direction (depending on the overlay with KNP borders). All data sets are freely available through the website of the Centre for Environmental Data Analysis (CEDA), and can be accessed as described in the Data Availability Statement.

Digital elevation models

As described earlier, two types of DEMs were calculated: DSMs and DTMs. Both products were filtered using the same parameters in the post-processing phase and cropped to the boundary of the KNP. Elevation spanned from 105 m to 832 m above sea level as well as 106 m to 833 m above sea level for the DTM and DSM, respectively. The median height value for the KNP DTM and DSM was 333.3 m and 333.4 m, respectively. From the surface height model, features below the horizontal extent of four pixels can be identified. This allows for a minimum mapping unit (MMU) up to a meter or less. Examples of objects which can be located in the height data set are smaller trees as well as distinctive bushes, which are typical for the savanna environment.

Figure 5 visualises both elevation models using an identical colour range as well as according to their orthomosaic subset. As visualised in the exemplary zoomed windows in the figure, the DSM allows the identification of individual trees and bushes, especially when they are foliated. Subset 2, which displays the village of Skukuza, demonstrates the precise delineation of single houses and other artificial structures. The DTM displays a very good representation of the terrain height after the removal of objects above the ground. Because of the large extent of the study area (~19 500 km²), some artefacts may not have been removed by the DTM filtering algorithms completely, and thus remain in the final products partially or to their full extent.

FIGURE 5: Orthomosaic (upper), digital surface model (middle) and digital terrain model (lower) generated from aerial imagery for four different sites across the Kruger National Park at a scale of 1:5000.

Orthorectified aerial imagery

The orthomosaics were created using CIR ortho images as well as the previously described DSM at 0.25 m cell size. After the mosaicking process, the ortho images were combined to form a seamless orthomosaic of the original Leica DMCIII tiles. Thus, the results enable the analysis of savanna ecosystems at the scale of an individual tree or shrub. An exemplary visualisation of the great level of detail of the 0.25 m orthomosaic, located along a central stretch of the Limpopo river close to the northern KNP border, displaying different scales at which the data can be used is shown in Figure 6.

FIGURE 6: Orthorectified colour-infrared aerial imagery at different scales enabling local to regional analysis; single bushes and trees can easily be delineated.

Orthorectified imagery of this scale can be used in a variety of scientific analysis, including vegetation cover analysis (Le Bris, Tassin & Chehata 2013; Lucas et al. 2000), biomass estimation (Meng, Pang & Zhang 2014; Shin et al. 2016), or the evaluation of the effects of severe flood events (Ichim & Popescu 2019).

Validation

Digital terrain model

The results of the validation of the 0.25 m DTM are shown in Figure 7. Here, the relation between the dGNSS-derived ground points and the modelled surface height from VHR aerial imagery is visualised. For the various locations across which dGNSS points were collected, the agreement with the estimated elevation was found to be very strong, with a coefficient of determination (R²) of 0.99. Local variations were found between the different acquisition areas. The median of the LE90 for all in situ points was 1.02 m with a site-based (averaged for each location) minimum of 0.12 m, a maximum of 2.06 m, and a 95th percentile of 1.79 m. Strongest deviations from the DMC-based terrain estimation were found near some of the reservoir sites. These sites were either dried out or at least only partially filled with water at the time of the dGNSS data acquisition. The best agreement with the dGNSS height estimations was found around the dried-out Hartbeesfontein reservoir, north of the Timbavati River, as well as in the sites that were located in and around Skukuza. Here, flat surfaces, such as a sports field or a runway, were measured during the field campaigns.

FIGURE 7: Validation of the modelled terrain (digital terrain model) elevation from stereoscopic aerial imagery with ground measurements for entire area and Lugmag Reservoir sub-site. More than ten study sites were included in the validation (see Figure 2 for the distribution of the sites).

Digital surface model

Similar to the validation of the DTM, the DSM was validated using a pixel posting of 0.25 m. The results are visualised in Figure 8. The validation of the DSM was carried out using VHR LiDAR data that were acquired in January and February 2020, representing the end of the wet season. Across all sample sites (with 250 validation points at each site), a trend of small underestimation of vegetation height was found in the modelled DSM. While this effect varied by location, it was present among all sites. Similar to the derived DTM, the calculated DSM exhibited a very good overall agreement with the reference data exhibiting a coefficient of determination (R²) value of 0.99. However, similar to the 95th percentile of height deviation of 5.61 m, the LE90 was found to be larger than the DTM error reference value with 2.58 m averaged (median) over all UAV acquisition locations that are visualised in Figure 2. Here, the median of the deviation varied between 0.5 m and 2.6 m across the various LiDAR validation sites. The largest deviation between the UAV-based reference data and the modelled surface height was found at the Kambeni Experimental Burn Plots (EBP), which are located near Numbi Gate, with a LE90 height difference of 6.8 m. The best fit was found at the Makhohlola exclosure plot near Crocodile Bridge, close to the Mozambique border, exhibiting a deviation between reference and modelled data of 1.6 m.

FIGURE 8: Validation of the modelled surface elevation (digital surface model) from stereoscopic aerial imagery with unmanned aerial vehicle-based lidar measurements for entire area and Shabeni Experimental Burn Plot sub-site. Seven different locations were used for the validation.

Orthomosaic

The visual inspection of the final orthorectified aerial images revealed very high agreement between existing and freely available ortho products. As described in the Orthomosaics section, the geolocation accuracy of the orthomosaic was assessed using the precise location of several ground control points (GCPs) (trig stations) that are distributed across the KNP. The trig stations were solely used for a horizontal accuracy validation, not for the correction of the aerial images. In total, 58 validation points were used to analyse the geolocation. As visualised in Figure 4, the individual circular error measures represent different levels of confidence (Circular Precision Probability [CPP]). Table 2 summarises the individual levels of error probabilities as described in the Orthomosaics section. All validation points were found to be within a 2.52 m radius of the associated pixel in the orthomosaic with a probability of 90%. The accuracies at different in situ locations varied significantly because of multiple reasons (explained in Discussion section).

Discussion

The elevation models derived in this study were found to exhibit a very strong agreement with the reference data sets described in data description section. Both models exhibited R² values of 0.99 and SD of 0.76 m and 2.03 m for the DTM and DSM, respectively. Thus, the different ecosystems and landscapes of the KNP did not significantly influence the correlation between the reference data and the estimated elevation models. Visual interpretation as well as statistical analysis of the vertical error in the derived DEMs showed a fairly good agreement between the high-resolution height information with the available reference data sets. While the in situ samples used for the validation of the elevation models were well-distributed across the entire KNP, some areas were not covered sufficiently (e.g. mountainous areas). To further improve on the characterisation of the performance of the height estimations, more reference data would be desirable in the future, especially in areas with steep slopes.

Comparing the different types of elevation models, the DTM was found to feature a lower vertical deviation from the in situ dGNSS data with a median LE90 of 1.02 m. In contrast, the derived DSM, which was validated using UAV-based LiDAR data, displayed a median LE90 of 2.58 m. Especially in areas with no or only sparse vegetation cover, the estimated DTM performed fairly well against the ground points. Likewise, the DSM effectively captured the open savanna vegetation present in the Makholola exclosure plot in the southeast of the KNP, exhibiting a LE90 of 1.6 m (open savannas on basaltic soils). The largest deviations of the DTM from the reference ground points were found around a reservoir south of the Olifants river (central KNP), with a LE90 of 2.06 m. A number of possible sources of uncertainty underlie the accuracy of the DTM. For example, one source of uncertainty could be seasonal variation occurring between the acquisition date of the aerial imagery (end of the dry season) and the collection of terrain dGNSS points (end of the wet season). Also, seasonal events such as floods could have altered ground height between aerial and dGNSS acquisitions. In open areas, such as those located around the KNP water reservoirs, the presence of vegetation (e.g. grasses) can also negatively influence the comparison between the two acquisitions, as some objects above the ground may have not been filtered out properly in the DTM conversion process. The conversion from the DSM to the DTM was carried out using several filtering techniques and smoothing parameters, including gradient thresholds and moving-window bump and pit filters. The type and aggressiveness of filtering grids was tested and selected after visual and statistical interpretation. However, large structures, such as long walls, buildings, or bridges, might not have been completely removed from the DSM, thus transforming the DTM ground height as well. Such remaining artefacts were likely to introduce differences in height values compared to the original GNSS points. In addition to affecting the vertical accuracy, the filtering mechanisms also altered the horizontal dimensions of the data sets. As a result, the pixel posting does not necessarily equal the spatial resolution, as some of the filters calculated neighbourhood-based statistics, averaging the height values of single pixels in relation to their neighbouring pixels using the thresholds and parameters described earlier. Further, the location of some dGNSS samples on the tops or edges of dams possibly led to incorrect values in the estimated DTM, because of misinterpretations in the epipolar matching algorithms.

Savanna ecosystems are vulnerable to shifts in vegetation, as well as to the loss of biodiversity because of natural or anthropogenic causes (Smit et al. 2010). Thus, savannas can be attributed as dynamic environments, where seasonal differences, often reinforced by wildfires, hugely impact the growth of mostly woody vegetation. In contrast to the pure terrain elevation, the height of vegetated areas can change drastically over time, such as after extreme droughts, exceptionally moist wet seasons or bushfires. These seasonal events combined with herbivory effects (e.g. elephant impacts) foster resprouting, regrowth, and coppicing from stems and roots (Asner & Levick 2012; Clarke et al. 2005; Mograbi et al. 2017). Consequently, seasonal differences can lead to errors in the deviation estimation between reference data and calculated height models, as data were collected 18 months apart. Thus, the growth of woody vegetation, which varies between different acquisition dates of the reference data and the aerial imagery, can lead to height differences. The seasonal differences become even more visible when analysing Δh in the estimated DSMs that were located on or around the EBP, which are scattered across the KNP. Since these sites are frequently burned (burnings every one to six years and for some plots up to decades), vegetation patterns here are highly dynamic (Chown 2010). However, the impact of these management practices also depends on environmental conditions, such as varying amounts of fire fuel or rainfall (Biggs et al. 2003).

The quality of the DSM could also have been affected as a result of the aerial imagery being acquired during the late dry season, when vegetation cover is generally sparse. Open canopies, which are common for the deciduous plants of KNP at that time of the year, can create challenges for the SGM matching algorithms. The algorithms will occasionally miss elevated pixels from open canopy vegetation during leaf-off season, instead capturing the ground underneath sparse canopies. Thus, the algorithm may fail to capture accurate height information for certain trees during leaf-off season. These circumstances led to larger height reference deviations in the DSM, especially in the structurally more complex sites such as the Kambeni burn plot, where the LE90 was found to be highest with 6.8 m. While a data acquisition campaign during the wet season would be preferable in terms of image-matching success and, consequently, vegetation structure and plant vitality mapping (e.g. through chlorophyll content – NIR interactions), wet season acquisitions are often not feasible as cloud cover heavily impedes flight missions carrying optical sensors (Hoekman & Varekamp 2001).

Besides artefacts that originated in the DSM filtering process, other effects, such as visible cutlines between entity borders or flight path edges, are still visible in some parts of the mosaic. This is especially true in areas with topographic features such as hills (or koppies) in the elevated parts of the KNP, where some image paths can cause height differences between neighbouring pixels (but these differences rarely exceed 0.5 m). Because of increased slope conditions in those areas, illumination effects, such as shadows from mountain ridges or hilltops, can severely impact the performance of the process of matching epipolar image pairs. These effects are only visible in very small portions of the elevation models and are, in most cases, negligible. However, in the case of fine-scale analyses, such as local watershed modelling, these small errors should be taken into account. Further filtering steps should also take into account the spatial variability and textural heterogeneity of the KNP.

To provide the basis for the analysis of the high-resolution aerial imagery without terrain effects, orthomosaics at the original cell size of 0.25 m were derived. The comparison to the geolocation of trig stations, as described in the Orthomosaics section, revealed a very good overall agreement between the orthorectified aerial imagery and the reference ground information. Major sources of uncertainty for the validation of the orthorectified aerial imagery arise from the way trig stations were identified in the 0.25 m DMC aerial imagery. Since the locations were estimated from the trig station metadata in the 0.25 m imagery, locating these landmarks was a challenging process that could not be automated. The difficulty of detecting reference points in the orthomosaics could be mainly attributed to the appearance of the triangulation stations, which range from windmill towers to platforms, vary in height, and can be with or without beacons. While locating tall trig stations with large poles that cast visible shadows was fairly simple, the identification of smaller stations was challenging and subject to a certain amount of subjectiveness. However, these circumstances were kept to a minimum during the selection of GCPs, which were mostly visible in the DMC imagery, which is why only 58 of the 89 available trig stations were used for the validation.

Conclusion

With the data sets created and published through this study, we present the first wall-to-wall spatial VHR elevation models and orthomosaics derived from aerial imagery for the KNP, South Africa. We produced elevation data sets (DTM and DSM) that are publicly and freely available at three different pixel resolutions (0.25 m, 1 m and 5 m). Orthomosaics (CIR) were generated at the original cell size of 0.25 m. These data sets can be used as reference information for a broad variety of savanna ecosystem analyses, including vegetation structure analysis (vertical and horizontal), biomass estimation, hydrological modelling (including erosion and flood modelling), habitat suitability modelling, or as a reference data set for other elevation models of lower resolution. To further improve data quality, we will continuously review and investigate extraction and filtering algorithms for photogrammetric elevation data in the KNP. Our goal is to derive an optimal top-of-canopy model, especially in open savanna woodlands. For future validation work, more extensive elevation data sets would be beneficial to further assess the accuracy of the models in different contexts and to improve the matching precision in the SGM algorithm. In order to make the products presented in this study even more statistically comparable to existing height derivatives, greater quantities of better stratified in situ data sets need to be utilised. This would allow for more standardised validation strategies, such as those proposed in the American Society for Photogrammetry and Remote Sensing (ASPRS) Positional Accuracy Standards for Digital Geospatial Data (ASPRS 2015). To this point, collaborative efforts in conjunction with the greater savanna science community can assist in evaluating the existing products and further refining future products.

Acknowledgements

The authors would like to thank Geospace Intl. for the personal provision of Digital Mapping Camera (DMC) data and thus enabling an early start of this project. The authors also acknowledge the support of the Chief Directorate: National Geo-spatial Information (CDNGI) in providing the validation data for the orthomosaics. Thanks also to the information technology (IT) support team of CATALYST led by Guillaume Morin, who supported us during the setup of our high-performance cluster-based processing workflow. We would also like to express our thankfulness towards CATALYST for financially supporting the presentation of results of this study at international conferences. Furthermore, this study was supported by the SPACES II funded projects Ecosystem Management Support for Climate Change in Southern Africa (EMSAfrica) and South Africa Land Degradation Monitor (SALDi).

The authors declare that they have no financial or personal relationships that may have inappropriately influenced them in writing this article.

Processing and evaluation of the presented products as well as the preparation of the manuscript were carried out by K.H. and M.U.. Major technical support to derive the final datasets provided by J.S.B., J.B., P.B., A.D., E.H. and J.Z. supported the validation strategy that was accomplished. W.L. and I.S. provided necessary methodological input for the study. B.J. provided the initial aerial imagery. M.N.-R. initiated the data delivery process. W.L. and C.S. contributed to the conceptualisation of the study. All co-authors contributed towards and/or commented on the manuscript.

This article followed all ethical standards for research without direct contact with human or animal subjects.

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

The data sets derived in this study are freely accessible through the Centre for Environmental Data Analysis (CEDA), and are available for download via the following link: http://dx.doi.org/10.5285/deab4235f1ef4cd79b73d0cbf2655bd7. Updates of these products will be continuously added and accordingly announced on the platform. Users of the data products should reference this article as source of the data.

The views and opinions expressed in this article are those of the authors and do not necessarily reflect the official policy or position of any affiliated agency of the author.

References