Short Communication

Woody cover mapping in the savanna ecosystem of the Kruger National Park using Sentinel-1 C-Band time series data

Marcel Urban, Kai Heckel, Christian Berger, Patrick Schratz, Izak P.J. Smit, Tercia Strydom, Jussi Baade, Christiane Schmullius

Received: 13 Jan. 2020; Accepted: 02 June 2020; Published: 04 Aug. 2020

Copyright: © 2020. 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 savanna ecosystems in South Africa, which are predominantly characterised by woody vegetation (e.g. shrubs and trees) and grasslands with annual phenological cycles, are shaped by ecosystem processes such as droughts, fires and herbivory interacting with management actions. Therefore, monitoring of the intra- and inter-annual vegetation structure dynamics is one of the essential components for the management of complex savanna ecosystems such as the Kruger National Park (KNP). To map the woody cover in the KNP, data from European Space Agency’s (ESA) Copernicus Sentinel-1 radar satellite (C-Band vertical–vertical [VV]/vertical–horizontal [VH]) for the years 2016 and 2017, at 10 m spatial resolution and repeated acquisitions every 12 days, were utilised. A high-resolution light detection and ranging (LiDAR) data set was reclassified to produce woody cover percentages and consequently used for calibration and validation. Woody cover estimation for different spatial resolutions was carried out by fitting a random forest (RF) model. Model accuracy was assessed via spatial cross-validation and revealed an overall root mean squared error (RMSE) of 22.8% for the product with a spatial resolution of 10 m and improved with spatial averaging to 15.8% for 30 m, 14.8% for 50 m and 13.4% for 100 m. In addition, the product was validated against a second LiDAR data set, confirming the results of the spatial cross-validation of the model. The methodology of this study is designed for savanna vegetation structure mapping based on height estimates by using open-source software and open-access data, to allow for a continuation of woody cover classification and change monitoring in these types of ecosystems.

Conservation implications: Information about the state and changes in woody cover are important for park management and conservation efforts. Both increasing (e.g. because of atmospheric carbon fertilisation) and decreasing (e.g. because of elephant impact) woody cover patterns will have cascading effects on other ecosystem processes such as fire and herbivory.

Keywords: woody cover; earth observation; LiDAR; radar; machine learning.

Introduction

Savanna ecosystems are dominated by different densities of grasses and woody plants with inter-annual changes because of dry and wet seasons. They cover half of the African continent and around 20% of the global land surface and are of great significance for ecology (e.g. living environment) and economy (e.g. fuelwood, timber) (Kutsch et al. 2008; Main et al. 2016; Scholes & Walker 1993). The savanna ecosystems of South Africa are shaped by disturbance processes such as droughts, fire and herbivory (Druce et al. 2008; Scholes & Archer 1997; Stevens et al. 2016), as well as anthropogenic impacts like climate change (e.g. increase in atmospheric CO₂) or management actions. It is therefore not surprising that these savannas have been undergoing various changes during the last decades (Buitenwerf et al. 2012; Skowno et al. 2016). Information about changes in woody cover and above-ground biomass (AGB) in national parks (e.g. Kruger National Park [KNP]) is important for park management and conservation efforts, as changes in woody cover are likely to have effects on other ecosystem patterns and processes. For example, an increase in woody cover will lead to a reduction in grass and herbaceous biomass (Berger et al. 2019), which will have cascading effects on herbivores (i.e. favouring browsers to grazers [Smit & Prins 2015]) and fire regimes (i.e. reducing fire frequency [Smit et al. 2012]). Furthermore, woody cover maps are useful for a range of other park management applications. For example, the KNP has used earlier versions of woody cover maps for adjusting rhino census estimates based on visibility (Ferreira et al. 2015) and understanding how woody encroachment may impact the large mammal communities of the park (Smit & Prins 2015). Also, although remote sensing monitoring of woody cover will not be able to replace extensive field monitoring of woody vegetation composition (Kiker et al. 2014), it does provide a cost-effective and wall-to-wall approach for monitoring fractional cover, which is considered as a satellite remote sensing essential biodiversity variable (Pettorelli et al. 2016).

The recent launch of new satellite missions, such as the Copernicus Sentinel fleet of the European Space Agency (ESA), has led to a tremendous increase in the availability of freely available earth observation (EO) data that can be utilised for woody cover-monitoring applications. Sentinel-1A/B, two C-Band Synthetic Aperture Radar (SAR) satellites, have a revisit time of up to 12 days in South Africa and acquire images with a spatial resolution of 10 m (Attema et al. 2009).

Estimation of woody cover by using EO data from different sources across different wavelengths was investigated in various studies (Bucini et al. 2010; Higginbottom et al. 2018; Main et al. 2016; Urbazaev et al. 2015). However, the reproducibility of these applications might be limited because of the utilisation of EO data, which are not systematically acquired or cost intensive (e.g. Japan Aerospace Exploration Agency - Advanced Land Observing Satellite [JAXA’s ALOS], aerial images and so on). The goal of this article was to introduce a workflow to derive woody cover information from freely available Sentinel-1A/B time series and light detection and ranging (LiDAR) data (Smit et al. 2016). This workflow was applied to produce a wall-to-wall woody cover map in different spatial resolutions for the KNP for 2016 and 2017. The estimation of woody cover by using C-Band data has shown great potential in the South African Lowveld region (where the KNP is located) because of deeper penetration of the radar signals into the open savanna vegetation (Mathieu et al. 2013). The workflow is publicly available to allow the continuation of woody cover monitoring, transferability to other savanna regions and future woody cover change mapping.

Data and methods

Sentinel-1 data

The wall-to-wall woody cover estimation for the KNP was carried out by using dense time series information from Sentinel-1A/B C-Band data between January 2016 and April 2017. The Sentinel-1 data were retrieved via the Copernicus Open Access Hub (https://scihub.copernicus.eu/ – Copernicus Sentinel data [2016, 2017]). The data sets were acquired in Ground Range Detected (GRD) format from relative orbit 145, which covers the entire KNP within one overpass. The temporal resolution of the time series is 12 days and its spatial resolution is 10 m. Sentinel-1 SAR backscatter in VV (vertical–vertical/co-polarisation) and VH (vertical–horizontal/cross-polarisation) were used as inputs for the woody cover prediction model.

The Sentinel-1 data were pre-processed by using GAMMA routines (Wegmüller et al. 2016) implemented in pyroSAR, a free Python framework for large-scale SAR satellite data processing (Truckenbrodt et al. 2019), which provides a user-friendly solution for the pre-processing of SAR satellite data from recent and historical EO missions. PyroSAR also offers the possibility to use open-source software packages (e.g. ESA’s Sentinel Application Platform [SNAP]), resulting in comparable backscatter information in comparison with GAMMA.

The pre-processing steps for the Sentinel-1 data include (1) radiometric calibration to convert from digital values to radar backscatter; (2) orthorectification, utilising height information from the freely available Shuttle Radar Topographic Mission (SRTM) (United States Geological Survey [USGS] 2017) with a spatial resolution of 30 m and the precise orbit state vectors (precise orbit ephemerides [POE]); and (3) terrain flattening utilising SRTM to correct for radiometric differences caused by local incidence angles (Small 2011).

Light detection and ranging data

A LiDAR data set with a spatial resolution of 2 m representing vegetation heights was reclassified to woody cover estimates. These were used as input for the random forest (RF) modelling. The data set is available for downloading through the Carnegie Airborne Observatory (CAO) maps server (Smit et al. 2016). The LiDAR data were acquired by the CAO (Asner et al. 2007) at the end of wet season in April/May 2010, 2012 and 2014, respectively, and covered a rectangular footprint of approximately 52 km² (2 km × 26 km) (Smit et al. 2016). For this study, we used the data from 2014, to reduce the time gap between the training data and Sentinel-1 time series and thus diminish the influence of land cover changes. It needs to be mentioned that abrupt disturbances (e.g. fire, herbivory and drought) during that time are influencing the woody cover estimation. An independent validation was carried out using a LiDAR data set collected in May 2012, which covers the areas (500 m of each riverfront) along the Sabie, Olifants and Letaba River in the KNP (Milan, Heritage & Tooth 2018). The conversion from LiDAR vegetation heights to woody cover percentages was carried out for both data sets identically, which is dealt with in the next section.

Methodology

Preparation of training and prediction data

As a first step, both LiDAR Canopy Height Model (CHM) was converted to percentage woody cover. We assumed that each pixel with a height value below 0.5 m represents ground or grassland pixels, whereas all other pixels represent woody cover pixels (Smit et al. 2016; Urban et al. 2018). The resulting binary mask (0 = non-woody, 1 = woody) was then used to calculate the woody cover percentage for a regular grid with a cell size congruent to the spatial resolution of the Sentinel-1 (10 m) (sum of the woody pixel divided by the amount of LiDAR pixel fitting in one Sentinel-1 cell and multiplied by 100). The individual Sentinel-1 images were co-registered, mosaicked for the area of the KNP and combined into a layer stack for each polarisation (VV and VH), respectively. Both Sentinel-1 stacks (predictor variables) were then fused with the LiDAR woody cover (response variable) and utilised as input data for the training of the RF model (Figure 1).

FIGURE 1: The methodological framework for deriving woody cover information from Sentinel-1 time series by using light detection and ranging information (Scripts 1–4 are freely available via zenodo.org).

Random forest model

The woody cover estimation was performed utilising the decision tree classifier RF (Breiman 2001) through an implementation of the statistical software R (Machine Learning in R [MLR] package [Bischl et al. 2016]). The modelling workflow included (1) tuning, (2) spatial cross-validation, (3) training and (4) prediction (Figure 1). During model fitting, the independent decision trees are controlled by two main parameters: mtry, which describes the number of available predictor variables to split each node, and ntrees, the number of regression trees (Breiman 2001). These so-called ‘hyperparameters’ were optimised (tuning) before model fitting, as they were found to be the most sensitive variables with respect to accuracy and computation time, to identify the best parameter set for the spatial cross-validation and training (Heckel et al. 2020; Huang & Boutros 2016).

The tuning was accomplished for mtry [1, 2] and ntrees [10, 50, 100, 300, 700] with five repetitions and internal spatial cross-validation for each parameter combination. With an increasing mtry value, the runtime of the model rises likewise because of the large number of observations. Therefore, we conducted a grid search approach by using a priori knowledge of the hyperparameter search space utilising only a limited number of mtry values. The tuning revealed that the optimal parameter set, resulting in the highest accuracies and best model computation duration ratio, is mtry = 1 and ntrees = 300.

According to this, spatial cross-validation (Brenning 2012) and training of the model were performed utilising these parameters. We used a spatial cross-validation with five folds and 100 repetitions to ensure the derivation of stable results, for example, consistent accuracies. The trained model was utilised for the prediction on the Sentinel-1 time series for the entire KNP to produce the first wall-to-wall woody cover map for the year 2016/2017 at spatial resolutions of 10 m, 30 m, 50 m and 100 m. The R scripts as well as example data and the final woody cover product are freely available through zenodo.org for future usage.

Ethical consideration

I confirm that ethical clearance was not needed/required for the study.

Results

The woody cover map, derived from Sentinel-1 time series data acquired between 2016 and 2017, is shown in Figure 2. The estimations for woody cover range from 0% to 100%. The highest woody cover is reached in the northern part and the mopane-dominated landscapes on granite in the mid-west of the KNP as well as the mixed thorn and woodland areas in the southern part (30% – 50% and for some regions even higher than 50% on average). Lowest woody cover values were found in the open grasslands on basalt in the eastern part (below 20% on average). These broad-scale spatial patterns indicate that woody cover in KNP is significantly controlled by the underlying geology and the north-to-south rainfall gradient (Venter et al. 2003). Finer-scale contrasts in woody cover are also clearly discernable on the woody cover map. For example, the contrast between the lower woody cover gabbro areas within the higher woody cover granitic areas and the more open sodic areas as opposed to higher woody cover riparian zones are clearly visible.

FIGURE 2: Wall-to-wall woody cover map of the Kruger National Park with a spatial resolution of 10 m (right). The footprints of the light detection and ranging data used for training (red, Smit et al. 2016) and independent validation (orange, Milan et al. 2018) are shown in the woody cover map. Comparison of the woody cover information to high-resolution imagery from Google Earth for three selected regions (A–C) as well as a scatterplot visualising the independent validation of the woody cover map with a woody cover estimate based on light detection and ranging data from Milan et al. (2018) (lower left).

The RF approach revealed that the Sentinel-1 time series data from the dry season of 2016 were the most important predictors for deriving woody cover information. Furthermore, the cross-polarised scenes between July and September 2016 have been identified to have the highest ranking in the variable importance, followed by the co-polarised scenes of the same time span. This implies that Sentinel-1 predictors representing the dry season are key variables for retrieving woody cover in heterogeneous savanna ecosystems. This is most certainly caused by a lower impact of surface moisture on the C-Band backscatter as well as a better separation between woody plants and grasses during the dry season (Bucini et al. 2010; Heckel et al. 2020; Urbazaev et al. 2015).

The spatial cross-validation, which was performed on the training data, resulted in a root mean squared error (RMSE) of 22.8% for the product with a spatial resolution of 10 m, and improved with spatial averaging to 15.8% for 30 m, 14.8% for 50 m and 13.4% for 100 m as has been found in other radar studies such as Santoro et al. (2013). The spatial cross-validation revealed coefficient of determination values of R² = 0.3 (10 m), R² = 0.5 (30 m), R² = 0.6 (50 m) and R² = 0.6 (100 m), respectively. These error estimates are slightly lower in comparison with the accuracy of other woody cover products with spatial resolution ranging between 50 m and 120 m (Bucini et al. 2010; Higginbottom et al. 2018; Main et al. 2016; Urbazaev et al. 2015), with RMSE of around 10%. In comparison with this study, these products are lacking a spatially weighted cross-validation (except Higginbottom et al. 2018), resulting in modelling accuracies that do not take into account effects of spatial autocorrelation (Brenning 2012). An additional independent validation was carried out by using woody cover estimates based on the LiDAR data from Milan et al. (2018), confirming the results of the spatial cross-validation of the model (Figure 2 – scatterplot in the lower left for the spatial resolution of 10 m). Furthermore, the following RMSE and R² are found for the other spatial resolutions: 30 m (RMSE = 19%, R² = 0.5), 50 m (RMSE = 18%, R² = 0.6) and 100 m (RMSE = 16%, R² = 0.6).

Discussion

The presented workflow was applied to produce a wall-to-wall woody cover map in different spatial resolutions (e.g. 10 m, 30 m, 50 m and 100 m) for the KNP for 2016 and 2017, which are freely available via a data repository on zenodo.org. The user can choose between different spatial resolutions based on their needs. The woody cover maps provide insightful information for sustainable park management and future planning and conservation-related activities. The overestimation of woody cover below 20% is a typical phenomenon in radar retrievals (Mathieu et al. 2013; Santoro et al. 2011) caused by surface contributions to the signal, such as roughness. Underestimation for woody cover above 60% has been observed likewise in other studies (Bouvet et al. 2018) and is explained by saturation of the C-Band backscatter in dense canopies.

The presented workflow shows the capability of Sentinel-1 radar time series to retrieve fine-scale woody cover estimates in savanna ecosystems at regional scale. We found similar error estimates when comparing the accuracies to other woody cover products (e.g. Bucini et al. 2010; Higginbottom et al. 2018; Main et al. 2016; Urbazaev et al. 2015) at comparable spatial resolutions. This article describes the approach to derive woody cover from Sentinel-1 time series based on open-source software and open-access data to foster the reproducibility of the framework for individual research interests as well as to facilitate the potential of transferring the presented methodology to other savanna ecosystems.

It is worth mentioning that radar data are sensitive to environmental conditions, e.g., surface roughness and surface moisture, which need to be considered in future applications utilising this approach. Moreover, the definition of a height threshold for woody cover estimation, to exclude the non-woody vegetation component, might require an adjustment in other regions. We are currently focussing on comparing woody cover maps of different years to build up a woody cover change-monitoring system, where the mentioned environmental factors are analysed. In the near future, such a monitoring system will likely benefit from open access LiDAR data from the Global Ecosystem Dynamics Investigation (GEDI) (Coyle et al. 2015) on board of the International Space Station (ISS) and holds the potential to become an essential component of park management in the KNP and other national parks in South Africa.

Acknowledgements

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

The processing, analysis, interpretation of the used data and the preparation of the manuscript were completed by the first author M.U. K.H., C.B. and P.S. supported the analysis and interpretation of the data as well as proofreading of the manuscript. I.S. and T.S. supported during writing of the manuscript and proofreading. J.B. and C.S. supported the interpretation of the results and revised the manuscript.

The study was funded by the project South African Land Degradation Monitor (SALDi) sponsored by the German Federal Ministry of Education and Research (BMBF) in the framework of the Science Partnerships for the Assessment of Complex Earth System Processes (SPACES II) under the grant 01LL1701A. It was also partly funded by the European Union (EU) Horizon 2020 Research and Innovation Program under the Grant Agreement No. 640176 (BACI – Towards a Biosphere Atmosphere Change Index). Co-authors received funding from the ARS AfricaE and EMSAfrica research projects, which are part of the SPACES and SPACES II initiatives and funded by the BMBF (grant numbers 01LL1303D and 01LL1801D) and the project Munich Center for Machine Learning (MCML) from the BMBF (grant number 01IS18036A). Substantial scientific and logistical support in the KNP was received through our two SANParks projects ARS Africae EO (earth observation) (ref. no. SCHCC1235) and CO-LD-EMS (ref. no. BAAJ1519).

The Carnegie Airborne Observatory (CAO) LiDAR data collection from Smit et al. (2016) was funded by the Andrew Mellon Foundation. The CAO was supported by the Avatar Alliance Foundation, John D. and Catherine T. MacArthur Foundation, Andrew Mellon Foundation, Mary Anne Nyburg Baker and G. Leonard Baker Jr., and William R. Hearst III. The LiDAR data for independent validation of the woody cover map were provided by G. Heritage (University of Salford, UK) and D. Milan (University of Hull, School of Environmental Sciences, UK).

The programming code (statistical software R), Sentinel-1 and LiDAR training data as well as data sets for model prediction of three selected regions in the KNP are freely available via the following link: https://doi.org/10.5281/zenodo.3728186.

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 authors.

References