*Article* **The MODIS Global Vegetation Fractional Cover Product 2001–2018: Characteristics of Vegetation Fractional Cover in Grasslands and Savanna Woodlands**

#### **Michael J. Hill 1,2,\* and Juan P. Guerschman <sup>1</sup>**


Received: 24 December 2019; Accepted: 24 January 2020; Published: 28 January 2020

**Abstract:** Vegetation Fractional Cover (VFC) is an important global indicator of land cover change, land use practice and landscape, and ecosystem function. In this study, we present the Global Vegetation Fractional Cover Product (GVFCP) and explore the levels and trends in VFC across World Grassland Type (WGT) Ecoregions considering variation associated with Global Livestock Production Systems (GLPS). Long-term average levels and trends in fractional cover of photosynthetic vegetation (*F*PV), non-photosynthetic vegetation (*F*NPV), and bare soil (*F*BS) are mapped, and variation among GLPS types within WGT Divisions and Ecoregions is explored. Analysis also focused on the savanna-woodland WGT Formations. Many WGT Divisions showed wide variation in long-term average VFC and trends in VFC across GLPS types. Results showed large areas of many ecoregions experiencing significant positive and negative trends in VFC. East Africa, Patagonia, and the Mitchell Grasslands of Australia exhibited large areas of negative trends in *F*NPV and positive trends *F*BS. These trends may reflect interactions between extended drought, heavy livestock utilization, expanded agriculture, and other land use changes. Compared to previous studies, explicit measurement of *F*NPV revealed interesting additional information about vegetation cover and trends in many ecoregions. The Australian and Global products are available via the GEOGLAM RAPP (Group on Earth Observations Global Agricultural Monitoring Rangeland and Pasture Productivity) website, and the scientific community is encouraged to utilize the data and contribute to improved validation.

**Keywords:** vegetation; grassland; savanna; fractional cover; trend; ecoregion; bare soil; livestock; production systems

#### **1. Introduction**

Vegetation Fractional Cover (VFC) is an important global indicator of land cover change, land use practice, and landscape and ecosystem function [1,2]. Seasonal dynamics and long-term trends in the fractional cover of photosynthetic vegetation (*F*PV), non-photosynthetic vegetation (*F*NPV), and bare soil (*F*BS) may identify changes in cropping cycles, impacts of livestock grazing, clearing or planting of woody vegetation, and ecosystem responses to climate shifts. The Australian Vegetation Fractional Cover Product (AVFCP) is derived from the 500 m MODIS (MODerate Resolution Imaging Spectroradiometer) NBAR (Nadir BRDF-adjusted Reflectance) product (MCD43A4) and has been comprehensively documented, validated, and improved over several versions [3–6]. In this study, we present the Global Vegetation Fractional Cover Product (GVFCP) based on the same methodology. Both the Australian and Global products are available through the GEOGLAM RAPP (Group on Earth Observations Global Agricultural Monitoring Rangeland and Pasture Productivity) website (https://map.geo-rapp.org). A list of abbreviations appears in Table A1.

A key feature of the GVFCP is the explicit sensitivity to *F*NPV enabling better discrimination of dry cellulosic vegetation and bare soil fractions than is possible with the family of greenness indices. For this reason, the GVFCP is particularly important for monitoring of the non-forested and non-desert global biomes, woodlands, savannas, grassland and shrublands, where the overstory canopy is open, and understory dynamics involve the relationships between cellulosic herbaceous vegetation and/or senescent arboreal leaf litter and bare soil. The first version of the AVFCP [3] used the relationship between the Normalised Difference Vegetation Index (NDVI) and the Short Wave Infrared Ratio (SWIR32; a ratio of the MODIS 2130 nm and MODIS 1640 nm bands) in a linear unmixing method to approximate the relationship between NDVI and the Cellulose Absorption Index (CAI) developed and demonstrated by [7–11]. An improved version of the AVFCP utilizes all seven MODIS NBAR bands plus band transformations and interaction terms in a sophisticated unmixing algorithm [4]. Recent recalibration of the AVFCP with MODIS Collection 6 inputs and an updated validation data set established improvements in the accuracy of the retrieved fractions of photosynthetic vegetation (*F*PV RMSE 0.112), non-photosynthetic vegetation (*F*NPV RMSE 0.162), and bare soil (*F*BS RMSE 0.130) [5].

Based on Australian studies [12–14], the AVFCP provides a consistent and verifiable estimate of the cover fractions for green canopy, senescent vegetation and surface litter, and the visible soil surface. Given the diversity of natural ecosystems across Australia that include tropical rainforest, tropical and temperate savanna, temperate grasslands, semi-arid and arid shrublands and grasslands, temperate eucalypt forests, and temperate rainforest, and the extent of agricultural and pastoral lands across Australia, it is reasonable to propose that the GVFCP should perform with similar levels of uncertainty across the diversity of global land covers and land uses, especially in the non-forested and non-desert biomes.

Many global analyses have focused on tree cover change [15–17] and land cover change [18]. The recent comprehensive study of [19] partitions cover into tall vegetation, short vegetation, and bare soil, and the greening analysis of [20] looks at trends in the annual values of the MODIS leaf area index product [21]. In addition, other studies have focused on the total cover and soil or bare fractions [22,23]. There are several global tree cover [24,25] and tree cover change [15] products at multiple resolutions and many FAPAR (Fraction of Absorbed Photosynthetically Active Radiation) [26,27] and green canopy cover or leaf area index products [28–30]. The global cover change study by [19] identifies changes in short vegetation (which includes shrubs under 5 m in height) and bare ground. The cover fractions are mapped from peak growing season canopy [19] meaning that deciduous savanna trees and shrubs are probably fully foliated, and understory grasses are actively growing. By contrast, the GVFCP at full temporal and spatial resolution follows the phenology of overstory and understory throughout the seasons and explicitly retrieves *F*NPV representing senescent understory.

Studies focusing on global scale cover, dynamics, and change in savanna landscapes are less common (e.g., [31,32]), but a study based on persistence of GIMMS (Global Inventory Modeling and Mapping Studies) NDVI signals indicated larger areas of increased than decreased vegetation persistence [33]. However, there has been considerable discussion on the nature of stable states, methods of detection, and drivers of change in savannas (e.g., [34–37]). Due to the specific methodology used to derive it, the GVFCP provides a unique measure of global changes in vegetation fractional cover especially important for woodlands, savannas, grasslands and shrublands which are subject to heavy utilization and conversion for human food production and where cellulosic herbaceous cover is important.

It is instructive to provide context for the behaviour of remote sensing products using global datasets that are wholly or partially independent of remote sensing, and describe global vegetation in terms relevant to conservation, ecosystem function, and productivity. The Terrestrial Ecoregions of the World (TEoW; [38]) and the subset of these used to define World Grassland Types (WGT [39]) provides an effective framework within which to view the patterns of VFC for the grassy ecoregions of the world. Since these grassy ecoregions are subject to major utilization by humans for extensive and intensive

agriculture, the Global Livestock Production Systems (GLPS) classification [40] provides a means to sample variation within WGT Ecoregions due to differences in utilisation between production systems.

In this study we have three objectives:


#### **2. Materials and Methods**

The analysis undertaken here provides an overview of global VFC in map form; records benchmark average levels of VFC across production systems for Formation, Division, and Ecoregion areas of the grassy biomes of the world; and documents long-term trends in VFC from 2001–2018 for these areas. The analysis is designed to provide an overview framework and baseline for more detailed future studies.

#### *2.1. Data*

#### 2.1.1. Global Vegetation Fractional Cover Product

The GVFCP estimates the fractions of photosynthetic (green) vegetation (PV), non-photosynthetic (non-green) vegetation (NPV), and bare soil (BS) within the pixel using a spectral linear unmixing method [4,5]. The base version of the GVFCP, available on the GEOGLAM RAPP-MAPP website (https://map.geo-rapp.org) is calculated using reflectance values from the MODIS MCD43A4 product which is a rolling daily 16-day composite weighted to retrieve the value on the ninth day of the 16-day period [41]. The GVFCP utilizes every eighth day from the MCD43A4 product starting with day one of each year with an adjustment for leap years proving 46 dates in each calendar year. The data are stored in a geographic latitude-longitude raster with a WGS84 spheroid. A second monthly composite product (also available on the website) is created by aggregating all values falling within each calendar month using a medoid function (a multidimensional median) following [42]. Pixels flagged as water, snow, or low quality in the MCD43A4 product are ignored in the calculation of the GVFCP. The current product is the result of a decade of development and the most recent version derived from MODIS Collection 6 NBAR data was calculated and validated using an updated and expanded database of 3022 ground measurements of fractional cover across Australia [5] and had a Root Mean Square Error (RMSE) of 11.3%, 16.1%, and 14.7% in the PV, NPV, and BS fractions, respectively. The GVFCP is derived from spectral unmixing of all seven MODIS optical bands from the 500 m MODIS Nadir BRDF-adjusted Reflectance Product (NBAR, MCD43A4 Collection 6; see Appendix A Table A1 for abbreviations) using previously described methods [4]. The product consists of four derived layers corresponding to the decimal values of cover fractions *F*PV, *F*NPV, *F*BS, and the residual from the spectral unmixing in reflectance values (RNORM). The GVFCP is calculated from the daily NBAR data product producing a phased 16-day composite every eight days. For this overview study, the monthly version of the product was used and resampled to 0.05◦ resolution (approximately 5 km2) by averaging.

#### 2.1.2. Terrestrial Ecoregions of the World

The Terrestrial Ecoregions of the World (TEoW) provide a biogeographic realization of terrestrial biodiversity [38] (Figure 1a). They represent land areas that share species, ecosystems, seasonal dynamics, and environmental conditions, and therefore, in part, indicate boundaries for intrinsic capability or risk for human uses. This dataset is used to meet the first objective of this study by providing general context for the description of the GVFCP at TEoW Realm level. Since an important attribute of the GVFCP is the derivation of *F*NPV and *F*BS, additional datasets that focus on the characteristics and human utilization of the subset of ecoregions either partially or totally dominated by herbaceous vegetation were employed to examine behaviour in the grassy regions of the world.

#### 2.1.3. World Grassland Types

The World Grassland Types (WGT; [39]) represent a subset of terrestrial ecoregions that were combined with the International Vegetation Classification (IVC; [43]) to spatially define the major grassy vegetation systems of the world (Figure 1b; see Table A1 for abbreviations). They describe 75% of the IVC grasslands using ecoregions. The remaining IVC grasslands were absent from the mapping due to their fine scale distribution [39]. The WGT Ecoregions cover an area of 35.5 M km2. The Tropical Lowland Shrubland, Grassland and Savannas (TLSGS), Tropical Montane Shrubland, Grassland and Savanna (TMSGS), and the Warm Semi-Desert Scrub and Grassland (WSDSG) WGT Formations were selected to represent the savanna, woodland and scrub grasslands (SWSG) which are the focus of this Special Issue. These SWSGs occupy 46.9% of the WGT area (16.7 M km2) of which 9.5 M km2 are in the northern hemisphere in Africa and neo-tropical South America, and 7.16 M km<sup>2</sup> are in the southern hemisphere in Africa, Australia, and South America. The northern hemisphere contains 7.0 M km2 of TLSGS, 0.3 M km<sup>2</sup> of TMSGS, and 2.2 M km2 of WSDSG. The southern hemisphere contains 6.1 M km2 of TLSGS, 0.25 M km2 of TMSGS, and 0.82 M km<sup>2</sup> of WSDSG.

#### 2.1.4. Global Livestock Production Systems

The Global Livestock Production Systems (GLPS) data describes the terrestrial land surface using livestock density maps, crop type maps, land cover, and information on management practices [40] (Figure 1c). The classification describes 12 production systems combining rangelands (L—grazing only), mixed rainfed and mixed irrigated production (M—grazing and cropping) with hyper-arid (Y), arid (A), humid (H), and temperate (T) environments at a spatial resolution of 0.00833◦ (approximately 1 km2; Figure 1c). These are: Rangelands Hyper-arid (LGY); Rangelands Arid (LGA); Rangelands Humid (LGH); Rangelands Temperate (LGT); Mixed Rainfed Hyper-arid (MRY); Mixed Rainfed Arid (MRA); Mixed Rainfed Humid (MRH); Mixed Rainfed Temperate (MRT); Mixed Irrigated Hyper-arid (MIY); Mixed Irrigated Arid (MIA); Mixed Irrigated Humid (MIH); and Mixed Irrigated Temperate (MIT). The remaining land areas are classified as Urban or Other (such as forest, ice, rock, etc.). These classes describe land use with each WGT Ecoregion. The WGT Ecoregions are dominated by seven GLPS types with the six major production systems being: LGA (36.8%); LGT (17.6%); MRA (8.4%); MRT (14.5%); LGH (4.7%); and MRH (5.6%) and the remaining area being Other.

#### *2.2. Analysis Method*

#### 2.2.1. Trend Analysis for World Grassland Types

The long-term trend from 2001–2018 in the monthly *F*PV, *F*NPV, and *F*BS was evaluated for the WGT Ecoregions using the Mann–Kendall test [44,45] with a custom-built function in IDL. This non-parametric test detects monotonic trends in time series data where the least squares regression approach is invalid due the autocorrelation in time series data derived from remote sensing. We used the Mann–Kendall approach for trend analysis on the full monthly data set since these data provided more observations than annual means by a factor of 12. The time series was smoothed by the medoid procedure, spatially averaged from 500 m to 5 km pixels by the pixel aggregation, and completely continuous in the WGT Ecoregions. As a result, we did not consider it necessary to apply the seasonal Kendall method [46]. The significance (p-value) of the slope in each pixel was recorded. Global maps of pixels within the WGT Ecoregions where trends were significant at p < 0.1 were produced. The same spatial patterns and trends were evident at p < 0.05 but patterns were clearer and more contiguous at p < 0.1 similar to [20].

**Figure 1.** (**a**) Terrestrial Ecoregions of the World (TEoW) at Realm level. (**b**) World Grassland Type Ecoregions (WGT; legend shows formations). (**c**) Global Livestock Production Systems (GLPS). System codes as follows: Rangelands Hyper-arid (LGY); Rangelands Arid (LGA); Rangelands Humid (LGH); Rangelands Temperate (LGT); Mixed Rainfed Hyper-arid (MRY); Mixed Rainfed Arid (MRA); Mixed Rainfed Humid (MRH); Mixed Rainfed Temperate (MRT); Mixed Irrigated Hyper-arid (MIY); Mixed Irrigated Arid (MIA); Mixed Irrigated Humid (MIH); Mixed Irrigated Temperate (MIT).

#### 2.2.2. Summarizing Levels and Trends across World Grassland Types and Savanna-Woodlands

The levels and trends in VFC were examined at increasing levels of detail across the WGT. The WGT and GLPS layers were combined to create a composite layer with individual classes for each GLPS-ecoregion combination. Class mean and standard deviation values for long-term average and trend in *F*PV, *F*NPV, and *F*BS were extracted using ArcGIS®zonalstats. For trends, zonal statistics were collected only for slope values > 0.1 or < −0.1 using a mask that restricted extractions to only those pixels with significance at p < 0.1. The areas of significant trend by Division and Ecoregion were calculated on geographic grids using a weighted area algorithm to correct for pixel distortion with latitude away from the equator. The variation across GLPS types within WGT Divisions was displayed using the fence box plot. The fence box analysis plots the median and divides the observed values into quartiles with first and third quartiles defining the hinges, and the upper and lower fences being defined by multiplying the interquartile range by 1.5. The whiskers represent the minimum and maximum observations falling inside the fences and points falling outside this range represent outliers. The fence-box plots at Division level therefore show the range in mean values among GLPS types within that Division (whether it contains one or many Ecoregions) and the range in the standard deviation values from each GLPS type within that Division. Examination of the SWSG was carried out by the same methods at individual ecoregion level, and GLPS class behaviour was explored within selected ecoregions exhibiting large areas of significant trends.

#### **3. Results**

#### *3.1. Global Patterns of Average VFC*

Global geographical distribution of long-term average *F*PV and *F*BS exhibit the typical associations with forests and deserts and the gradients in between (Figure 2). The distribution of long-term average *F*NPV tends to highlight the large grasslands of North America, South America, central Eurasia, Tibet, Mongolia, and the savannas of southern Africa and Australia (Figure 2). Based on TEoW definitions, there are 19.5 M km2 of tropical and subtropical grasslands, shrubland, and savannas representing 13.26% of the terrestrial land surface (Table 1). In addition, the global coverage of Temperate Grasslands and Shrublands, Montane Grasslands and Shrublands, and Flooded Grassland and Savannas brings the total area of grassy biomes to 35.439 M km<sup>2</sup> representing 24.05% of the terrestrial land surface. Across these grassy biomes, the Montane Grasslands and Shrublands have lower average *F*PV and higher average *F*BS than the other grassy biomes (Table 1). Temperate and Montane Grasslands and Shrublands have higher average *F*NPV than the Tropical and Subtropical Grasslands, Savannas and Shrublands and the Flooded Grassland and Savannas. The global latitudinal variation in average FPV, FNPV, and FBS exhibits typical peaks in FPV for tropical forests between 20◦ S and 20◦ N and boreal forests above 50◦ N and in FBS in northern and southern hemisphere arid lands between latitudes 20 and 30◦ (Figure 3). However, there is a notable rise in FNPV between 20 and 40◦ N which may reflect the large expanses of steppe environments either limited by precipitation or converted to cereal cropping. With much less land at mid- to low-latitudes, the dynamics between the cover fractions in the southern hemisphere exhibit distinct narrow latitude zones of oscillation between higher and lower levels of *F*PV and *F*NPV between 30 and 45◦.

**Figure 2.** Global average fractional cover from 2001–2018. The ternary plot shows the correspondence between colour and the values fractions of photosynthetic vegetation (*F*PV), non-photosynthetic vegetation (*F*NPV) and bare soil (*F*BS). Each ternary axis represents colours corresponding to two-factor mixtures, the greater colour intensity indicates more dominance of a single cover type.

**Table 1.** Mean and standard deviation of fractional cover (%) from 2001–2018 across the Realms of the Terrestrial Ecoregions of the World (TEoW). Note that overall means do not sum to exactly 100% due to temporal and spatial averaging effects.


**Figure 3.** Global latitudinal variation in average fractional cover of *F*PV, *F*NPV, and *F*BS.

#### *3.2. Variation in Vegetation Fractional Cover within Global Livestock Production Systems*

The GLPS types exhibit wide variation in mean *F*PV across ecoregions with large ranges between upper and lower hinges and fences. There was less variation in mean *F*NPV and *F*BS, although the latter exhibited wide variation in the upper and lower fence values and more extreme (outlier) ecoregions than *F*PV and *F*NPV (Figure 4). Median and hinge values for *F*PV increase from hyper-arid (Y) to arid (A) to humid (H) for both grazing (L) and mixed (M) land uses, while median and hinge values for temperate (T) systems were similar to arid systems for grazing and mixed land uses. The average *F*BS also exhibits a wide range of values across ecoregions within a GLPS type with values > 70% and ≤ 5% for LGA and LGT suggesting differences in interactions between the vegetation and livestock type and management. Average *F*BS is notably low across ecoregions for both LGH and MRH. Average *F*NPV exhibits a lower range across ecoregions within a GLPS type, and between GLPS types than *F*PV, with median values of ecoregion means across all GLPS types lying between 25%–45%. A gross indication

of the impact of the GLPS on vegetation cover can be gained from comparing the medians and ranges of the Other class with those of the production systems. For *F*PV, the median value is higher and the range is wider than all GLPS types except for the LGH, MRH, and MIH types where tropical pasture and cropland conversion results in higher productivity. Similarly, *F*BS for the Other class is lower than all GLPS types except for the LGH, MRH, and MIH types.

**Figure 4.** Fence box plots showing themedians, upper and lower hinges and upper and lower fences ofmean long-term fractional cover of ecoregions over within Global Livestock Production System (GLPS) classes.

#### *3.3. Variation in Average Vegetation Fractional Cover in World Grassland Type Divisions*

The WGT Ecoregions are organized by Formations and Divisions (Table A2) with the Formations representing a broad structural type, and the Divisions representing the regional geographical representations of these Formations. The WGT Divisions exhibit a wide range of long-term average levels of *F*PV, *F*NPV, and *F*BS and substantial variation in the levels of the upper and lower hinges and fences (Figure 5) There is also major variation in the standard deviations of VFC among GLPS classes with a WGT Division. Across the WGT, *F*PV exhibits greater variation in upper and lower hinge values than *F*NPV and *F*BS. The Divisions within the TLSGS and TMSGS Formations (lowland and montane savanna woodlands) exhibit relatively high average *F*PV and low average *F*BS except for the NSSDSG Division which covers the Sahelian region of northern Africa. Among the other Divisions, the MBDG, MSACSDSG, PCSDSG, TACSDSG, and AWSDSG within the CSDSG (semi-desert) and MSGFM (Mediterranean) Formations are

notable for very low average *F*PV, median levels of *F*NPV close to 50%, and relatively high levels of *F*BS. The WEGS Division in western Eurasia, and the montane AMGS, and the alpine CAASFMG Divisions exhibit wide variation in average *F*PV but WEGS exhibits a lower median and much lower range between hinge values for standard deviation than AMGS and CAASFMG.

**Figure 5.** Fence box plots showing the medians, hinges, and inner and outer fences of average long-term vegetation fractional cover of ecoregions within World Grassland Type Divisions (listed in Table A2. (**a**) *F*PV; (**b**) *F*NPV; and (**c**) *F*BS.

#### *3.4. Trends in Vegetation Fractional Cover*

#### 3.4.1. World Grassland Types

The spatial distribution of pixels exhibiting significant trends at p < 0.1 is shown in Figure 6. Across the WGT Ecoregions, the extent of significant trends was greater for *F*BS than for the vegetation fractions since vegetation cover trends were split between change in *F*PV and change in *F*NPV and these vary in different ways for different geographies. There are significant positive and negative long-term trends in *F*PV, *F*NPV, and across the WGT Ecoregions (Figure 7). The WGT Ecoregions are displayed over the global map of countries, showing that the significant trends affect many African and Eurasian countries whilst the USA, Brazil, Argentina, and Australia with their large areas of grassland and savanna experience change over very large areas. There are positive trends in *F*PV in the northern great plains of North America, parts of China, and parts of eastern and southern Africa (Figure 7a). There are notable negative trends in *F*NPV in Patagonia, across Sahelian and Sudanian Africa, in Mongolia and on the Tibetan Plateau, and in the Mitchell Grasslands of Northern Australia (Figure 7b). However, the largest areas of significant trends occur for *F*BS with negative trends across the Great Plains of North America, across central China, across northern and southern Australia, and areas of positive trends in Sahelian, Sudanian and East Africa, western Mongolia, Ukraine and southern Russia, and the Mitchell Grasslands of Northern Australia (Figure 7c).

Examination of the variation due to GLPS in trends for WGT Divisions show that although the median values within Divisions seldom exceed ± 0.2%, the lower hinge values are as low as −0.5% and the upper hinge values are as high as 0.4% for *F*PV, as low as −0.4%, as high as 0.4% for *F*NPV, and as high as 0.4% for *F*BS and as low as −0.5% for *F*BS (Figure 8). It is particularly notable that certain Divisions (CAASMG, EECSDSG, WECSDSG, EEGS, NAGS, and WEGS) contain individual GLPS-Ecoregion combinations with positive and negative trends in *F*NPV as great as ± 0.5% and in *F*BS as low as −1.0% and as high as 0.5% that sit well beyond the upper and lower fences (outliers). These Divisions are in the northern hemisphere and are in the alpine ASFMG, semi-desert CSDSG, Mediterranean MSGFM, and temperate TGMS Formations suggesting major land use effects between GLPS types. Among the Divisions in the TLSGS and TMSGS Formations (savanna woodlands), most northern hemisphere Divisions show small increasing trends in *F*PV and relatively insignificant changes in *F*NPV and *F*BS. However, in the southern hemisphere the lower hinge and fence values for *F*PV, upper and lower hinge and fence values for *F*NPV, and upper hinge and fence values for *F*BS in the mopane and bushveld savannas (MS), and the montane WCAMWS, AMGS, and BPMSG identify GLPS classes with negative trends in *F*PV and positive trends in *F*NPV and *F*BS (Figure 8).

The fence box analysis describes the range in variation in VFC trends across GLPS types within WGT Divisions, however the Divisions vary in size and this may mask localised strong trends in large Divisions (such as the Sahelian Acacia Savanna which is a single enormous Ecoregion) and emphasize trends in small Divisions with more uniform climate, terrain, and edaphic features (such as the Ethiopian Montane Grasslands and Woodlands). Hence, it is important to define the area and area percentages where significant positive and negative trends in VFC are occurring (Figure 9). The area analysis identifies eight Divisions where areas in excess of 200,000 km<sup>2</sup> show significant positive trends in *F*PV (EBGS, EEGS, GPGS, WCAMWS, EAXSG, AMS, ATS, and BPLSGS) and seven Divisions where areas in excess of 100,000 km<sup>2</sup> show significant negative trends in *F*PV (EBGS, EEGS, GPGS, EAXSG, PGS, BPLSGS, and AWSDSG; Figure 9a). Note that five Divisions have both large areas of positive and negative trends. By contrast, some of the smallest Divisions in area, have significant positive or negative trends in *F*PV over a large proportion of their area: greater than 60% of IMW and AAFSMG and about 50% of NGMM have positive trends in *F*PV, while greater than 20% of CGM, CVFMWMS, PGS, AMGS, and MMGS have significant negative trends in *F*PV.

**Figure 6.** (**a**) The spatial distribution of pixels exhibiting significant long-term trends at p < 0.1 from 2001–2018 for *F*PV; (**b**) The spatial distribution of pixels exhibiting significant long-term trends at p < 0.1 from 2001–2018 for *F*NPV; and (**c**) The spatial distribution of pixels exhibiting significant long-term trends at p < 0.1 from 2001–2018 for *F*BS.

**Figure 7.** (**a**) Geographical pattern of significant (p < 0.1) positive or negative long-term trajectory of vegetation fractional cover for: *F*PV; (**b**) Geographical pattern of significant (p < 0.1) positive or negative long-term trajectory of vegetation fractional cover for *F*NPV; and (**c**) Geographical pattern of significant (p < 0.1) positive or negative long-term trajectory of vegetation fractional cover for *F*BS.

**Figure 8.** Fence box plots showing the median, upper and lower hinges and upper and lower fences of the average slope of the regression trend line from 2001–2018 for fractional cover within WGT vegetation divisions (listed in Table A2). (**a**) *F*PV; (**b**) *F*NPV; and (**c**) *F*BS.

**Figure 9.** *Cont*.

**Figure 9.** Area (km2) and percentage of area of WGT Divisions exhibiting significant positive or negative trends in (**a**) *F*PV, (**b**) *F*NPV, and (**c**) *F*BS. Key to Formation and Division acronyms in Table A2.

There are also large areas of significant change in *F*NPV with 13 Divisions having areas in excess of 100,000 km<sup>2</sup> showing significant positive trends in *F*NPV and seven Divisions having areas in excess of 200,000 km2 show significant negative trends in *F*NPV (Figure 9b). As with *F*PV, there are several Divisions with large areas of positive and negative trends. In addition, greater than 20% of the area of MBDG, NAWDSG, PGS, and AWSDSG show positive trends and greater than 20% of the area of MBDG, NSSDSG, EAXSG, AASFMG, PCSDSG, BPFMWMS, NZGS, and NGMM show a negative trend in FNPV. Percentage areas of significant trends tended to be greater for FBS than for FNPV and FPV, since FBS reflects combined changes in FPV and FNPV (Figure 9c). There were 15 Divisions with greater than 200,000 km2 showing negative trends in FBS.

Only five Divisions had areas greater than 200,000 km<sup>2</sup> showing positive trends in FBS. However, 21 Divisions had greater than 20% of their area showing positive trends in FBS (WECSDSG, CVFMWMS, CGM, MBDG, WEGS, NSSDSG, SSDS, EAXBG, MSACSDSG, PCSDSG, TACSGSG, AMS, NZGS, PGS, SAMG, ATS, ESADSW, MS, AMGS, BPMSG, AWSDSG, MMGS). Large areas of both positive and negative trends in FBS only occurred in three Divisions (NSSDSG, SSDS, AWSDSG). There 19 Divisions with greater than 20% of their area showing negative trends in FBS.

#### 3.4.2. Savanna Woodland and Scrub Grasslands

Within the SWSG (see Table A3 for a list of ecoregion names and associated Divisions), individual ecoregions exhibit wide variation in the percentage of their area with significant positive and negative trends in *F*PV, *F*NPV, and *F*BS (Figure 10). In the northern hemisphere parts of Africa, greater than 50% of the area of the Saharan flooded grasslands and the Kinabalu montane alpine meadows show significant positive trends in *F*PV (Figure 10a). By contrast, more than 20% of the areas of the Masai Xeric Grasslands and Shrublands and the Northern and Southern Acacia-Commiphora Bushlands show significant negative trends in *F*PV. In North and South America, the Chihuahan desert, Llanos, Guianian savanna, Rio Negro campinarana, and Pantepui have more than 20% of their areas exhibiting significant positive trends in *F*PV. In the southern hemisphere, the Australian savannas and scrublands and South American lowland and montane savannas exhibit significant positive trends in *F*PV, with 50% or more of the area of the Brigalow Tropical Savanna, Central Range Sub-Alpine Grasslands, Cordillera de Merida píramo, and Northern Andean píramo exhibiting positive trends *F*PV. However, in southern Africa, the Victoria Basin forest-savanna mosaic, Madagascan ericoid thicket, Angolan montane forest-grassland mosaic, and the Rwenzori-Virunga montane moorlands all have greater than 20% of their areas exhibiting negative trends in *F*PV.

**Figure 10.** *Cont*.

**Figure 10.** Area (%) of savanna, woodland and shrubland WGT Ecoregions showing significant monotonic trends in (**a**) *F*PV; (**b**) *F*NPV; and (**c**) *F*BS.

The Masai Xeric Grasslands and Shrublands and the Northern Acacia-Commiphora Bushlands also had large areas of negative trends in *F*NPV, while the Saharan Flooded Grasslands, Victoria Basin forest-savanna mosaic, Madagascan ericoid thicket, Angolan montane forest-grassland mosaic, and the Rwenzori-Virunga montane moorlands had large areas of positive trends in *F*NPV (Figure 10b). However, there are several ecoregions showing large areas of significant trends in *F*NPV not strongly associated with area patterns for *F*PV. Large areas of positive trends in *F*NPV occur in the Chihuahan desert or North America and the Great Sandy Desert of Australia, and large areas of negative trends in *F*NPV occur in the Mitchell Grasslands and Central Range Sub-Alpine Grasslands of Australia (Figure 10b).

There are very large areas of significant positive and negative trends in *F*BS in many savanna ecoregions (Figure 10c). Across much of the Australian and South America savannas much larger percentages (> 30%) of the ecoregion areas exhibit negative trends in *F*BS than exhibit positive trends; the Brigalow savanna has a very large area of negative trends in *F*BS. However, the Mitchell Grasslands are the exception with around 50% of the area of this very large ecoregion exhibiting a positive trend in *F*BS. There are some sharp distinctions between African regions in the northern hemisphere, with Masai Xeric Grasslands and Shrublands and the Northern, Southern and Somali Acacia-Commiphora Bushlands having greater than 40% of their areas exhibiting significant positive trends in *F*BS while West African and montane East African ecoregions show much smaller and balanced areas between positive and negative trends. About 60% of the Chihuahan desert ecoregion exhibits a negative trend in *F*BS that corresponds to the areas of positive trends in *F*NPV and *F*PV. Across the SWSG Formations, 3.13 M km<sup>2</sup> exhibited significant positive trends in *F*BS and 2.79 M km2 exhibited negative trends in *F*NPV (Table 2). These areas represented 18.1% and 16.8%, respectively of the total area of SWSG Formations (Table 3). However, 2.31 M km<sup>2</sup> exhibited significant positive trends in *F*PV, 3.89 M km2 exhibited negative trends in *F*BS, and 1.69 M km<sup>2</sup> exhibited positive trends in *F*NPV representing 13.8%, 10.1%, and 23.3% of the total area of SWSG Formations. The TLSGS Formation contained large areas where trends exceeded <sup>±</sup> 0.1 FC units yr−<sup>1</sup> but were not significant at p <sup>&</sup>lt; 0.1.


**Table 2.** Area of Savanna Woodland and Scrub Grassland (SWSG) Divisions exhibiting significant trends in Vegetation Fractional Cover with magnitudes <sup>&</sup>gt; 0.1 or <sup>&</sup>lt; <sup>−</sup>0.1 % yr <sup>−</sup>1.

**Table 3.** Percentage of area of SWSG Divisions exhibiting significant trends in VFC <sup>&</sup>gt; 0.1 or <sup>&</sup>lt; <sup>−</sup>0.1 % yr <sup>−</sup>1.


#### *3.5. Variation in Vegetation Fractional Cover in Example Ecoregions*

While extensive examination of VFC behaviours within ecoregions is beyond the scope of this study, it is worthwhile to examine several notable regional trends in *F*PV, *F*NPV, and *F*BS and compare these across GLPS types. The east African Acacia-Commiphora Bushlands and Masai Xeric Grassland and Shrublands contained large areas of significant trends, but relativities among GLPS types were different (Figure 11a). In the Northern and Somali ecoregions, the LGH and MRH systems showed positive trends in *F*PV, while the MRA system showed a positive trend in *F*BS. By contrast, all GLPS systems across the Southern Acacia-Commiphora Bushland exhibited positive trends in *F*BS and negative trends in *F*PV and *F*NPV. In the Masai Xeric Grassland and Shrubland, the main features are a positive trend in *F*BS for LGA and a negative trend in *F*PV for the Other class.

The Chihuahan Desert exhibits small to moderate positive trends in *F*PV and *F*NPV across all GLPS types and large negative trends in *F*BS especial for LGA, LGT, MRA, and MRT (Figure 11b). In the Mitchell Grasslands, LGA and MRA systems dominate and both show positive trends in *F*BS and negative trends in *F*NPV. In the Llanos, positive trends in *F*PV across LGA, LGH, MRA, and MRH types are matched by small negative trends in both *F*NPV and *F*BS.

**Figure 11.** Trends in vegetation fractional cover among Global Livestock Production Systems (GLPS) areas within selected Ecoregions in (**a**) Africa; and (**b**) North America, South America, and Australia. The histogram indicates the average and the whiskers indicate the standard deviation among GLPS pixels within the Ecoregion.

#### **4. Discussion**

This study has introduced a unique remotely-sensed global vegetation product derived from MODIS reflectance data, the GVFCP, which is available as both an eight day and monthly 500 m resolution product. We have mapped global averages for *F*PV, *F*NPV, and *F*BS for data derived from MODIS for 2001-2018. The study then explored the behaviour of *F*PV, *F*NPV, and *F*BS across GLPS types within the WGT Divisions and Ecoregions including analysis of long-term trends in vegetation fractional cover with a focus on the SWSG Formations. The analysis illustrated the variation in average long-term levels and long-term trends of VFC associated with different GLPS types, with vegetation Formations and with geographically separate Divisions within the same Formation. The variation observed emphasized the strong interaction between land use management and vegetation characteristics. Although other products have previously provided global quantitative retrievals of *F*PV and *F*BS,

the GVFCP is new and important since it provides an explicit, spectrally functional measure of *F*NPV. This provides effective discrimination of both global regions with high average levels of *F*NPV and high seasonal fluctuations in *F*NPV, such as in the grassy understorey or tropical savannas, and the massive leaf litter associated with the deciduous forests of the eastern USA. It also highlights global regions where certain GLPS types or cropland expansion produce large amounts of stubble and crop residue such as southern and central China, south-eastern and south-western Australia, and the US corn belt.

The study has focused more detailed evaluation of levels and trends in *F*PV, *F*NPV, and *F*BS on the grassy biomes of the world represented by WGT Ecoregions, and for this Special Issue on the SWSG Formations within the WGT. There are large differences between average levels and long-term trends in *F*PV, *F*NPV, and *F*BS across the WGT. Fence box graphical analysis identified specific Divisions where average *F*BS was low, where average *F*PV was high, and where the upper and lower hinges, fences, and outlier points indicated massive variation in average levels within regional Formations (i.e., in different hemispheres and continents) indicative of possible climate or land use effects. Recent studies have identified global "greening" with an increase in the *F*PV attributable to developments in India and China [20]. Our analysis agrees with [20] in terms of the "greening" of north-eastern China, the US Great Plains, and eastern Australia. It also shows that the SWSG Formations, much of which experience an LGH production system, and contain all the major grassy savanna, woodland and scrub ecoregions, tend to exhibit more areas of positive to neutral trends in *F*PV than negative trends in *F*NPV or positive trends in *F*BS. Although it is difficult to make direct comparisons with the study of [19], there is broad agreement on the magnitude of areas (SWSG here, tropical dry forest and tropical shrubland in [19] experiencing negative trends in bare ground, and a positive trend in short vegetation (SV) in that study, most of the *F*PV and *F*NPV in SWSG in this study). These tropical humid systems either have reliable seasonal rainfall to support land cover changes (e.g., South American tropical savannas), or no conversion potential (e.g., Australian tropical savannas). However, our analysis also identifies specific Divisions and individual Ecoregions in the WGT where "browning" is occurring in Argentina, Australia, East Africa, Southern Africa, and Eurasia, some of which were also identified by [19,20]. Here, the GVFCP is able to associate the "browning" with explicit trends in both *F*BS and *F*NPV providing more insight into the changes that are occurring.

The analysis at WGT Divisional level identifies "hot spots" of change in VFC based on both area and percentage area of significant positive or negative trends. The total area is important because of the implications for the productivity of the ecosystem, the implications for livelihoods and food supplies, and the impacts on human and wildlife populations. The percentage areas are important because they identify levels of stress on regional and unique representatives of global vegetation Formations with consequences for biodiversity, endangered species, and extinction risks. When WGT Divisions are ranked by percentage of their area with significant positive trends in *F*BS, the first five Divisions are the Patagonian Cool Semi-Desert Scrub and Grassland, the Pampean Grassland and Shrubland, the East African Xeric Scrub and Grassland, the Californian Grassland Meadow, and the Madagascan Montane Grassland and Shrubland with between 38% and 50% of their area affected. Potential for land degradation arising from increasing *F*BS is occurring in large regions such as Patagonia, where aridity and grazing pressure are reducing the cover of palatable grasses [47] and the pampas of South America where conversion from woodland and pasture to cropland has occurred [48]. It can also arise in smaller regions such as the Central Valley Grasslands of California with recent land use change [49], and the very small Madagascan ericoid thicket ecoregion subject to the same pressures on land cover more widely evident across the island [50,51].

However, some of the largest areas of positive trends in *F*BS and negative trends in *F*PV and *F*NPV occur in the largest WGT Divisions such as the North Sahel Semi-Desert Scrub and Grassland which is so large that it crosses many national boundaries and contains a patchwork of areas with both positive and negative trends in *F*PV, *F*NPV, and *F*BS. Analysis found that although there was about 200,000 km2 of positive trends in *F*PV, there were massive areas of almost 800,000 km<sup>2</sup> showing a negative trend in *F*NPV and about 600,000 km2 showing a positive trend in *F*BS. The drivers and manifestations of

environmental change in the Sahel have been scientifically contested for some time [52]. The mixed trends in VFC across the NSSDSG tend to reinforce the contested debate and emphasize that both "greening" and "browning" are occurring at the same time in different locations [53,54]. For example, agricultural practices that reduce vegetation cover, i.e., produce a positive trend in *F*BS, increase potential for wind erosion, dune migration, and soil loss [55]. Hence, negative trends in *F*NPV may indicate removal of both dry grass from rangelands, and residue from croplands. However, relationships with precipitation support areas of greening (positive trends in *F*PV) which may be due to agroforestry [56] and recovery of woody vegetation [54], and promotion of tree cover around Sahelian farms contrasts with reduced woody cover on farmlands in the sub-humid zone such as the WCAMWS [57].

For the SWSG Formations (TLSGS, TMSGS, and WSFDSG), the trend analysis revealed a mosaic of significant positive and negative trends in *F*BS, *F*NPV, and *F*PV covering 42.1%, 26.9%, and 21% of the area, respectively. This represents a massive change in vegetation fractional cover in the grassy savanna, woodland, and scrub biomes of the world over the past 20 years. When the analysis is focused on the individual Ecoregions, trends in VFC can be associated with particular factors unique to specific Ecoregions. In East Africa, in the Acacia-Commiphora Woodlands (EAXSG), positive trends in *F*BS may be associated with increased logging of woodlands for charcoal production [58–60], recent high frequency of droughts [61] due to a decline in long-season rains [62,63], and long-term impacts of expansion of populations, cropland and pastoralism ([64–66]. Land use changes from forest cover to cultivation, pastoralism, and plantations has led to increased surface run-off but variability exists due to site-specific factors [67]. In the Chihuahan Desert ecoregion, invasive grasses and shrubs may be changing the ecosystem dynamics leading to the observed negative trends in *F*BS and positive trends in *F*NPV with consequent increase in fire risk [68–70]. Since 2000, there has been a major expansion of cropping, exotic pastures, and oil palm plantations especially in the western Llanos of Colombia leading to the moderate positive trends in *F*PV [71], with recent potential for acceleration post the Peace Agreement [72,73]. The Mitchell Grasslands of Northern Australia and especially the western Barkly Tableland exhibit a strong positive trend in *F*BS that may relate to long term interactions between cycles of precipitation and grazing pressure [74].

The GVFCP has explicit quantitative uncertainty estimates for the eight-day 500 m product based on the multiple comprehensive cycles of calibration and validation undertaken for Australia [4,5]. In this study we used the aggregated monthly 5 km2 resolution product which has values that represent the medoid of the eight-day values for each month averaged across a hundred 500 m pixels. The global scale for comparison of WGT Formations, Divisions, and Ecoregions and the size of even the small Ecoregions leads to the values presented in the fence box plots being derived from hundreds to hundreds of thousands of pixels with variation within and between Ecoregions well in excess of the published RMSE values of between 11% and 16% for *F*PV, *F*NPV, and *F*BS. If the magnitude of trends, is examined for Ecoregions with trend values beyond the upper and lower fences, values of between ± 0.5% to <sup>±</sup> 1.0% yr−<sup>1</sup> result in changes in average VFC across Divisions and Ecoregions of between 9% and 18% over the 18 years of the time series. Upper fence values for standard deviations of pixel trend values within ecoregions were as high as 1.2% yr−<sup>1</sup> and as low as 0.1% yr−<sup>1</sup> indicating that there was major heterogeneity in responses within Ecoregions, related to their size, population distributions, land use and land cover, and between Ecoregions where vegetation Formation and geographical location, as well as anthropogenic factors drive variation.

There is extensive scope for much more detailed examination of the trends in the GVFCP across the globe, and especially within the WGT, and at finer scale for smaller grassy vegetation types. Specifically, there is a need to segment the time period from 2001–2018 and explore discontinuities and changes in trends that may be attributable to changes in climate cycles, land use, and production systems. In addition, there is a need to explore with spatial explicitness, the association between changes in *F*PV, *F*NPV, and *F*BS; and examine the trends in total cover (*F*TC) as an indicator of wind and water erosion potential [13]. There is also a need for the product to be widely used beyond current Australian applications and be subject to further field validation analysis especially in particular

ecosystems with very different vegetation structures, arboreal vegetation phenology, and soil surface colours and conditions. Since the GVFCP is freely available on the GEOGLAM RAPP-MAP website, the authors are encouraging scientists and agencies to explore and utilize these data.

#### **5. Conclusions**

This study has presented the GVFCP and explored the levels and dynamics of *F*PV, *F*NPV, and *F*BS across the grassy ecoregions of the world. It has documented benchmark long-term average levels, and areas and magnitudes of significant trends in VFC at Divisional level for WGT Ecoregions. Focused exploration of the savanna woodlands showed that:


The GVFCP is made freely available to the global scientific community in the hope that they will use it and provide further validation studies. There is great potential for development of a Landsat/Sentinel 2 global product using the same methodology. This will become necessary upon the demise of the MODIS sensors since the Visible Infrared Imaging Radiometer Suite sensor does not carry the 2105–2155 nm short wave infrared channel and sensitivity to *F*NPV and *F*BS will likely be diminished.

**Online Resources:** Global and Australian products are freely available along with analytical tools at https: //map.geo-rapp.org.

**Author Contributions:** Conceptualization, J.P.G. and M.J.H.; Methodology, J.P.G.; Validation, J.P.G. and M.J.H.; Formal analysis, M.J.H. and J.P.G.; Investigation, M.J.H.; Resources, J.P.G.; Data curation, J.P.G.; Writing—original draft preparation, M.J.H.; Writing—review and editing, M.J.H.; Visualization, J.P.G.; Project administration, J.P.G.; Funding acquisition, J.P.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research spans 15 years. The initial stage was partially funded by the NASA Earth Science Enterprise Carbon Cycle Science research program (NRA04-OES-010). Development of the current Australian and Global products was funded by the Australian Government's National Landcare Program and CSIRO. GEOGLAM RAPP is an initiative of CSIRO, GEO, and GEOGLAM.

**Acknowledgments:** The GVFCP is maintained for GEOGLAM RAPP by Biswajit Bala.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

**Table A1.** Acronyms and abbreviations used (in addition to WGT acronyms documented in Tables A2 and A3).



**Table A1.** *Cont*.




**Table A2.** *Cont*.

**Table A3.** Savanna, woodland and shrubland ecoregions with World Grassland Types (Dixon et al., 2014). (See Table A2 for description of Formations and Divisions; Continent codes as follows: AF—Africa; NA—North America; SA—South America; AU—Australasia (includes New Guinea and New Zealand)).



**Table A3.** *Cont*.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **Terrestrial Laser Scanning for Vegetation Analyses with a Special Focus on Savannas**

**Tasiyiwa Priscilla Muumbe 1,\*, Jussi Baade 2, Jenia Singh 3, Christiane Schmullius <sup>1</sup> and Christian Thau 1,4**


**Abstract:** Savannas are heterogeneous ecosystems, composed of varied spatial combinations and proportions of woody and herbaceous vegetation. Most field-based inventory and remote sensing methods fail to account for the lower stratum vegetation (i.e., shrubs and grasses), and are thus underrepresenting the carbon storage potential of savanna ecosystems. For detailed analyses at the local scale, Terrestrial Laser Scanning (TLS) has proven to be a promising remote sensing technology over the past decade. Accordingly, several review articles already exist on the use of TLS for characterizing 3D vegetation structure. However, a gap exists on the spatial concentrations of TLS studies according to biome for accurate vegetation structure estimation. A comprehensive review was conducted through a meta-analysis of 113 relevant research articles using 18 attributes. The review covered a range of aspects, including the global distribution of TLS studies, parameters retrieved from TLS point clouds and retrieval methods. The review also examined the relationship between the TLS retrieval method and the overall accuracy in parameter extraction. To date, TLS has mainly been used to characterize vegetation in temperate, boreal/taiga and tropical forests, with only little emphasis on savannas. TLS studies in the savanna focused on the extraction of very few vegetation parameters (e.g., DBH and height) and did not consider the shrub contribution to the overall Above Ground Biomass (AGB). Future work should therefore focus on developing new and adjusting existing algorithms for vegetation parameter extraction in the savanna biome, improving predictive AGB models through 3D reconstructions of savanna trees and shrubs as well as quantifying AGB change through the application of multi-temporal TLS. The integration of data from various sources and platforms e.g., TLS with airborne LiDAR is recommended for improved vegetation parameter extraction (including AGB) at larger spatial scales. The review highlights the huge potential of TLS for accurate savanna vegetation extraction by discussing TLS opportunities, challenges and potential future research in the savanna biome.

**Keywords:** terrestrial laser scanning (TLS); savanna; Above Ground Biomass (AGB); 3D point cloud; vegetation structure

#### **1. Introduction**

#### *1.1. The Role and Importance of Savannas*

Savannas are heterogeneous ecosystems characterized by mixed physiognomies of woody and herbaceous layer [1]. Savanna ecosystems span one-fifth of the global land [2,3] which makes up approximately 20% of the terrestrial vegetation and contribute 30% of all terrestrial ecosystem gross primary productivity [4,5]. Savanna ecosystems also contribute significantly to the global carbon cycle with a net primary productivity

**Citation:** Muumbe, T.P.; Baade, J.; Singh, J.; Schmullius, C.; Thau, C. Terrestrial Laser Scanning for Vegetation Analyses with a Special Focus on Savannas. *Remote Sens.* **2021**, *13*, 507. https://doi.org/10.3390/ rs13030507

Academic Editor: Michael J. Hill Received: 28 December 2020 Accepted: 28 January 2021 Published: 31 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of 1 to 12 t C ha−<sup>1</sup> year−<sup>1</sup> and a carbon sequestration rate averaging 0.14 t C ha−<sup>1</sup> to 0.39 Gt C year−<sup>1</sup> [4].

Savannas are vital biomes that play a major role in the provision of ecosystem services. Moreover, they contribute significantly to the socio-economic needs of most rural households, especially in Africa, where many are dependent on natural resources for their livelihood [6,7]. Savannas support approximately 20% of the world population and account on average 28% of household income in rural areas of developing countries [8,9]. In west and central Africa, income from bushmeat is equivalent to approximately \$1000 per year [7].

Since savannas are home to 20% of the world's population, humans have caused changes in these ecosystems through land clearing for agriculture and livestock production, infrastructure, fire lighting or suppression and fuelwood collection [9,10]. The alteration of savanna structure and function has direct effects on woody vegetation and results in woody encroachment as a result of the interaction of drivers mainly herbivory, fire and atmospheric carbon dioxide [11,12]. Anthropogenic land use and land cover change affects vast areas of the savanna ecosystem and the consequences can be substantial for the regional climate cycle [9]. Against this background, it is important to regularly map and monitor savannas and to quantify the carbon pools and dynamics of these ecosystems because they contribute significantly to global productivity [13].

#### *1.2. Monitoring of Savannas Using Remote Sensing*

Much of the current understanding on savanna vegetation structure is derived from traditional field-based methods, such as forest inventories, permanent sample plots and transects [14,15]. Field-based estimates are essential in the validation of remote sensed data [16], though they are limited in that they are time consuming, expensive, labour intensive and often require destructive sampling [17,18]. Remote sensing techniques are more appropriate for monitoring when compared to traditional methods [19]. Especially for consistent analyses at different spatial and temporal scales and when considering inaccessible areas, remote sensing with limited fieldwork has been proposed as an efficient method [17,19].

The main types of Earth observation data used for vegetation studies are based on optical (multispectral, hyperspectral), Radio Detection and Ranging (RADAR) and Light Detection and Ranging (LiDAR) sensors [20].

Multispectral data can be acquired from spaceborne (satellites) or airborne (aircrafts, Unmanned Aerial Vehicles (UAVs) platforms and from the ground imagery (spectrometers). Multispectral data has the advantage of higher spatial resolution [21] and is publicly available for most regions in the world e.g., Landsat data. It is limited in terms of spectral resolution [21] and spaceborne optical data often suffer from cloud cover and mixed pixels [20,22].

Hyperspectral data can be acquired from spaceborne and airborne sensors and hyperspectral cameras [23]. Hyperspectral data is acquired in hundreds of narrow bands, and thus it has a higher spectral resolution enabling distinguishing of vegetation types [24]. It is limited though in terms of spatial resolution [21] and suffers from band redundancy where bands are strongly correlated [24].

RADAR data have the advantage of being largely independent of atmospheric and illumination conditions. However, they are limited because data-specific issues need to be compensated for or at least considered when using RADAR for vegetation monitoring [25]. Among others, these issues comprise the geometry of Synthetic Aperture RADAR (SAR) imagery, the dielectric properties of the surveyed surfaces, RADAR speckle complications due to the separation between signal and noise in the data, RADAR shadow and saturation problems [25].

To date, LiDAR systems are considered to be the most accurate and recommended remote sensing technology for vegetation parameter retrieval by producing high-resolution 3D data called point clouds [26,27]. They are available as space and airborne, mobile and UAV as well as terrestrial LiDAR sensors [27]. Spaceborne LiDAR provides the opportunity to make large-scale assessments of vegetation structure for global to regional biomass estimation [28]. Airborne LiDAR is also able to cover relatively large areas and measures vegetation structure with high accuracy [29]. Both spaceborne and airborne LiDAR systems are limited in detecting smaller trees and shrubs. The latter are key attributes in savanna ecosystems [30]. By using airborne or spaceborne LiDAR these attributes and are restricted to landscape level estimates. Mobile and UAV LiDAR sensors offer a faster approach for acquiring 3D point clouds. However, mobile LiDARs are constrained in extracting canopy and height attributes [31], whereas UAV LiDARs have difficulties in capturing tree trunks [32] and are also limited in flight time due to batteries [33]. An overview of the main advantages and disadvantages of the described monitoring technologies i.e., (field surveys, multi-spectral, hyperspectral, RADAR and LiDAR) platforms are also discussed in Béland et al. [27], Eisfelder et al. [34] and Galidaki et al. [35].

#### *1.3. Terrestrial Laser Scanning (TLS)*

Terrestrial Laser Scanning (TLS) has grown in interest because of its ability to produce 3D point cloud data for characterizing vegetation structure [36]. Accordingly, over the past decade, a considerable amount of research has been directed towards the 3D quantification of vegetation structure using TLS [37,38]. Vegetation structural parameters are difficult and laborious to measure in the field, especially in complex ecosystems like the savanna. However, TLS has proven efficient in deriving these parameters because it produces consistent and objective measurements with high precision. Furthermore, TLS measurements are relatively easy to interpret and parameter extraction can be automated [39,40].

A terrestrial laser scanner or ground-based LiDAR is a LiDAR system from which detailed 3D information called point clouds can be collected [41]. Using a laser beam, terrestrial laser scanners provide a high-resolution 3D view of the objects [36]. Depending on their range measurement principle, terrestrial laser scanners are classified into two types: phase shift or time of flight scanners [36,42]. Phase shift scanners measure distances by analysing the shift between the emitted and the received laser wave. In contrast, time of flight scanners compute distances based on the difference between the emission and reception of the laser pulse [42].

TLS measurements can be collected based on single or multiple scans. Using single scans, the scanner is placed at a single location within the study area whilst multiple scans are obtained when the scanner is placed at several positions within the test site [43]. For rapid and efficient extraction of vegetation parameters with reduced fieldwork and faster post-processing, single scans are preferred [44]. However, this approach is prone to occlusions (e.g., a tree being in the scan shadow of another tree or any other object) reducing the accuracy of vegetation structure retrieval and the ability to get a full 3D representation of an object [44]. For the best plot coverage and a complete 3D description of the objects under consideration, a multi-scan TLS setup is necessary [31,36].

Multi-scan TLS data allow for reliably retrieving common tree attributes such as Diameter at Breast Height (DBH), tree height, canopy structure and vertical height profiles [45,46]. In addition to these parameters, more detailed and complex vegetation attributes such as leaf water content [47], basal area and stem density [48], Above Ground Biomass (AGB) change [49], tree crown change [50] stem curve profiles [51], stem shrinkage [52], have been derived from TLS point clouds.

TLS provides unmatched 3D information on vegetation structure right down to branch and leaf scales, making it far superior to data obtained through traditional field work [38,53]. Using the intensity of the returned laser pulses, TLS can be used to map pigment concentration, photosynthetic capacity and water content of the vegetation [27]. It is also possible to separate wood from leaf material [54,55]. TLS offers a non-destructive approach to quantify canopy and stem volume with high accuracies, thereby enabling the estimation of AGB with reduced uncertainties as well as the development and improvement of reliable allometric equations [56,57]. TLS also enables the estimation of individual

tree components [58] and it can also be employed as a calibration and validation dataset because of its ability to acquire data from the ground that is similar to traditional measures of vegetation structure [59,60].

The inherent heterogeneity of savannas presents a challenge regarding accurate vegetation characterization. To overcome this challenge, 3D data on the horizontal and vertical structure of the savanna are needed. TLS is able to account for savanna heterogeneity (Figure 1) by providing detailed information on all savanna vegetation types and layers [61,62]. It therefore holds much potential for making holistic AGB assessments, including long and short stature vegetation [63–65].

Different drivers can result in the loss of large trees in the savanna landscape [66] and changes in vegetation structure. The non-destructive nature of TLS makes it more so suitable for quantification of change in vegetation structure and reduction of uncertainties when estimating AGB [50,67]. The use of repeated TLS scans allows for monitoring these and other losses (and potentially gains) in savanna AGB over time. Multi-temporal TLS therefore holds the potential to support carbon accounting programmes such as REDD+.

While TLS is able to capture the 3D configuration of vegetation at high spatial resolution, data can only be acquired for a limited areal extent and data acquisition usually requires extensive field work [68]. This problem is solved by the integration of TLS with other remote sensing data to ensure the coverage at landscape scales. For example, landscape studies were conducted in the savanna by the integration of TLS with airborne LiDAR [69], TLS with MODIS satellite data [70] and TLS with ALOS PALSAR L band data [62]. Moreover, TLS scans frequently feature occlusions ("shadows") caused by objects close to the scanner which, in turn, block the view of other objects further away. Although occlusions are not a major problem in open savannas, they can be minimized by the use of multiple TLS scans [30,64]. Another shortcoming of TLS systems is that the point density decreases with distance from the scanner, an effect called beam divergence [68,71] which causes the point cloud to be unreliable at a certain distance.

#### *1.4. Overall Goal and Objectives of the Review*

This paper aims to review the scientific literature on the use of TLS for vegetation analyses. It puts a special but not exclusive geographic focus on savannas in order to establish the current status quo in savanna vegetation parameter extraction. In contrast to existing reviews in the TLS literature [72,73], this study presents the current TLS methods in practice and identifies those suitable for accurate biomass assessments in the savanna biome. To address this gap, our meta-analysis takes a detailed look at the geographic distribution of TLS studies, the parameters retrieved from TLS point clouds, known retrieval methods as well as the main opportunities and challenges arising from the studied literature. The remainder of the paper is structured as follows: Section 2 presents the approach and Section 3 the main findings of the review. The main discussion points and conclusions derived from the review are presented in Sections 4 and 5, respectively.

#### **2. Materials and Methods**

A systematic literature review was conducted in Web of Science [74]. Access to its catalogue was granted through the Thüringer Universitäts- und Landesbibliothek Jena (ThULB). The advanced search functionality in Web of Science was used to develop and apply a set of predefined queries. The search included all journals in the Web of Science Core Collection™ database. Peer-reviewed conference papers were also included. Searching for keywords in the title (TI) and in some instances in all fields (ALL) (when a query returned fewer than ten articles) was performed with a time frame up to May 2020. The time span for the search was all years that is from 1945 to 2020. In the event that a search query returned more than 100 articles, the search was refined by the research areas forestry or remote sensing to narrow down on specific studies employing TLS for vegetation parameter extraction. Each query consisted of different keywords and related synonyms (Table 1) which were combined through logical operators (i.e., "AND", "OR").


**Table 1.** Search terms used to build the search queries for this review.

An example of one of the search queries developed is: TI = ((TLS OR Terrestrial LiDAR OR Terrestrial Laser Scann \*) AND (Tree Structur \*)). This search query would return all journal articles that have either the word or terrestrial LiDAR or terrestrial laser scanner or terrestrial laser scanning and the word tree structure or tree structural in their title.

Following the above search scheme, a total of 158 research articles were initially retrieved from Web of Science. In a second step, certain papers were retained for the review according to a prescribed set of criteria. The criteria used to refine the initial search results were as follows:


This resulted in a final number of 113 reviewed publications. To qualitatively and quantitatively assess the final selection of papers, a comprehensive meta-analysis was conducted based on 18 predefined study attributes. Besides author(s), year of publication and publishing journal, the list of reviewed meta-parameters included: study area(s) (country, continent, biome); type of laser scanner used (TLS model/device); range finder of the scanner (phase shift or time of flight); type of TLS data used (xyz and intensity); scanner setup (TLS settings for scanning); field setup (use of reflectors, GPS/dGNSS); number of scans (single vs. multiple scans); parameters retrieved (e.g., DBH, vegetation height, crown diameter); retrieval method (algorithm used); number of studied trees and shrubs (total number of trees/shrubs measured in the research); performance metrics such as coefficient of determination (R2), Root Mean Square Error (RMSE) and bias; main discussion points (e.g., advantages and limitations of a used approach); and lastly research gaps and potential avenues for future research. The above search efforts were undertaken to ensure that all relevant TLS literature was collected. The full list of papers considered and finally selected in this review is available in Table S1 (supplementary material).

#### **3. Results**

#### *3.1. General Overview*

Of the 113 publications used for the detailed review, 80% of the publications focused on the extraction of vegetation parameters from TLS point clouds, 12% used TLS as auxiliary data (e.g., TLS data used to validate airborne LiDAR) [78–80] and 8% applied TLS for technical experiments (e.g., reflectance modelling using TLS intensity data) [81–83]. Among the 113 papers, 97 were publications from peer-reviewed journals whilst 16 were peer-reviewed proceedings. The most popular journals were Remote Sensing (20 studies), Remote Sensing of Environment (9 studies), Agricultural and Forest Meteorology (8 studies) and ISPRS Journal of Photogrammetry and Remote Sensing (6 studies). Most proceedings papers stem from the International Geoscience and Remote Sensing Symposium (IGARSS), with 9 studies.

#### *3.2. Temporal and Spatial Patterns*

Figure 2 shows the trend observed in the number of TLS studies. Beginning in the early 2000s, TLS technology evolved slowly. New systems with a range of capabilities and applications emerged and this explains the rise in the number of studies from the early 2000s to date. TLS technology had a real breakthrough in vegetation studies around 2013 and the upward trend seems to continue in 2019. A further significant increase in peer-reviewed studies can be observed for 2018 and 2019 (25 and 22 articles, respectively). During this time, the number of published papers per year has more than doubled when compared to the period after the breakthrough (on average 10.8 studies per year for 2013–2017).

The geographic distribution of the reviewed studies by country is shown in Figure 3. To date, most studies were conducted in the United States, China and Finland. From a continental perspective, Europe has the highest number of studies (38%) [84–86] followed by North America (24%) [46,59,87] and Asia (19%) [88–90]. Africa (6%), Australia (5%) and South America (7%), with more than half of their landmass belonging to the savanna biome [1,2], had few TLS studies [64,91,92].

The geographic distribution of studies by biome (after Olson et al. [93]) and the percentage of reviewed studies per biome are shown in Figures 3 and 4. Studies that were conducted in more than one biome were counted multiple times and four laboratory studies were excluded from this analysis, resulting in a sample size of n = 114. The reviewed papers covered a total of 11 biomes. While most studies were conducted in temperate forests (47%) [75,94,95], much less papers were published in boreal (16%) [96–98], tropical and subtropical forests (10%) [99–101] as well as grasslands and savannas (9%) [70,91,92]. The cumulative total for the savanna biome was 13%, including 4% of savanna research conducted in temperate regions ("Temperate Grasslands Savannas and Shrubland").

**Figure 2.** Number of published studies per year and continent. Over the past 15 years, research on TLS for vegetation assessments has increased significantly and knowledge has been generated especially for the European, North American and Asian region. Please note that the low number of studies in 2020 is an artefact of the review period ending on 31 May 2020.

**Figure 3.** Geographic distribution of reviewed studies. Biome classification according to Olson et al. [93].

**Figure 4.** Percentage of reviewed studies per biome. Biome classification according to Olson et al. [93].

#### *3.3. TLS Instruments and Data Used*

A variety of TLS instruments can be purchased nowadays. They are offered by different manufacturers, with the main companies being RIEGL Laser Measurement Systems GmbH [102], FARO Technologies Inc [103] and Leica Geosystems [104]. Among others, TLS instruments differ in their, measurement range (the maximum distance that can be scanned e.g., 600 m), sample rate (pulses per second e.g., 122,000), wavelength (e.g. 1550 nm (near-infrared)), angular sampling (angle measurement resolution e.g., 0.02◦), waveform (discrete or full-waveform), range measurement principle (time of flight or phase shift), size, mass and power efficiency. Based on our review, Table 2 provides an overview of frequently used TLS instrument models. Among the most popular scanners were Riegl VZ 400 (26 studies) [90,105,106], Leica HDS6100 (10 studies) [107–109], Riegl VZ 1000 (9 studies) [110–112], Faro Focus 3D (7 studies) [50,67,113] and Faro Focus 3D X 330 (also 7 studies) [83,114,115].

Figure 5 shows the number of scan positions and the area covered in the reviewed studies. It becomes clear that a larger area requires a higher number of scan positions, the higher the number of scan positions the fewer the studies and most studies use up to 10 scan positions and cover an area of less 5000 m2. The number of scan positions and area covered also depend on the vegetation type and the level of detail required for the specific research. For high detail data, some studies employed a high number of scan positions on a relatively small area. Li et al. [63] used 18 scan positions on a 400 m<sup>2</sup> plot on shrubs and Lau et al. [99] used 16 scan positions on each selected tree in a 1200 m<sup>2</sup> plot in a tropical forest. Most inventories follow a field setup with multiple scans, with five

(28 studies) [121–123], four (23 studies) [124–126] and three (13 studies) [113,127,128] scan positions being the most popular. However, a single scan experimental design is also very common (28 studies [83,89,129]).


**Table 2.** Commonly used terrestrial laser scanners.

**Figure 5.** Number of scan positions and area covered (m2) in the reviewed studies.

TLS instruments provide two types of data that can be used for forestry and vegetation studies. These are 3D coordinates (XYZ data, point clouds) and the intensity of the reflected signal. The signal can either be discrete or full-waveform. Discrete data provide single or multiple returns for each laser pulse [76,130]. Full -waveform contain measurement of the reflected laser beam as a function of range [76,130]. These two data types can be integrated with auxiliary data from the field, including RGB images from digital cameras. Most studies used XYZ data (85 studies) [131–133] whereas intensity information was considered much less (8 studies) [129,134,135]. A triple combination of both data types and RGB images (2 studies) [121,136] as well as a combination of XYZ and RGB (3 studies) [75,111,137] was not commonly observed. Instead, a number of studies took advantage of XYZ and intensity data (15 studies) [127,138,139]. No paper presented results based on a combination of TLS intensity data and RGB images.

#### *3.4. TLS Point Cloud Pre-Processing*

The scanner in use usually comes with software for pre-processing the point cloud; for example, Riegl scanners come with Riscan Pro, Faro with FARO Scene and Leica with Leica Cyclone. The first stage in prepossessing if a multi scan setup was used is to conduct coarse registration and georeferencing to the local coordinate system using the retroreflectors as tie points [43,96]. Fine registration can also be achieved using multi-station adjustment through an iterative closest point algorithm in Riscan Pro software [43,140]. Multi-station adjustment requires a minimum of 3 TLS setups with a minimum overlap of 30%. The algorithm uses surface and point normals to compute the adjustment [141]. The registered point cloud is then clipped to the radius of the plot in the case of circular plots or to the extent of the plots or to remove low density points at the plot edges [140,142].

The final pre-processing stage involves filtering the point cloud to remove noise from the data [43,121,142]. Noise and outlier removal is done to increase accuracy and reduce errors when extracting information from the point cloud data [115]. The noise can be classified into internal and external noise, internal noise being the inherent noise from the scanning device and external noise from the reflective nature of the scene [143,144]. The common noise points from vegetation are leaves, extremely high or low points, and from occlusions [115,144]. The filtering is performed using various approaches such as statistical, neighbourhood, projection, signal and partial differential based filtering techniques [144]. For example, in the case of statistical outlier removal, the filtering reduces the error on the point cloud by computing the average distance to its neighbours and rejects points that are further away from the average distance [115,144]. Point cloud filtering can be performed in software packages such as Cloud Compare, LAStools, MeshLab and programming libraries such as Point Cloud and Open Graphics Library [143]. After pre-processing various algorithms are then applied to extract the vegetation parameters.

#### *3.5. Methods Used with TLS Point Clouds*

Of the 113 papers reviewed, 87 papers focused on the extraction of vegetation parameters from TLS point clouds. Of the 87, only 15 studies worked in the savanna biome [39,145,146]. The most frequently employed methods to derive vegetation attributes are the RANdom SAmple Consensus (RANSAC) algorithm for extracting DBH, stem curve profiles and the detection of stems [97,147], the use of the highest z coordinate of the point cloud for estimating heights [61,148], voxel-based and radiative transfer models for assessing LAI [149,150], Canopy Height Models (CHMs) for delineating crown attributes [91,92], and Quantitative Structure Models (QSMs) for computing tree volume and branch parameters [151,152]. Popular software packages to implement and work with these methods are Lastools [153], Matlab [154], R [155] and Python Programming [156], Cyclone [157], FARO Scene [158], RiscanPro [159], CompuTree [160] and Cloud Compare [161]. For evaluating the above extraction techniques, performance metrics such as R2, RMSE and bias are typically used [51].

#### *3.6. Vegetation Parameters Extracted from TLS Data*

The vegetation parameters that can be retrieved from TLS data can be broadly classified into primary and secondary parameters [22]. Primary parameters are those directly derived from the TLS point cloud, such as DBH, height, volume and tree canopy characteristics [22]. Secondary parameters can only be estimated indirectly and are derivatives of primary parameters. They are usually obtained through modelling [22]. Examples of secondary parameters are AGB, basal area, stem density and LAI [22]. Besides the common vegetation parameters e.g., DBH and height, some studies focused on extracting unique and complex vegetation parameters such as tree planar projection [162], the elevation of angle and azimuth of stem inclination [52], clumping index [163], ground cover [92] and structural loss [164].As shown in Figure 6, most of the reviewed studies retrieved primary parameters, including DBH (41), height (32), tree canopy features (e.g., crown diameter, canopy cover) (20), trees or tree stems (20) and tree volume (15). With only five and six articles, respectively, branch parameters and stem curve profiles gained the least attention. Of the secondary parameters, LAI was the most extracted parameter (21), followed by AGB (11), basal area (5) and stem density (3). Details regarding specific parameters and the methods to retrieve them are provided in the following subsections. Particular emphasis is put on savanna ecosystems.

**Figure 6.** Primary and secondary vegetation parameters retrieved in the reviewed TLS studies. Studies were counted multiple times if two or more parameters were retrieved.

#### *3.7. Primary Vegetation Attributes* 3.7.1. DBH

From the reviewed papers, only one study derived DBH in the savanna [92]. Muir et al. [92] employed circle fitting to single-scan TLS data at heights between 1.28 m and 1.32 m above the Digital Elevation Model (DEM) in an open woodland site in Central Western Queensland Australia. They compared the resulting stem diameter measurements with data collected in the field and obtained an R2 value of 0.77 [92]. Slightly better performances are often observed for DBH extraction in temperate, tropical and boreal forests, with R2 values ranging from 0.81 to 0.99 [94,165,166] when compared to field reference data. For a pine woodland test site in France, least-squares circle fitting was reported to be the most accurate method (RMSE = 0.019 m) for estimating DBH when compared to least squares cylinder fitting and circular Hough transformation [85]. Similarly, Calders et al. [40] and Chen et al. [167] obtained R2 values of 0.97 when using least squares circle fitting for DBH estimation with a multi-scan TLS scan with 5 positions in a eucalypt open forest and boreal forest respectively. Yurtseven et al. [148] obtained an R<sup>2</sup> value of 0.99 with a multi-scan with 8 positions in a Mediterranean forest employing randomised hough transformation with least square regression. In general, multiple TLS scans and a lower stem density yielded the best results [114,168,169].

#### 3.7.2. Vegetation Height

Tree height estimation in the savanna was reported in five studies [30,62,91,92,145]. Odipo et al. [62] and Muir et al. [92] employed the use of CHMs to estimate tree height from multiple and single TLS scans (30 positions and single scan) and obtained R2 values of 0.64 and 0.98 respectively when compared to ground referenced data. Singh et al. [30] observed an RMSE of 0.25 to 1.45 m based on multi-scan Long Range -TLS (LR-TLS) with 4 positions. They computed heights by measuring the vertical distance between treetops and the stem base. Kirton et al. [145] employed cylinder fitting to extract tree height and obtained a standard error of ±0.002 m using a TLS single scan. It was observed in this study that it is more challenging to derive the height of savanna shrubs than it is for savanna trees. Muir et al. [92] also reported that shrub heights are difficult to calculate from TLS data, with R<sup>2</sup> values of 0.381 and 0.08 for maximum and average shrub height, respectively. In contrast to this observation, the extraction of height from short stature vegetation was achieved with reasonable accuracy in montane areas as well as desert and xeric shrubland biomes [170]. Vrlingie et al. [170] applied a 2-D spatial wavelet analysis on a TLS derived CHM to extract shrub height and achieved an R2 of 0.94 when compared to field measured shrub height.

#### 3.7.3. Stem Detection

The detection of stems from TLS point clouds is an important first step in processing, to enable extraction of other vegetation parameters [171]. Only one study [147] from the review was conducted in the savanna biome. Burt et al. [147] employed treeseg algorithm (which extracts tree point clouds from larger point clouds) on a tropical and on an open vegetation cloud using TLS multi scans with 5 positions. They achieved 70% and 96% stem detection accuracy, respectively [147]. The detection of stems was achieved with high accuracies in temperate and boreal forests [138,171]. Kelbe et al. [138] used voxel stem segment modelling on a single TLS scan dataset acquired using low resolution TLS in various forest types in New England. Stems were detected with accuracies of (R2 = 0.99; RMSE = 0.16 m) [138]. In a boreal forest Zhang et al. [171] employed connected component segmentation using both single and multi- scan with 5 positions. They achieved stem detection accuracies of 93% for single scan and 99% for multi-scan TLS setup when compared to reference data [171]. To date, most of the commonly used algorithms (e.g., density-based clustering [172], Hough transformation [85] and RANSAC cylinder fitting [48,173] used for stem detection achieved stem detection accuracies of 70% to 99% when applied in temperate, tropical and boreal forests.

#### 3.7.4. Crown Attributes

Among the tree crown attributes dealt with in the scientific literature are canopy cover, canopy height/length, crown width/diameter and tree crown change. In the savanna biome, extraction accuracies for trees were better as compared to shrubs. Singh et al. [30], Odipo et al. [62] and Muir et al. [92] employed CHMs to determine canopy cover. Singh et al. [30] achieved an RMSE of 5.7–15.9% at distances up to 500 m with 4 scan positions, whilst Odipo et al. [62] and Muir et al. [92] achieved R2 of 0.56 and 0.93 as compared to the ground referenced data with 30 and single scan positions respectively. Lower accuracies were obtained when estimating shrub canopy cover, with Muir et al. [92] reporting an R<sup>2</sup> of 0.23 when compared to field measurements. Yet, the application of 2D spatial wavelet analysis achieved an R<sup>2</sup> of 0.51 for 36 shrubs in a desert and xeric shrubland environment [170]. In this study, shrub crown area was calculated using spatial wavelet analysis using the detected radius while the actual shrub crown area was determined by measuring the minor and major axis and then calculating the crown as an ellipse [170]. The use of different algorithms for accurate extraction of canopy cover was also investigated by Yurtseven et al. [148] in a mediterranean woodland and scrub. They employed convex and concave hulls and they concluded that the concave planar underestimates canopy cover area as compared to the convex hull. A method to quantify canopy cover change was considered by Olivier et al. [50] in a sugar maple and balsam fir forest in Canada. They used bi-temporal multi-scan TLS data to identify vegetation boundaries, extract new material formed or displaced between the two-time period [56].

#### 3.7.5. Tree Volume and Branch Parameters

Volume of standing trees or above ground volume is total stem wood and branches [174]. Branch parameters are parameters like branch total length, branch diameter, branch insertion angle, branch depth into the crown and branch radius. Over the years, a range of algorithms has been developed to compute tree volume and branch parameters from TLS point clouds [151,175]. However, none of the corresponding studies was conducted in the savanna biome. QSMs are the most common algorithms used in literature [96,151,175]. They have the ability to provide reliable volumes estimates of stem and branches of trees, but also have limitations in reconstructing leaves [176]. Figure 7 shows an example of a 3D point cloud and the resulting QSM. Branching structure and branch lengths were determined in boreal and tropical forests using QSM [151,175]. Kaasalainen et al. [175] used a multi scan with 2 positions to reconstruct branching structure and concluded that branching structure can be reproduced with ±10% accuracy when compared to laboratory reference measurements. Lau et al. [151] used a multi-scan with 13 positions and successfully reconstructed 87% of the branch lengths. They compared the reconstructed branch lengths to manually measured branches from destructively harvested trees [151]. For the estimation of shrub total and green biomass, Olsoy et.al [54] employed voxel-based and 3D convex hull algorithms to first determine shrub volume in a xeric sagebrush ecosystem in Southern Idaho, USA using multi scan with 2 positions. They then modelled shrub total and green biomass by subsetting the point cloud into green and non-green woody components based on the difference in their reflectivity. Their findings showed that convex hull algorithms (R<sup>2</sup> = 0.92, R<sup>2</sup> = 0.83) performed slightly better than the voxel-based method (R2 = 0.86, R<sup>2</sup> = 0.73) in estimating total and green biomass respectively [54]. The TLS derived shrub biomass was compared to the destructive samples of shrubs [54]. These and other volume-related algorithms (e.g., bounding cuboids, outer hull models, square-based columns) have also been compared for test sites in the temperate and tropical regions [101]. It was found that convex hull techniques tend to overestimate while the square-based column method underestimates volume. Moreover, voxel-based estimates are mainly affected by the point density and level of occlusions in the point cloud [101].

**Figure 7.** (**a**) TLS point cloud coloured by height (lowest values in blue and highest values in red). (**b**) Resulting QSM model from the algorithm by Raumonen et al. [177]. Data Source: (Madhibha, 2016) [178].

#### 3.7.6. Stem Curve Profiles

Stem profiles or stem taper is the "stem diameter measured as a function of height" [179]; that is, how the diameter of the tree changes at different heights along the stem. Six studies used TLS point clouds to determine stem curve profiles, five of the studies were conducted in boreal forests [48,51,52,162,166] and one was conducted in temperate broadleaf and mixed forest [58]. The methods used to extract stem curve profiles from TLS point clouds are either circle or cylinder fitting at different height bins [48,51,52,58,162,166]. Maas et al. [58] extracted stem curve profiles using circle fitting at multiple heights and obtained an RMSE of 4.7 cm when compared to a harvested stem. Liang et al. [51] observed that most algorithms used for the extraction of stem curve profiles results in an RMSE of 1.3 and 6.0 cm for single scan data and 0.9 and 5.0 cm from multi scan data for three different stem densities (i.e. low, medium and high stem densities) in a boreal forest. The RMSE values obtained are a function of both the stem density and the number of scan positions employed [51].

#### *3.8. Secondary Vegetation Attributes*

#### 3.8.1. AGB

AGB is defined as "the total amount of biological material (oven-dried) present above the surface in an area" [180]. AGB is the most important secondary attribute derived from TLS point clouds. The accurate extraction of all other parameters is a necessity in improving tree biomass estimation. A total of five studies estimated AGB in the savanna using TLS [62,64,70,90,126]. Zimbres et al. [64] predicted plot AGB of three vegetation types in the Brazilian cerrado. The TLS derived AGB was validated against the AGB

values obtained from allometric equations. The R<sup>2</sup> values obtained were 0.92 and 0.88 for the woodland savanna dry and rainy season respectively and 0.58 for the forested savanna. Due to occlusions, no conclusive results were obtained for the gallery forest [70]. Cuni-Sanchez et al. [70] estimated AGB change in five different vegetation types using vertical plant profiles extracted from a multi-scan TLS over a period of 20 years in Gabon. However, the savanna plots were not assessed because all the trees had a diameter of less than 10 cm [70]. Odipo et al. [62] also estimated AGB change in a south african savanna using TLS derived canopy cover and height metrics for plot level AGB and then further extrapolating to a wider scale using multi-temporal L band SAR. The TLS derived metrics were validated with field measured parameters. Significant relationships were obtained between the TLS derived predictors (canopy cover & height) with field biomass (R2 = 0.93). Seven variables extracted from TLS data were used to build regression models with field measured AGB to estimate grass biomass in China [90]. The regression methods used were Simple Regression (SR, Stepwise Multiple Regression (SMR), Random Forest (RF) and Artificial Neural Network (ANN). Results obtained showed that the highest prediction of biomass was obtained using SMR model (R2 = 0.84), followed by SR (R<sup>2</sup> = 0.80), RF (R<sup>2</sup> = 0.78) and lastly ANN (R<sup>2</sup> = 0.73) [90]. Cooper et al. [126] investigated the application of Structure from Motion (SfM) photogrammetry and TLS for the estimation of grass biomass and compared this to destructive harvest grass biomass. The R<sup>2</sup> values for the TLS derived grass biomass was 0.46 and total biomass was 0.57, whilst for the SfM it was 0.54 for the grass biomass and 0.72 for the total biomass [126]. TLS data have also been used to derive tree and grass AGB change [49,115]. Srinivasan et al. [49] extracted various tree parameters (including AGB) from loblolly pines in Texas, USA by the use of multi-temporal TLS data. The best results (R<sup>2</sup> = 0.50, RMSE = 10.09 kg) were obtained by directly estimating AGB change from the point cloud when compared to field measured AGB change. Guimarães-Steinicke et al. [115] also investigated grass biomass temporal change in Thuringia, Germany. They used 25 TLS derived metrics to model biomass [115].

While most of the reviewed studies focused on the estimation of tree biomass, only few papers dealt with shrub biomass. Shrub biomass was investigated by Greaves et al. [59] in the arctic tundra. They employed a leave one out cross-validation technique of the canopy volume against the harvested shrub biomass to establish predictive relationships between the TLS canopy volume and harvested shrub biomass and between Airborne Laser Scanning (ALS) canopy volume and TLS derived shrub biomass estimates. The TLS produced more accurate predictions of shrub biomass as compared to the ALS (TLS R<sup>2</sup> = 0.78 and ALS R2 = 0.62). Li et al. [181] merged TLS and ALS to estimate the biomass of sagebrush in southwest Idaho, USA by building predictive models. TLS derived biomass was used for calibration of a biomass predictive model with an R<sup>2</sup> of 0.87. The TLS derived shrub volume was compared to the destructive shrub biomass measurements (R2 = 0.88) [181].

#### 3.8.2. LAI

LAI is defined "as the ratio of leaf area to per unit ground surface area" [39,130]. A study was conducted in the savanna biome by Béland et al. [39]. They applied a voxelbased approach making use of the contact frequency method on six trees. To this end, they employed two multi-scan setups with 2 and 3 positions and compared the TLS derived estimates against direct measurements obtained from leaf harvesting. The leaf area ranged from 30 to 530 m<sup>2</sup> whilst crown LAI was 0.8 and 7.2 with a general bias of 14%. To ensure an accurate extraction of LAI, the green and non-green components were distinguished in the TLS data based on a separation threshold applied to the return intensities [39]. Olsoy et al. [182] estimated the LAI of sagebrush from TLS data with 2 scan positions and found that the estimated LAI was highly correlated with their reference field data (R<sup>2</sup> = 0.78).

#### 3.8.3. Basal Area

Basal area is defined as a "cross-sectional area of all stems in a plot measured at breast height and usually expressed per unit of the land" [183]. Five studies estimated basal area

from TLS point clouds, with four being conducted in a temperate broadleaf and mixed forest [75,78,85,184] and one from the boreal zone [48]. No study has been published yet for the savanna biome. Yrttimaa et al. [48] and Tansey et al. [85] derived basal area from TLS point clouds by first calculation basal area for a single stem, then scaling it to one hectare. The TLS derived basal area was compared to the field derived basal area which was derived by considering the cross-sectional area of a tree to be circular [48]. The R<sup>2</sup> observed were 0.87, 0.86 and 0.63 for a plot radius of 6, 11 and 16 m respectively [48]. Côté et al. [184] used an architectural model called LiDAR to Tree Architecture (L-Architect) by reconstructing forest stands with scans of individual trees. The L-Architect derived basal area was compared to the measured values of basal area to assess accuracy. A relative difference of −0.01% was obtained between the L-Architect and field measured basal area [184]. A methodology based on voxel data structure was applied to derive the basal area of individual trees in a heterogeneous vegetation structure in Washington Park Arboretum [75]. The TLS derived basal area was derived by fitting a 3D cylinder primitive to the trunk assuming that the trunk was a geometric circle. The TLS derived basal area was compared to the manually measured basal area. The TLS derived basal area was underestimated compared to the field-based methods for conifers, deciduous and mixed stands (Conifer: Field = 12.3116 m2, TLS = 6.9217 m2; Deciduous: Field = 1.4275 m2, TLS = 1.0172 m2 & Mixed: Field = 1.8431 m2, TLS = 1.575 m2) [75]. Aijazi et al. [78] used a super voxel segmentation method to estimate basal area and concluded that parameter estimation error ranged from 1.6 to 9% from TLS point clouds.

#### 3.8.4. Stem Density

Stem density is defined as "the number of trees per hectare" [185]. Three of the reviewed studies estimated stem density, with one being performed in a savanna ecosystem [92]. Stem density is derived from the number of stems identified in the scanned area and comparing it to the field derived stem density [48,85,92]. Yrttimaa et al. [48] varied the plot size in order to investigate its effects on stem density. The smaller the plot size, the higher the correlation of the TLS-derived stem density with field measurements (r=6m|R<sup>2</sup> = 0.77; r = 11 m | R<sup>2</sup> = 0.67; r = 16 m | R<sup>2</sup> = 0.38). Muir et al. [92] determined stem density from TLS point clouds as well. They achieved an R<sup>2</sup> of 0.1 and reported a large RMSE of 130.9 trees/hectare that was attributed to plot complexity.

#### **4. Discussion**

#### *4.1. Geographical Trends in the TLS Literature*

From a global perspective, most TLS studies have been conducted in Europe, North America and China (Figure 3). Despite their large geographical extent in Africa, South America and Australia, savannas have rarely been investigated in these areas. In the early 2000s, TLS-based vegetation studies were mainly conducted in Europe and North America (Figure 2). From 2009 onwards, TLS studies expanded into the rest of the world, including Africa, Asia, Australia and South America. An increasing collaboration among TLS research groups from different parts of the world can explain this observation [57,100,147,171]. To further advance TLS research in savannas and to learn more about their ecology, building scientific capacity in financially constrained countries is key.

Savannas cover ca. 50% of Africa's landmass. Yet, only 6% of the reviewed studies were conducted in Africa. Of those studies, the majority focused on South Africa. Only a few articles presented results for other African countries, including Mali [39], Ethiopia [186] and Gabon [70]. One reason for this observation is that TLS research in Africa is mainly constrained by financial limitations. The acquisition and processing of TLS data usually require expensive equipment and software [22]. Moreover, it necessitates trained personnel, which is limited in most third-world countries [176]. With the advent of low-cost terrestrial laser scanners [65,187,188], more TLS research can be expected in the future for Africa.

Since savannas cover large geographical areas, TLS data can usually not be used by itself to study an entire region due to its limited footprint. Instead, an integration of TLS

point clouds with other geospatial data is necessary for large-scale assessments. In the past, the creation of wall to wall maps of savannas has been achieved by combining TLS with L-band SAR [62] and airborne LiDAR [69]. Currently, the existence of freely available satellite data from Sentinel-1/2 and Landsat 5–8, as well as coarser resolution imagery from MODIS [70] can resolve upscaling issues from fine spatial resolution from TLS to medium resolution from Landsat and Sentinel then to coarse resolution using MODIS. The planned launch of BIOMASS by the European space agency and Multi-footprint Observation LiDAR and Imager (MOLI) by JAXA after Global Ecosystems Dynamics Investigation (GEDI) and Ice Cloud & Elevation Satellite-2 (ICESat-2) will present an opportunity for landscape-level studies of the savanna [60]. Based on these and similar data sets, TLS point clouds can be used in wide-area mapping for calibration and validation purposes [60].

#### *4.2. Challenges & Opportunities*

Single and multiple scans (five, four, three and two positions) were found to be the most commonly used TLS setups in the reviewed literature. Single scan approaches have been used successfully to estimate DBH [138,166], height [92] and to detect stems [189]. Single scan inventories received a lot of attention because they allow for fast and less labour intensive data acquisitions and do not require the registration of multiple point clouds [109,112]. However, single scan approaches are also more prone to occlusions [43]. Data from multi-scan TLS campaigns were shown to yield higher accuracies and reduced error margins [45,139]. Therefore, multiple TLS scans with larger plot sizes (>15 m in radius) provide the only opportunity to adequately capture the heterogeneous nature of savannas and to characterize the full spectrum of savanna vegetation structure [92]. Depending on the parameters to be extracted, the requirements to be met and the available resources, the choice between a single- or multi-scan TLS setup should be made carefully. For example, long range single scans give reliable data in the dry season (leaf-off) and not in the wet season (leaf-on) [30].

DBH is one of the most reliable proxies for biomass [190]. This is why it is also the most frequently measured quantity not only by foresters, but also with regard to the reviewed literature. It has been shown that DBH can be delineated automatically based on multiple TLS scans [40]. However, these techniques have not yet been adjusted or tested for more complex vegetation structures [88]. To date, most algorithms to retrieve DBH (and also stems and stem curve profiles) assume a circular or cylindrical stem geometry. However, these rather simple geometries are not well-suited to describe the short stature vegetation and deformed trees that are typically found in savannas [40,191]. While a number of DBH extraction techniques (e.g., RANSAC algorithm [186], point cloud slicing at 1.3 m [75], least-squares circle fitting [52] and Hough transformation [86] exist), they have only been developed and demonstrated in temperate and boreal forests. For savannas, only the circle fitting method has been explored yet [92].

Height was extracted from TLS point clouds in the savanna, mostly through the employment of CHMs to achieve the goal. Muir et al. [92] obtained an R2 of 0.97 and 0.38 for TLS-based estimates of tree and shrub heights, respectively. Small stature woody vegetation is often not well captured by CHMs, which leads to an underestimation of its height when interpolation is used to create the canopy surface [170]. This issue has not yet been adequately addressed, and more accurate height assessment techniques still need to be developed [170]. Potential ways to achieve this would be to elevate the scanner above the vegetation layer or place the scanner on a suitable vantage position [30,92], to capture TLS data of shrubs with higher point density (e.g., 16 pts per m2, [170]) and to distinguish between trees and shrubs in the TLS point cloud and to create two corresponding tree- and shrub-related CHMs that can be analysed separately [92].

Other methods for vegetation height extraction are available, including the combined analysis of TLS and airborne (also UAV) LiDAR data [140,192], the highest z-coordinate (height is derived by the highest value of the return) [51], the vertical distance between the highest point and stem base at the ground (the vertical distance is measured from the base of the stem to the highest point on a segmented point cloud) [30] and cylinder fitting (a cylinder is fitted on a segmented point cloud and the height is denoted as the highest point in a vertical cylinder) [58,145]. While they are also not without errors, these techniques have neither been developed, adjusted nor evaluated for savannas. Known problems, e.g., uncertainties from treetop shading [86], are especially common in the riparian zones of savanna ecosystems.

The extraction of LAI using TLS point clouds is not very common across all biomes [22] including the savanna. Calders et al. [79] argue that caution is needed when using LAI estimates and the variable that is measured from ground-based sensors is effective Plant Area Index (ePAI) and effective Wood Area Index (eWAI). They recommend that TLS measurements are ideal for estimation of gap fraction, ePAI and eWAI as compared to instruments that are dependent on certain illumination conditions [79]. Gap fraction is an indirect way of estimating LAI and TLS is recommended for the estimation of gap fraction to derive ePAI and eWAI [79]. It has been used successfully in the past to estimate LAI [193]. TLS intensity data can also be used to compute LAI. They allow for differentiating between photosynthetically active (green) from non-active (woody) vegetation parts [39,54]. However, they can also be compromised by instrumental, atmospheric and object properties as well as the scanning geometry [82,129]. In order to reliably make use of TLS intensity data for vegetation applications, conversion of instrumental intensity into the proportional or equal target reflectance is required [129,135].

While tree crown attributes have been successfully derived from TLS point clouds, more research needs to be directed towards crown change and shrub parameters. With regard to the former, the method by Olivier et al. [50] needs to be tested in the savanna biome. Vierling et al. [170] suggested the use of vantage points, high cliffs or tall tripods in order to achieve steep viewing zenith angles. These enable the acquisition of enough returns from both vegetation and ground. Following this procedure, reliable construction of DEMs and characterization of vegetation structure is ensured. Other ways to capture more accurate canopy details are the acquisition of full-waveform TLS data [170], the separation of trees and shrubs on the basis of height thresholds and noise filtering [92,148,170] and also separation through the application of reflectance thresholds and LiDAR return deviations [69].

Various reconstruction algorithms are available to extract plant volume from TLS data [48,54,96,101]. A limitation common to many of these techniques is that they rely on simple geometries (e.g., cylinders) to model vegetation components. This usually leads to an overestimation of volume due to the tapering of branches and the complexity of the tree architecture [176]. Another challenge lies in the exact evaluation of TLS-based volume estimates. True validation requires destructive sampling which is seldomly possible and desired [176]. Despite these limitations, volume reconstructions from TLS point clouds are key to improving tree allometries, non-destructively estimating and monitoring biomass change as well as assessing structural gains and losses in savanna ecosystems [30,175].

AGB is a key parameter in the study of vegetation structure and crucial indicator of terrestrial carbon pools and productivity [140,190]. Accordingly, it is also a popular secondary vegetation attribute in the TLS literature. Volume estimates form the basis of AGB assessments with TLS data, which are less prone to errors when compared to the simple use of allometric equations [40,176]. QSMs are most frequently employed for estimating AGB from TLS point clouds. However, reconstruction quality still suffers from certain errors and uncertainties. These are related to the modelling of complex tree crowns and smaller branches, the optimal use of wood density values (which vary both within and between species and regions), the use of simple reconstruction geometries such as cylinders, the separation of leaf and woody components as well as the validation of AGB values [71,176]. Despite these limitations, volume-based AGB estimates do not require information on DBH and show better agreement with reference data. As a nondestructive proxy for estimation of AGB, they should be further explored in the savanna context [40,176].

Using TLS point clouds, branch level parameters from a much larger percentage of a whole tree can be extracted [194]. Combined with other tree attributes, they hold the potential to calculate whole tree carbon and water use [194]. Especially reconstruction algorithms like QSMs [177] have enabled their retrieval. To date, there are very few studies on this topic and none of them focused on savannas. QSMs have been reported to yield accurate results when modelling thick branches [195]. Accuracies usually drop for smaller branches, particularly in cases where leaf-on scans are employed [195]. In order to improve QSMs with regard to branch parameter retrieval, they should be able to correctly separate the point cloud into leaf and woody components [71].

Shrub volume and AGB has also been estimated using TLS-based reconstruction algorithms. Greaves et al. [142] employed volumetric surface differencing and voxel counting in the arctic tundra whereas Olsoy et al. [54] used voxel-based and 3D convex hull techniques to assess the volume of shrubs. The existing methods hold the potential to be successfully transferred and adjusted to savanna ecosystems. This would allow for quantifying shrub layer biomass and help making more accurate AGB assessments overally. Moreover, these integrated biomass estimates could be used as calibration and validation reference data for mapping and monitoring savannas at larger spatial scales [142].

#### *4.3. Future Outlook*

The accurate extraction of vegetation parameters from TLS data is a challenging task. This review has shown that there are still a lot of research gaps that need to be addressed. This is especially true for savannas, which are in many ways much different from those biomes that have so far gained more attention in the literature. Savannas are heterogeneous and dynamic ecosystems that feature unique and complex mixtures of trees, shrubs and grasses. This requires novel developments, (re-) adjustments as well as extensive tests and experiments when it comes to TLS algorithms and applications.

Savannas are regularly subject to disturbances from fire, herbivory, droughts and humans. These disturbances lead to vegetation structural changes (including the loss of large trees [66,196] that can affect savanna ecological functioning [196]. As a future avenue for mapping and monitoring savanna vegetation, the acquisition and analysis of multitemporal TLS point clouds should be considered. While research on this subject is still in its infancy (e.g., [49,50,115]), there is no doubt that it holds a large potential for characterizing and understanding savanna dynamics [70], including phenological processes (i.e., leaf-on vs. leaf-off conditions) [54]. Among others, this is because of the non-destructive nature of TLS [49,50] and its ability to repeatedly capture 3D vegetation information with great detail [58]. Levick and Asner et al. [196] quantified the loss of large savanna trees and observed declines of more than 30% for some sites. However, since they were using ALS, below-canopy treefall and canopy-related losses (e.g., branch breakages) could not detected [196]. Such limitations can be addressed by acquiring and investigating time series of TLS scans. By developing suitable change detection algorithms, multi-temporal TLS will enable accurate and regular accounting of carbon losses and gains in savannas.

To date, woody AGB mapping in the savanna has mainly been focusing on tree biomass. Despite being a key component of savannas, shrubs have not been considered in this regard. Future research should therefore focus on an accurate and more holistic assessment of woody AGB in savannas, including trees, small understory trees and shrubs [65]. One pathway to achieve this would be a combination of TLS and airborne LiDAR data. TLS is well-suited to characterize understory vegetation structure in terms of both trees and shrubs [64,91]. Especially regarding the latter, existing algorithms that have initially been developed for sagebrush and arctic shrubs [54,110,142] could be tested and refined. Airborne (including UAV) LiDAR offers 3D point clouds with a bird's eye perspective of the savanna. While it has been found to be of limited use in shrub-dominated ecosystems [170], it is superior for tree height and canopy structure mapping [27]. By combining TLS and airborne LiDAR data, synergies between both technologies are exploited while their individual limitations can be overcome. To increase spatial coverage

and to create wall-to-wall AGB maps, optical (e.g., Landsat and Sentinel-2) and SAR (e.g., Sentinel-1, ALOS PALSAR-1/2, BIOMASS) satellite time series should be integrated into such a monitoring scheme.

From a technical perspective, TLS-based assessments of volume and AGB have not adequately taken into account the different types of woody vegetation in the savanna. An opportunity to improve them would be the separation of TLS point clouds into individual trees and shrubs. Afterwards, specific 3D models tailored to each savanna life form as well as dedicated noise filtering algorithms could be applied to yield more accurate reconstruction results [64,92,170]. For instance, QSMs could be employed to extract volume information from tree point clouds while convex hulls or voxels could be used in the case of shrub point clouds. To this end, methods exist that have been tested on trees in temperate [37], boreal [96,175] and tropical forests [101,151] as well as shrubs in the temperate regions [54,181]. Future studies in the savanna need to transfer, train and validate these techniques.

TLS intensity data have already been demonstrated to hold potential for vegetation analyses in a few studies. However, only limited overall evidence can be found in the literature so that a number of unresolved scientific questions still remain. In particular, more research is needed on the improved retrieval of the LAI with TLS intensity [80,150]. Among others, this could be achieved by a better separation of leaf and wood material [176]. Other challenging tasks based on TLS reflectance data comprise the study of seasonal changes (i.e., leaf-on vs. leaf-off vegetation properties) [54] and leaf biochemical parameters to understand the photosynthetic capacity of savanna vegetation [27].

#### **5. Conclusions**

3D data from TLS can make a vital contribution to mapping the current state of savanna vegetation, monitoring and quantifying its dynamics over space and time and, in this way, provide a better understanding of the underlying ecological factors and processes as the very basis for sustainable savanna management and long-term conservation. This study reviewed the scientific literature on the use of TLS for extracting vegetation parameters with a special geographic focus on savannas. The current state of the art was presented and methods relevant for TLS research in the savanna biome were identified and discussed. To achieve this, a comprehensive quantitative and qualitative meta-analysis of 113 research articles was conducted using 18 attributes. It was found that only few studies have dealt with TLS data for vegetation analyses in the savanna. Among others, research gaps exist with respect to a better separation of leaf and woody components in TLS point clouds, the accurate 3D reconstruction of savanna vegetation volume and biomass, more holistic assessments of woody AGB by considering both trees and shrubs as well as the application of TLS time series data for multi-temporal assessments of vegetation structural change, including studies quantifying seasonal differences in vegetation dry and green biomass. Future research should thus concentrate on extending previous experiments to savannas as well as developing novel information extraction routines for these.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2072-429 2/13/3/507/s1. Table S1: Full list of reviewed research articles.

**Author Contributions:** Conceptualization and methodology by T.P.M. and C.T. Meta-analysis by T.P.M. with the help of C.T. Draft writing by T.P.M. and C.T. Reviewing and editing by C.T., J.B., J.S., C.S. and T.P.M. Supervision by C.T., C.S. and J.B. TLS data used in Figure 1 was acquired and processed by J.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was made possible by funding from the Deutscher Akademischer Austauschdienst: DAAD Ref No. SPACES II.2 CaBuDe 57531823, Bundesministerium für Bildung und Forschung (BMBF): South African Land Degradation Monitor: Grant No. 01LL1701A, Bundesministerium für Bildung und Forschung (BMBF): Ecosystem Management Support for Climate Change in Southern Africa (EMSAfrica) Project: Grant No. 01LL1801D.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data is contained within the article or supplementary material. The data presented in this study are available in Supplementary Material Table S1.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


## *Article* **Exploring the Variability of Tropical Savanna Tree Structural Allometry with Terrestrial Laser Scanning**

**Linda Luck 1,2,\*, Lindsay B. Hutley 1, Kim Calders <sup>3</sup> and Shaun R. Levick <sup>2</sup>**


Received: 31 August 2020; Accepted: 19 November 2020; Published: 27 November 2020

**Abstract:** Individual tree carbon stock estimates typically rely on allometric scaling relationships established between field-measured stem diameter (DBH) and destructively harvested biomass. The use of DBH-based allometric equations to estimate the carbon stored over larger areas therefore, assumes that tree architecture, including branching and crown structures, are consistent for a given DBH, and that minor variations cancel out at the plot scale. We aimed to explore the degree of structural variation present at the individual tree level across a range of size-classes. We used terrestrial laser scanning (TLS) to measure the 3D structure of each tree in a 1 ha savanna plot, with coincident field-inventory. We found that stem reconstructions from TLS captured both the spatial distribution pattern and the DBH of individual trees with high confidence when compared with manual measurements (R<sup>2</sup> = 0.98, RMSE = 0.0102 m). Our exploration of the relationship between DBH, crown size and tree height revealed significant variability in savanna tree crown structure (measured as crown area). These findings question the reliability of DBH-based allometric equations for adequately representing diversity in tree architecture, and therefore carbon storage, in tropical savannas. However, adoption of TLS outside environmental research has been slow due to considerable capital cost and monitoring programs often continue to rely on sub-plot monitoring and traditional allometric equations. A central aspect of our study explores the utility of a lower-cost TLS system not generally used for vegetation surveys. We discuss the potential benefits of alternative TLS-based approaches, such as explicit modelling of tree structure or voxel-based analyses, to capture the diverse 3D structures of savanna trees. Our research highlights structural heterogeneity as a source of uncertainty in savanna tree carbon estimates and demonstrates the potential for greater inclusion of cost-effective TLS technology in national monitoring programs.

**Keywords:** allometry; biomass; carbon; cost-effective; LiDAR; TLS

#### **1. Introduction**

Savanna vegetation structure and biomass are shaped by the scale-dependent interaction of resource and disturbance drivers. Whilst limited by abiotic resource factors such as climatic conditions, moisture availability [1] and topoedaphic controls [2,3], savanna vegetation is also modulated by disturbance factors such as fire regime (frequency, severity) and storm damage, herbivory and human utilization [4–7]. The interaction and relative importance of these drivers differ at contrasting scales and their impact on above-ground biomass (AGB) can be localised and patchy, varying even at individual tree level [5,8]. A commonly referenced estimate approximates savanna vegetation at 30% of global production; however, the range of 1–12 t C ha−<sup>1</sup> yr−<sup>1</sup> net primary productivity [9] provided with this estimate indicates large unresolved uncertainty and considerable opportunity to improve

those estimates [10]. Global carbon models have long been challenged by poor representation of the structural heterogeneity in tropical savannas [11] and our understanding of savanna carbon dynamics remains limited [12–14].

At the individual tree level, resource constraints and disturbance drivers may result in strong variability of individual tree architecture and therefore AGB. Fire and herbivory are the two main drivers dramatically altering the vertical structure of savanna trees in African savannas [15,16]. In the tropical savannas of northern Australia however, termites in particular are a significant source of biomass consumption [17] and the activity of wood-eating termites results in stem and branch hollowing, increasing the susceptibility of woody vegetation to storm and fire damage. The resulting asymmetric tree architecture is particularly evident in the crown structure of large trees [18,19]. Tree crowns in north Australian savanna vegetation are estimated to account for up to nearly 45% of total tree AGB, where destructive sampling of a small number of trees (*n* = 48) has shown considerable variation in the proportion of AGB stored in the crown (Eucalypts: branch biomass 17.1% ± 7.6% and leaf biomass 3.6% ± 1.2%; non-Eucalypts: branch biomass 26.6% ± 11.7% and leaf biomass 5.5% ± 1.5%) [20]. However, allometric models commonly rely only on stem diameter (DBH), or DBH and height, to estimate AGB as these are relatively easy to measure in the field. Recent studies have shown that adding crown diameter as a predictive variable to estimate AGB provides more robust allometric models [21]. However, estimating crown dimensions by means of measuring vertical crown projection is generally omitted from field inventories due to time constraints and the inherent inaccuracy of such ground-based measures. Not accounting for variations in individual tree architecture could be particularly problematic in tropical savannas, where tree crown structure has the potential to vary significantly for a given DBH [22], however, the extent of heterogeneity in savanna crown architecture is yet to be established.

LiDAR scanning has become a well-recognised and rapidly developing approach to quantify forest structure in great detail, as it enables us to produce detailed 3D reconstructions of individual trees and their structure [23,24]. Of the range of LiDAR platforms, terrestrial laser scanning (TLS) has been increasingly used for investigating heterogeneity in tree and stand structure in space and time, due to its high point densities and scan collections from beneath the canopy [25,26]. TLS data allows us to extract traditional field measured variables such as DBH, stand basal area and tree canopy height with a high degree of precision [27]. In addition, TLS also enables new metrics, such as detailed vertical profiles [28] of tree and canopy structure [29–31] and individual crown-morphology to be quantified [32], thus providing a powerful tool to better quantify heterogeneity in tree crown structure in tropical savannas.

When applied in temperate Eucalypt open forest, TLS data was able to show a weakening allometric relationship between DBH and AGB with increasing DBH [33]. We hypothesise that in tropical savanna, given the cumulative effect of disturbance on tree architecture, variance in crown size is a major driver of this weakening relationship and a similar trend can therefore be observed between DBH and crown area. TLS now allows us to measure large numbers of individual trees to investigate the relationship between crown size and DBH as a potential source of uncertainty impacting the stability of savanna tree allometric relationships.

A major limiting factor to the adoption of TLS for vegetation surveys is the budgetary implications of purchasing high-quality scanners (>USD 100,000) generally used for geomorphological and ecological surveys [34]. Such scanners generally have small beam divergence, high signal-to-noise ratio and in-built GPS systems. These scanners also have a minimum scanning range of several hundred meters. However, for the purpose of vegetation surveys, scan points are generally located less than 50 m apart to avoid data loss through occlusion. Smaller, less precise and lower-powered TLS units generally used in architecture can technically obtain point clouds within that range, however, their efficacy in surveying savanna vegetation structure has not yet been investigated. While increasingly adopted in environmental research, the acquisition cost is a major barrier to the adoption of TLS scanning

in vegetation monitoring programs, and the advantages of TLS need to be replicable using more affordable options to enable wider adoption [35].

As such, the aim of this study is to explore the potential of lower-cost (~USD 20,000) TLS in capturing diverse tree architecture in tropical savanna ecosystems. Using this technology, our primary research objective is to determine the degree of variability present in the tree crown area for a given DBH across many hundreds of individual trees. We approached this by coupling a lower-cost, light-weight TLS scanning system with traditional field-inventory surveys at a research site in the tropical savanna of northern Australia.

#### **2. Materials and Methods**

#### *2.1. Study Site*

The study was carried out in August 2018 in open woodland savanna at the Terrestrial Ecosystem Research Network (TERN) Savanna Supersite at Litchfield National Park in the Northern Territory, Australia [36] (Figure 1). Litchfield National Park (13.17◦S, 130.79◦E), located 100 km south of Darwin, is subject to a wet–dry tropical climate with the majority of annual rainfall (1750 mm average annual rainfall at Buley rockhole, 7.4 km [37]) occurring between November and April. The survey area was a 1 ha plot and is dominated by *Eucalyptus miniata* with moderate stem density of 492 trees ha−<sup>1</sup> (>0.05 m DBH) and DBH values (up to 0.49 m). Most crowns are not overlapping, and leaves of local eucalypt species are largely near-vertically aligned, allowing good TLS coverage of tree crowns.

**Figure 1.** Study area. Left: The north Australian tropical zone spans across Western Australia (WA), the Northern Territory (NT), and Queensland (Qld). Right: canopy height model of the 1 ha study area from TLS (scale bar in m, color scale ranging from blue for lowest points through to red for highest points) showing heterogeneous distribution of tree size classes.

#### *2.2. Data Acquisition*

A manual field inventory survey was conducted at the 1 ha site where DBH, tree species and a tree health rating were recorded for every tree with DBH > 0.05 m. The trunks of multi-stemmed trees were measured individually if DBH exceeded 0.05 m. DBH was measured over-bark at 1.3 m above ground level. Tree health was assessed by a single assessor on a scale of 1 to 5 with 1 representing undamaged and 5 representing dead [38]. To correlate individual measurements to those obtained by the TLS survey, a stem map was produced using a real-time kinematic (RTK) positioning system (Leica GS16 with SmartLink). TLS point cloud data was collected during the same field campaign using a Leica BLK360 (wavelength 830 nm, maximum range 60 m at 78% albedo, beam divergence 0.4 mrad, range accuracy 4 mm at 10 m and 7 mm at 20 m [39]) in high point density collection mode (resolution 5 mm at 10 m, scan size approx. 65 mio points) taking less than four minutes for each scan (Figure 2). The plot was scanned in a regular grid fashion from 25 points with 25 m spacing between scan locations to minimise occlusion. A reflector was placed to the north of each scan to aid manual visual alignment of individual scans. Scanning took place in the morning and later afternoon to avoid strong midday breezes and was completed in one day.

**Figure 2.** (**a**) Leica BLK360 set up within the study area; (**b**) *Eucalyptus tetrodonta* sub-sampled to 2 points cm<sup>−</sup>2, point size 1. The crown is irregularly shaped with evidence of storm and/or fire damage.

#### *2.3. Point Cloud Processing*

The raw data was sub-sampled into a uniform point-spacing using an octree filter of 0.01 m. The individual point clouds were visually aligned using CloudCompare v 2.10.2 (Zephyrus), and co-registered using the Multi-station Adjustment module in RIEGL RiSCAN PRO v 2.9. The co-registered point cloud was georeferenced to the WGS84/UTM52S coordinate system using CloudCompare where RTK points of the flux tower and anchor points for its three guy ropes were used as a reference. Point cloud segmentation and extraction of tree attributes (DBH, tree height, crown area) was performed using LiDAR360 v 4.0. Spatial analysis was performed using QGIS v 3.8.0 and statistical analysis using RStudio. An overview of the workflow is shown in Figure 3 and described in detail below. Crown area was chosen as a representative measurement for crown size as it is readily replicated and commonly used in airborne LiDAR surveys (e.g., [40,41]) and is therefore scalable.

**Figure 3.** Workflow diagram for point cloud processing.

#### 2.3.1. Point Cloud Segmentation

The co-registered point-cloud was filtered for outliers and ground points were classified based on morphometric properties. A digital elevation model (DEM) was generated from the ground points and was then used to normalise the remaining point cloud to height above ground level. A mean shift algorithm, implemented within LiDAR360 TLS Forest Point Cloud Segmentation tool, was used for automatic segmentation of the normalised point cloud into individual tree clouds. The Individual Tree Editor tool was then used to manually edit falsely classified trees (Figure 4). The resulting tree attribute dataset (n = 976) was inspected for errors and a DBH to height ratio was used to flag any potentially false DBH measurements. Flagged entries were inspected in the segmented point cloud and trees where DBH could not be reliably established due to foliage present at 1.3 m above ground were excluded from the dataset. This process led to the exclusion of 169 unusable segments. Individual tree structural attributes were recalculated after the completion of the manual editing and quality control steps. The automated segmentation and manual editing tools allow for a streamlined and user-friendly segmentation process.

**Figure 4.** The study site (1 ha plot) segmented into 976 individual trees. Colours are assigned randomly; automated segmentation was followed by fine scale user editing.

#### 2.3.2. Matching of Field and TLS Survey

The manual and TLS survey data were overlaid as shapefiles in QGIS and corresponding trees were matched spatially. This was achieved by applying a 1 m buffer around points in the manual field survey data set. Offsets in the stem manual positioning due to signal drift was accounted for by moving TLS points into the buffer of corresponding manual points where a match could be confidently observed. Of the 492 trees manually surveyed, 450 could be matched in the TLS dataset. Attributes were then joined by spatial location and checked for duplicate records (indicating joins of two or more trees from one layer to one tree of the other).

To compare the spatial distribution of trees captured in both surveys, the plot was overlaid with a 10 m and a 20 m grid. For the TLS data, DBH values below 0.05 m DBH (less 0.07 m mean absolute error (MAE)) were excluded using the query builder. The number of trees located within each grid cell were counted for both surveys. The resulting polygons were converted to a raster and the difference between TLS and field tree count calculated using the raster calculator.

#### *2.4. Statistical Analysis*

All statistical analyses were performed in the R environment [42] v 3.6.0. Correlation of field and TLS-derived DBH was performed using linear regression (Figure 7). Visual inspection of the relationship between DBH and crown area suggested an inflection point beyond which the relationship between crown area and DBH changes. The location of this point was identified with a change-point regression using the chngpt package [43]. As the data are heteroscedastic (the variability of the outcome is not constant across the range of the predictor, violating the assumption of linear regression), a generalised least squares (GLS) was used to interpret the correlation before and after the change point (Figure 8). A paired sample *t*-test was then used to assess how closely spatial distribution of detected trees matched between the surveys. For each sub-plot the number of detected trees derived from each method was compared. For this analysis, trees with DBH smaller than 0.05 m were excluded from the TLS survey (less the mean average error from the field and TLS correlation), to match the sampling scope of the manual survey (Figure 6). All R code used is available in a github repository (see Supplementary Materials).

#### **3. Results**

#### *3.1. Tree Counts and Distribution Obtained from Field and TLS*

The TLS survey allowed for the capture of tree location and DBH measurements on a plot scale compatible with data collected using traditional field surveys. Visual inspection of the field inventory and TLS surveys showed similar patterns in the spatial distributions of DBH as captured by each survey method (Figure 5a,b) with a similar detection curve evident for both surveys (Figure 5c,d). Additionally, the TLS was able to detect a large number of small trees (DBH 0.02–0.04 m) which are often excluded from manual field inventory, as this sampling scale takes considerable time to undertake.

**Figure 5.** Spatial distribution and number of trees captured in both surveys conducted in August 2018. (**a**) distribution map of field survey; (**b**) distribution map of terrestrial laser scanning (TLS) survey, red denotes trees with stem diameter (DBH) < 0.05 m not captured in the manual survey (*n* grey = 449, *n* red = 358); (**c**) Histogram of field survey; (**d**) Histogram of TLS survey. TLS surveys are able to capture small trees that are impractical to include in manual field inventories.

Comparing the spatial distribution of tree detections in each survey highlights the potential impact of edge effect on small sub-plots (10 and 20 m2). Across all sub-plots there was no significant difference between the TLS and field inventory for tree detection (Paired samples *t*-test 10 m grid: t(99) = −1.815, *p* = 0.07, 20 m grid: t(24) = −1.97, *p* = 0.06) (Figure 6). However, for a number of sub-plots there was clear over- or under- estimation observed and are likely due to edge effects and/or GPS drift. Using sub-plots as simulated by the 10 m and 20 m grids may be disadvantageous, as individual trees are forced into small plots and thus not necessarily representative of tree density in the area.

**Figure 6.** Difference between TLS tree count and field tree count. The average tree count per grid cell in the field survey is not significantly different from the average tree count in the TLS survey (*p* > 0.05). For this analysis only trees with DBH > 0.05 m (−0.007 m mean absolute error (MAE)) were used.

#### *3.2. DBH Obtained from Manual Survey and TLS*

The DBH values of individual trees derived from the TLS survey were strongly correlated with those obtained manually in the field. The RTK GPS tagging of trees during the field survey ensured precise matching of a large sample of trees captured in the field survey to corresponding TLS measurements (two outliers were excluded from the analysis). There was a strong correlation between DBH measured in the field and DBH derived from the TLS point cloud (Figure 7), with R<sup>2</sup> = 0.98 and MAE = 0.007 m (Spearman's *ρ* = 0.93, *p* < 0.001). The MAE shows the overall average of errors while the RMSE (root mean square error) highlights strong outliers as it increases with variance in the frequency distribution of error magnitudes. For this model, an RMSE value higher than the MAE indicates that the main source of error likely stems from a small number of strong outliers, rather than a systematic over or underestimation of field DBH by the TLS derived DBH.

**Figure 7.** Correlation between field survey and concurrent TLS scan; ribbon showing 95% confidence interval. For every cm increase in TLS DBH there is an average 1.04 cm increase in field DBH (95% CI 1.03–1.05, *p* < 0.001) (*n* = 448).

#### *3.3. Variability in Crown Architecture*

Using TLS derived data, crown size (canopy area) was plotted as a function of DBH (Figure 8) which showed a strong positive relationship with heteroscedastic variance. Deviance from the mean increases along with maximum crown area, resulting in a non-linear relationship (Table 1). Our results align with previous observations made while calibrating the local allometric model [22]. A change-point regression shows a change point at DBH of 0.099 m (95%CI 0.056–0.14, *p* < 0.001). For trees below the change point, for every unit (1 m) increase in DBH there is a 88.55 unit increase in the crown area (95%CI 81.56–95.50); for trees above the change point, the increase in crown area per unit DBH increases to 240.87 and the 95%CI widens to 191.75–290. Trees with poor health ratings were not excluded, as this would obscure the true variability in crown size.

**Figure 8.** TLS derived crown area for a given DBH; ribbon showing 95% confidence interval. The spread of crown sizes increases with increasing DBH (*n* = 807).

**Table 1.** Variability in tree crown structural parameters as a function of DBH class, showing values of the mean and standard error of the mean.


#### **4. Discussion**

#### *4.1. Heterogeneity at Landscape/Stand Level*

Tropical savannas are disturbance-impacted ecosystems which result in stand-scale structural heterogeneity, introducing unquantified uncertainty to commonly used methods for estimating AGB that rely on small sub-samples of DBH-based allometry. When used in tropical savanna woodlands, TLS can provide accurate and detailed information that allows us to explore this source of uncertainty in allometric models. Furthermore, the large volume of samples obtainable using TLS allows us to negate further sources of uncertainty in traditional field surveys arising from the sampling limitations in the trees surveyed to represent a population when both calibrating and applying allometric models. Deriving stand size class distributions based on DBH measurements is commonplace in production

forestry through to ecological surveys of stand biomass. Manual tree measurement and obtaining the geolocation data required to study size class distribution requires considerable resources and few studies include geolocation of individual trees within a study site [44]. Instead, a common field method for estimating the size class distribution is to extrapolate from sub-plots that may be 30 m × 30 m or less in size [45]. However, in tropical savannas, the assumption of spatial homogeneity underlying this approach is unlikely to be met, given inherent spatial variability at distances of 20 to 30 m (Figure 6). TLS data allows us to extract geolocation along with DBH for individual trees and is increasingly used to investigate ecological processes such as random vs clumped tree distributions, structural change over time following disturbance [25,30,46] and over annual to decadal scales, or woody encroachment and thickening [47]. We show that TLS data obtained using a lower-cost scanner enables fast and accurate acquisition of both, full 3D tree measurements and geographic positioning required to investigate landscape scale heterogeneity in tree size class distribution (Figures 1 and 5).

#### *4.2. Accuracy of Lower-Cost TLS Scanner for Estimating DBH*

Direct comparison of DBH obtained from a manual field survey and concurrent TLS survey showed a very close fit with a mean absolute error (MAE) of less than 0.01 m, suggesting that lower-cost TLS data can provide traditional forestry and ecological metrics such as stand basal area and size class distributions with high precision. Our analysis of DBH obtained in the field as compared to TLS (Section 3.2) indicates random variation as the main source of error rather than systemic over or under-prediction of the model (Figure 7). Random errors are common in manual data collection [48,49], while errors in TLS measurements are commonly caused by wind, random errors in point acquisition, registration error or occlusion. Errors in data processing occur when measuring DBH, as the model assumes a perfect circle, which is often not the case. Our study design aimed to minimise TLS acquisition error and occlusion. Our data shows an RMSE of 1.05 cm, representing a 50% reduction compared to a similar study in temperate Eucalypt open forest where a minimum distance of 40 m between scan points was used [33]. However, the two studies also differ in stand structure, equipment and software used and the impact of either of those factors is yet to be quantified. The accuracy achieved in this study also reflects results from other vegetation types, such as mixed forests in Austria [50] and Germany [51], deciduous forest in Iran [52] and India [53], Chinese fir plantations [54], or tropical forest of Malaysia [55]. However, none of these are structurally representative of tropical savannas and the suitability of TLS for obtaining DBH was yet to be established. We show that lower-cost TLS can replicate field measured DBH and tree height, establishing the suitability of TLS to reliably capture traditional measures of vegetation structure in tropical savanna vegetation.

#### *4.3. Heterogeneity at the Individual Tree Level*

Variability in tree crown structure introduces poorly quantified and potentially significant uncertainty in allometric models predicting AGB based on DBH. Previous research shows that in north Australian savannas the proportion of biomass stored in tree crowns can vary between 12% and 45% of total tree AGB (mean ± SD of the branch and foliage components) [20]. However, this estimate is based on a very small sample size of 48 destructively harvested individuals and is unlikely to capture the true extend of variability. Our TLS analysis of over 800 individual trees confirmed and quantified a significant spread of crown area for a given DBH, increasing with DBH (Figure 8). This is in contrast to previous research published exploring the relationship between DBH and canopy diameter [56] where a linear relationship (R<sup>2</sup> = 0.63) was established across a much larger area. However, investigation of the raw data shows a weaker relationship for trees surveyed within 150 km of Litchfield NP, indicating that local-scale variability arising from disturbance might exceed that of regional-scale trends. Furthermore, structural variability is not only found in overall crown size, but also architecture (Figure 9). While the addition of crown diameter as a predictive variable improves allometric models [21], the use of canopy area as a function of canopy diameter assumes a circular crown shape. Using lower-cost TLS we were

able to highlight the impact of disturbance on tree architecture (Figure 9) and the potential effect on the stability of allometric relationships introduced by the assumption of homogeneity.

**Figure 9.** Two examples of typical mature *Eucalyptus miniata* individuals, both with a DBH of 0.38 m highlighting the dramatic difference in canopy structure, volume and biomass. Estimates of carbon storage of the left individual would be significantly overestimated using a manually measured DBH with allometry applied.

#### *4.4. Further Applications—Capturing Change and Irregularities*

The effect of fire as disturbance driver on the understory, including tree seedlings and saplings, and woody shrubs is a critical ecological process in savanna [57,58], however, detailed inventories of woody understory vegetation required to quantify this effect are not readily quantified using standard field inventory. TLS provides the potential for a more accurate inventory of understory density if high-resolution point clouds are generated. The ability to resolve the diameters of smaller stems (less than 0.05 m) provides crucial information on the recruitment potential of a stand and/or evidence of recent disturbance. This sampling detail is undertaken in some traditional inventory methods but can be extremely time consuming and costly. Furthermore, allometric equations have been calibrated for small to medium-sized trees with a small fraction of large individuals and are unlikely to be representative of seedlings and saplings.TLS allows us to capture understory vegetation such as seedlings, saplings and small shrubs (Figure 5) and use alternative volume-based analyses to estimate AGB. From the 3D point clouds, algorithms such as voxel or convex hull based approaches [59,60] can be used to estimate AGB of woody vegetation where metrics like DBH cannot easily be measured, or suitable allometric models are not available. When applied as successive scans, TLS provides

detailed observations of structural changes at a stand level, such as mortality and recruitment, or at individual tree level, such as damage to branching structure and crown development [61,62]. The ability to quantify recruitment and mortality opens up the opportunity to explore disturbance dynamics such as fire impacts or storm damage impacts within the regeneration niche at unprecedented detail.

Multi-stemmed trees and shrubs can be challenging to survey manually, especially when a forking height is located below or at the designated DBH measuring height. To avoid this issue, some studies opt to measure diameter closer to the ground, e.g., 0.1 m or 0.15 m above ground level [63,64]. However, such a strong deviation from the norm will prevent comparison or integration of data from multiple studies and requires further calibration of specialised allometric models. TLS is not restricted to single-point measurements and can be used to extract a range of data to estimate AGB. For instance, Quantitative Structure Modelling (QSM) has gained popularity as an alternative to DBH-based allometric models. Individual tree AGB is estimated based on wood density and point cloud generated tree volume, avoiding the need for allometric relationships and standardised measurements [33]. The TLS approach provides the freedom to adapt sampling strategies to include stand elements that will not conform to traditional survey methods and tree allometric models.

Accounting for loss from termite consumption in susceptible tree species remains a challenge for biomass estimates based on crown and stem volumes [65]. However, if hollowing functions accounting for this loss can be developed, the use of TLS offers considerable advantages over manually sampled inventories given the sampling speed, precision, the range of metrics that can be estimated and the ability to accurately quantify structural change over time.

#### *4.5. Negating the Cost Barrier*

TLS provides clear benefits over traditional manual field surveys and is increasingly utilised in the research community. However, capital cost has been a major barrier to the wider adoption of this technology, constricting the application of research outcomes by the wider community [34,35]. We show that rapid developments in commercial scanners now enable us to cost-effectively acquire point clouds of sufficient quality for detailed vegetation surveys. We demonstrate the benefits of TLS to vegetation monitoring programs by combining traditional and new measurements obtained through TLS to explore a major source of uncertainty in traditional allometric models when applied in frequently disturbed woodlands such as tropical savannas. We also highlight and discuss further advantages of TLS to capture stand elements traditionally difficult to capture using field inventories. Limitations apply to the use of lower-cost scanners, particularly in regards to their range and therefore ability to survey trees above a height of 50 m [35], however, we have demonstrated that such commercial options provide data of sufficient quality to make state of the art vegetation survey methods more accessible.

#### **5. Conclusions**

Tropical savannas have been challenging to quantify in global carbon models. We show how heterogeneity in savanna tree structure introduces uncertainty into established protocols for estimating AGB and therefore carbon content across spatial scales. At plot scale this heterogeneity leads to sample plots rarely being representative of a larger area, while at tree scale the strength of DBH as an indicator of tree AGB is diminished by strong variability in tree architecture. Furthermore, the traditional exclusion of small understory trees (e.g. <0.05 m DBH) from field inventories omits a potentially important indicator for vegetation and carbon dynamics, including recruitment and regeneration. TLS provides a more holistic way to retrieve quantitative structural measurements that directly relate to biomass properties in tropical savanna. Using TLS we were able to reliably obtain traditional field attributes such as DBH and tree height, but also access new meaningful tree attributes such as crown size and extent, and spatial distribution of understory vegetation. However, adoption of this technology in monitoring programs has been limited due to acquisition costs of high-quality TLS scanners, and skepticism remains as to the utility of TLS in vegetation surveys. We show that in

savanna vegetation a low-cost TLS scanner will produce high-quality LiDAR data and highlight the advantages and power in deriving meaningful tree attributes that will allow us to better quantify AGB at tree and plot scale. Establishing sampling protocols for highly detailed AGB surveys with reduced budgetary requirements may enable greater inclusion of TLS in national and international monitoring programs.

**Supplementary Materials:** The following are available online at https://github.com/LindaLuck/Luck2020Exploring, Analysis scripts for data analysis carried out in the R environment.

**Author Contributions:** Conceptualization: S.R.L. and L.L.; methodology: L.L., S.L; analysis: L.L.; resources: S.R.L., L.B.H.; data curation: L.L.; writing—original draft preparation: L.L.; writing—review and editing: S.R.L., L.B.H., K.C., L.L.; visualization: L.L.; supervision: S.R.L., L.B.H., K.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by a Charles Darwin University (CDU) postgraduate RTP scholarship and a Commonwealth Scientific and Industrial Research Organisation (CSIRO) top-up scholarship.

**Acknowledgments:** We acknowledge the Ecosystem Processes facility in Australia's Terrestrial Ecosystem Research Network (TERN) and NT Parks and Wildlife. Jodie Hayward and Jon Schatz are thanked for their field support. We also thank Mirjam Kaestli for statistical support. This manuscript was greatly improved by valuable comments from Keryn Paul and Garry Cook.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Leveraging TLS as a Calibration and Validation Tool for MLS and ULS Mapping of Savanna Structure and Biomass at Landscape-Scales**

**Shaun R. Levick 1,\*, Tim Whiteside 2, David A. Loewensteiner 2, Mitchel Rudge <sup>3</sup> and Renee Bartolo <sup>2</sup>**


**Abstract:** Savanna ecosystems are challenging to map and monitor as their vegetation is highly dynamic in space and time. Understanding the structural diversity and biomass distribution of savanna vegetation requires high-resolution measurements over large areas and at regular time intervals. These requirements cannot currently be met through field-based inventories nor spaceborne satellite remote sensing alone. UAV-based remote sensing offers potential as an intermediate scaling tool, providing acquisition flexibility and cost-effectiveness. Yet despite the increased availability of lightweight LiDAR payloads, the suitability of UAV-based LiDAR for mapping and monitoring savanna 3D vegetation structure is not well established. We mapped a 1 ha savanna plot with terrestrial-, mobile- and UAV-based laser scanning (TLS, MLS, and ULS), in conjunction with a traditional field-based inventory (n = 572 stems > 0.03 m). We treated the TLS dataset as the gold standard against which we evaluated the degree of complementarity and divergence of structural metrics from MLS and ULS. Sensitivity analysis showed that MLS and ULS canopy height models (CHMs) did not differ significantly from TLS-derived models at spatial resolutions greater than 2 m and 4 m respectively. Statistical comparison of the resulting point clouds showed minor overand under-estimation of woody canopy cover by MLS and ULS, respectively. Individual stem locations and DBH measurements from the field inventory were well replicated by the TLS survey (R<sup>2</sup> = 0.89, RMSE = 0.024 m), which estimated above-ground woody biomass to be 7% greater than field-inventory estimates (44.21 Mg ha−<sup>1</sup> vs. 41.08 Mg ha−1). Stem DBH could not be reliably estimated directly from the MLS or ULS, nor indirectly through allometric scaling with crown attributes (R<sup>2</sup> = 0.36, RMSE = 0.075 m). MLS and ULS show strong potential for providing rapid and larger area capture of savanna vegetation structure at resolutions suitable for many ecological investigations; however, our results underscore the necessity of nesting TLS sampling within these surveys to quantify uncertainty. Complementing large area MLS and ULS surveys with TLS sampling will expand our options for the calibration and validation of multiple spaceborne LiDAR, SAR, and optical missions.

**Keywords:** biomass; carbon; LiDAR; TLS

#### **1. Introduction**

Savanna vegetation structure is shaped by the interaction of climate, soils and a variety of disturbance agents acting at multiple spatio-temporal scales [1,2]. The actions of fire, herbivores, termites, and cyclones have marked effects on the structure of savanna tree crowns, leading to high degree of structural diversity [3–5]. Mapping and monitoring the structure of savanna ecosystems is therefore challenging not only because of their

**Citation:** Levick, S.R.; Whiteside, T.; Loewensteiner, D.A.; Rudge, M.; Bartolo, R. Leveraging TLS as a Calibration and Validation Tool for MLS and ULS Mapping of Savanna Structure and Biomass at Landscape-Scales. *Remote Sens.* **2021**, *13*, 257. https://doi.org/10.3390/ rs13020257

Received: 16 December 2020 Accepted: 11 January 2021 Published: 13 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

heterogeneous spatial patterning, but also because of disturbance driven variability in individual tree crown structure [6]. Structural attributes such as height, canopy diameter, and projected foliage cover are commonly assessed through field inventories—but these measures are often subjective and prone to sampling errors [7]. Robust assessment of structural change in tropical savannas requires methods that can account for spatial and temporal heterogeneity, and a high level of accuracy and precision in structural attribute measurement.

Remote sensing is very attractive from a monitoring and management perspective, offering large area assessment at ever increasing spatial, temporal and spectral resolutions [8,9]. Despite many advances in spaceborne imaging technology, savanna ecosystems remain challenging from an earth observation perspective due to their inherent mix of discontinuous woody cover and herbaceous vegetation [10]. Light-detection-and-ranging (LiDAR), especially airborne laser scanning (ALS), has emerged as a prominent technology for mapping and monitoring vegetation in a variety of ecosystems across the globe [11]. ALS is particularly well suited to the measurement of savanna ecosystem structure where it can delineate subtle topographic and vegetation characteristics [12]. The inclusion of ALS in ecological studies has advanced the field of savanna ecology, revealing new insights into the controls and drivers of carbon storage [13,14], and providing pathways for understanding how vegetation structure influences the ecology and diversity of various fauna [15,16].

Nonetheless, despite the many advantages of ALS, certain ecological processes operate at even finer spatial scales than can be assessed from traditional airborne platforms. Detecting stress, sub-canopy structural changes, and recruitment dynamics are all examples of processes that sit at the limits of airborne surveying. Terrestrial laser scanning (TLS) has emerged as the gold standard for fine-scale 3D reconstruction of above-ground elements, and is increasingly being used in forestry and ecological surveying [17,18]. One of the main constraints of TLS is the time-intensive nature of data acquisition and post-processing, and the relatively small spatial extent that can feasibly be covered [19]. Newer sensors are overcoming some of these barriers by increasing the speed and ranging distance, and enabling larger area coverage in suitable environmental conditions [20]. Reduced size and weight of modern LiDAR sensors has broadened the field of UAV-based laser scanning (ULS), offering a high degree of acquisition flexibility and increased spatial resolution due to lower and slower flying speeds [21–23]. In conjunction with these developments, handheld and mobile laser scanning (MLS) systems have also been gaining traction, enabling greater freedom of movement and increased speed of acquisition in comparison to TLS. MLS enables navigation in confined spaces, and the continuous mapping approach can reduce occlusion and provide thorough coverage in complex environments [24–26].

Given these recent advances in LiDAR sensor and platform technologies, choosing the right solution for a particular mapping and monitoring problem has become more challenging. The relative strengths and weaknesses of different systems needs to be explored across a range of ecosystems to best match system performance with the specific research or management objectives. The objective of this study was to assess the degree of divergence and complementarity of three different laser scanning systems for capturing the 3D structure of tropical savanna woodland structure. We explore the relative strengths and weaknesses of terrestrial laser scanning (TLS), mobile laser scanning (MLS) and UAV laser scanning (ULS) in a vegetation monitoring context, by addressing two key questions:


#### **2. Materials and Methods**

#### *2.1. Study Site*

This study was undertaken in northern Kakadu National Park in the Northern Territory of Australia (Figure 1). The climate is tropical, with maximum and minimum temperature averaging 34.1 ◦C and 22.6 ◦C respectively (Köppen climate classification = Aw). Rainfall is highly seasonal, with the vast majority of the annual 1557 mm yr−<sup>1</sup> falling predominantly in the December—March wet season [27]. The vegetation is representative of much of the Australian tropical savanna zone, with an over-storey dominated by *Eucalyptus tetrodonta*, *Eucalyptus miniata*, and *Erythrophleum chlorostachys*. The under-storey is characterized by annual grasses, particularly Sorghum species. Our data collection campaign took place on 23/24 July 2019, well into the dry season and one month after a fire passed through—the occurrence of which is typical every 1–2 years in this system [28]. The herbaceous grass layer was reduced by the fire, but small shrubs were still present.

**Figure 1.** Location of the study site in northern Australia (**a**), and an aerial view of vegetation cover within the 20 × 20 m subplots of the 1 ha study area (**b**).

#### *2.2. Field Inventory*

A series of long-term monitoring plots were established in the savanna woodlands of the Jabiru region in 2017 to provide reference sites for understanding natural system dynamics and informing the closure criteria of the adjacent Ranger uranium mine [29]. Each plot measures 100 m × 100 m, and the four corners and center are marked with permanent posts. We selected one of these sites for our study and each woody plant with a stem diameter > 0.03 m and that was > 1.5 m in height was inventoried in the field. Survey tapes were laid out in a 20 m × 20 m grid to guide the field crew, and each individual tree was: (i) tagged with a RTK dGPS (Leica GS16 with Precise Point Positioning); (ii) identified to species level; and (iii) measured with a DBH tape at 1.3 m above ground level. A total of 521 trees consisting of 572 stems were inventoried in the field.

#### *2.3. Laser Scanning*

The selected site was surveyed on 24 July 2019 with three different LiDAR systems, operating from terrestrial, mobile and UAV platforms (Table 1). Terrestrial laser scanning (TLS) was conducted from a surveying tripod at 1.8 m agl, using a Riegl VZ-2000i system with integrated RTK-GNSS (Figure 2). Sixteen scan locations were established, using a regular grid spacing of 30 m. The scanner was operated at 600 kHz with an angular sampling resolution of 30 mdeg. The scanner communicated with a RTK base station (Emlid RS2)

established at the site and operating over a LoRa network. Real-time positioning error of the scanner at each scan location was < 0.05 m.

**Figure 2.** Site conditions at the time of the LiDAR surveys. Occlusion from the grass layer was minimal due to the effects of a fire one month prior.

Mobile laser scanning (MLS) was conducted with a GreenValley LiBackpack D50. The unit consists of two Velodyne VLP-16 sensors, one operating vertically and one horizontally. The system uses Simultaneous Location and Mapping (SLAM) technology for co-registration. The MLS system was carried as a backpack by a field technician who walked slowly through the site ensuring consistent coverage and loop closure. The trajectory of the mobile unit was displayed in real time to the operator via a tablet to ensure full coverage of the site.

UAV laser scanning (ULS) point clouds were collected using the Nextcore® system (www.nextcore.co) which uses a Quanergy M8 discrete return LiDAR sensor integrated with a Spatial Dual (www.advancednavigation.com) INS including IMU and dual antenna RTK GPS, mounted on a DJI Matrice 600 Pro. The UAV was flown in a grid pattern at 3.5 ms−<sup>1</sup> at approximately 40 m agl with a line spacing of 18 m. The sensor pulse rate was 1200 kHz, which resulted in average point densities of 2000 points per m<sup>2</sup> within the one hectare plot.

**Table 1.** Instrument, surveying, and processing characteristics of the three LiDAR systems. Reported point densities are after the application of a 0.01 m sub-sampling spatial filter. Note: TLS = terrestrial laser scanning, MLS = mobile laser scanning, ULS = UAV laser scanning.


#### *2.4. Point-Cloud Processing*

Raw processing of the TLS dataset was conducted in Riegl's RiSCAN Pro software suite (v2.9). Scan data were projected into the WGS84 coordinate system (UTM Zone 53S) and filtered for noise based on reflectance and deviation characteristics. The Multi-Station Adjustment (MSA) plugin was used to finely co-register the individual scans, which were already well positioned given the integrated GNSS with real-time correction.

Raw MLS data were converted to a point cloud incorporating the IMU data from the backpack using proprietary SLAM algorithms in the LiBackpack software from GreenValley International. The point cloud was geo-referenced in LiDAR360 by detecting six pyramid structures that were placed throughout the site as ground controls. The 3D RTK-GNSS location of the peaks of each pyramid were recorded during the field campaign.

ULS point clouds were PPK geo-registered using the proprietary Nextcore® software. Each of the flight-lines underwent a statistical outlier removal (SOR) filter in CloudCompare. The number of neighbors used to compute mean distance was set to 10, and the standard deviation multiplier was set to 2. There was still some systematic misalignment of flightlines following the PPK correction. To correct for this, each flight-line was sequentially registered to its neighbor using the Iterative Closest Point (ICP) algorithm implemented in CloudCompare. The random sample limit was set to 250,000 and RMS error difference threshold set to 1 × <sup>10</sup>−08. After registration, flight-line point clouds were merged into a single point cloud.

Co-registered point clouds from each platform were exported to the ASPRS LAS format version 1.4 for further analysis. A 0.01 m spacing filter was applied to all three clouds to remove duplicate points and ensure an even distribution of points across the landscape.

#### *2.5. Structural Analysis*

#### 2.5.1. Canopy Characterization

The MLS, and ULS point clouds were finely co-registered to the TLS point cloud using the Iterative Closest Point (ICP) algorithm in CloudCompare [30]. The random sample limit was set to 500,000 and the RMS error difference threshold set to 1 × <sup>10</sup><sup>−</sup>05, with 100% overlap specified. Ground point classification was conducted with the *lasground\_new* tool in the LAStools suite [31], using -wilderness and -extra\_fine settings, and the point clouds were subsequently normalized to height above ground using *lasheight*. Canopy height models (CHMs) were generated with *lasgrid* at multiple spatial resolutions using the 95th percentile height. CHM resolution increased in 0.25 m increments from 0.25 m to 5 m resolution to enable a sensitivity analysis of the resolution at which the combinations of the three scanning systems were significantly different. Differences in canopy height distribution derived from the three laser scanning systems was tested across the range of CHM resolutions with the Kolmogorov–Smirnov test [32] in the *stats* package in R (version 4.0.3).

The 0.25 m resolution CHMs derived from the three scanning systems were used to assess the percentage of woody canopy cover at multiple height class thresholds that represented key strata in the savanna woodland (0.25 m, 3 m, 7 m, 12 m). For each height threshold, median woody canopy cover for 20 m × 20 m subplots embedded in the 1 ha site (n = 25) was compared between TLS-MLS and TLS-ULS with linear regression.

#### 2.5.2. Individual Tree Delineation and Measurement

Individual trees segmentation was applied to the point clouds derived from the three scanning systems using the approach described by [33] and implemented in Li-DAR360 (v4.1, GreenValley International) (Figure 3). Results were manually inspected post-segmentation and edited to correct any instances of over- or under-segmentation. Following this quality control checking, the point clouds representing individual trees were processed further in 3DForest (v0.5) to derive individual tree attributes—tree height, crown diameter, crown area, crown volume and stem DBH [34]. An overview of the full processing pipeline is provided in Figure 4.

**Figure 3.** Cross-section through the TLS point cloud, displaying an example of the individual tree segmentation. Color scheme is a randomised palette.

Above-ground woody biomass (AGB) was estimated at the individual tree level using a set of locally calibrated allometric equations that included the dominant species at the site [35]. AGB was derived from the field inventory as a function of stem DBH (Equation (1), and from the TLS data as a function of stem DBH and tree height (Equation (2)).

$$
\ln(AGB) = \beta 0 + \beta 1 \ast \ln(D) \tag{1}
$$

$$\ln(AGB) = \beta 0 + \beta 1 \ast \ln(D) + \beta 2 \ast \ln(H)^2 \tag{2}$$

Multiple linear regression and Random Forest modeling (randomForest [36]), conducted in R, were used to assess the relationship between individual tree biomass (as determined from DBH allometry) and additional measures of tree 3D structure (canopy height, crown area, crown width, crown volume).

**Figure 4.** Overview of the processing workflow and software packages used at different stages of analysis.

#### **3. Results**

#### *3.1. Woody Canopy Characterization*

All three LiDAR sensing systems provided high quality point-cloud data which captured the three-dimensional structure of the savanna vegetation in rich detail (Figure 5). Viewing the 1 ha plot from an oblique angle, the three point clouds showed very similar structure and it is evident that all three systems captured the canopy patterns well (Figure 5a–c). A cross-sectional slice through the point clouds (5 × 30 m) revealed the

narrower beam footprint and higher precision of the TLS instrument, with crisp trunk and branch detail, as well as good definition of material in the herbaceous layer (Figure 5d). The MLS and ULS point clouds were less resolute in comparison, exhibiting some minor noise around tree trunks and evidence of occlusion in the lower stems and herbaceous layer (Figure 5e–f).

The generation of canopy height model (CHM) raster layers confirmed that differences between the three scanning systems was minimal from a visual perspective and displayed very similar spatial patterns of canopy height (Figure 6a–c). Subtle differences were evident, however, and the canopy height distribution ((Figure 6d–f) differed significantly between all three systems at 0.25 m (KS-test, *p* < 0.001). Sensitivity analysis showed that these differences remained significant until a CHM resolution of 2 m was reached for the TLS-MLS comparisons, and 4 m for TLS-ULS and MLS-ULS comparisons (KS-test, *p* > 0.05).

Estimation of maximum canopy height was very consistent across the three systems, even at 0.25 m raster resolution. However, the MLS and ULS showed consistent bias in woody canopy cover estimation, over- and under-estimating the woody fraction respectively at multiple height intervals (Figure 7 and Table 2).

**Figure 5.** Oblique view of the point clouds derived from the three LiDAR platforms: (**a**) TLS, (**b**) MLS and (**c**) ULS. Cross-sectional view of the point clouds derived from: (**d**) TLS, (**e**) MLS, and (**f**) ULS.

**Figure 6.** Canopy height models derived from the three co-incident LiDAR datasets, plotted at 0.25 m resolution (**a**–**c**) and their height class distributions (**d**–**f**).

**Figure 7.** Canopy cover comparisons between MLS and ULS at different height class intervals, against a TLS reference. Solid black is the 1:1 fit.


**Table 2.** Summary statistics of woody canopy attributes collected by the three LiDAR systems at the 1 ha plot scale.

#### *3.2. Individual Tree Characterization*

Assessment of 3D structure at the individual tree level highlighted some further distinctions between TLS, MLS and ULS characterization of 3D structure. The TLS capture provided a very clean representation stem, branch and leaf structure, with points spaced close enough together to give the impression of a continuous surface (Figure 8a). Extracting a point-cloud slice from 1.25–1.35 m (to encompass the traditional DBH measurement height) revealed a clear ring of points that allowed for DBH measurement via cylinder fitting (Figure 8d).

The MLS scan captured the general size and shape of individual trees very well, including small under-storey plants, but loss of detail was evident compared to TLS (Figure 8b). The MLS point cloud exhibited more noise around tree stems, which is particularly clear in the thick ring of points extracted at DBH measurement height (Figure 8e). The ULS also captured the general shape of individual trees very well, and although many returns were received from tree stems the density of points in the 1.25–1.35 m range was much lower than from TLS and MLS (Figure 8c,f).

**Figure 8.** Example of a single large tree as captured by the three LiDAR systems (**a**–**c**). Horizontal slice cut through each stem at 1.25–1.35 m above ground level, showing the pattern and density of points available for DBH fitting from each system (**d**–**f**). .

Considering this, we used a broader portion of the tree stems (1–2 m) for cylinder fitting and DBH estimation. Comparison of the DBH distributions against those obtained from the field inventory showed that the TLS estimates closely mirrored the field-measured values (Figure 9a,b), but MLS and ULS estimates (Figure 9c,d) were skewed to larger values and were more normally distributed than the inverse-J patterns evident in the field and TLS data.

**Figure 9.** Distribution of DBH estimates from field inventory (**a**) and the three LiDAR approaches (**b**–**d**).

Linear regression between field-measured DBH and TLS-estimated DBH confirmed that the TLS could reliably extract individual tree DBH for a broad range of species (Figure 10). Taller tree species with larger sample sizes in the plot such as *Eucalyptus tetradonta*, *Corymbia porrecta* and *Xanthostemon paradoxus* showed the strongest fits (Figure 10), while smaller statured species with narrower DBH ranges exhibited greater error (e.g., Figure 10).

Based on the DBH representation findings discussed above, we restricted our analysis of above-ground woody biomass to the field and TLS datasets. Our field inventory returned a total of 572 stems, corresponding to 41.08 Mg ha−<sup>1</sup> of above-ground woody biomass and 20.54 Mg C ha−<sup>1</sup> from DBH-based allometry (Equation (1)). Segmentation of the TLS dataset identified almost 100 more stems (n = 669) than measured in the field survey, which were mostly in the smaller size classes (Figure 9b). TLS-estimated biomass, using DBH allometry according to Equation (2), was 7/% greater than the field-inventory estimate at 44.21 Mg ha−<sup>1</sup> of above-ground woody biomass and 22.10 Mg C ha<sup>−</sup>1.

Mapping out the distribution of individual tree biomass spatially showed that field records and TLS predictions were very closely co-aligned (Figure 11a–b). The TLS tree location mapping showed very few instances of omission from field-mapped locations, and a number of commissions which were predominantly located at the very boundary of the plot (Figure 11c). Omissions in the field-inventory mapping along the plot boundary (e.g., south-western corner) are indicative of positioning errors associated with the manual boundary line placement in the field.

**Figure 10.** Comparison of TLS-derived DBH with field-measured values for the dominant woody species.

#### *3.3. Allometric Scaling*

Given that the MLS and ULS point clouds from our acquisitions were not suitable for direct DBH/AGB retrieval (Figure 9), we explored the possibility of predicting AGB in the MLS and ULS datasets through allometric scaling of tree crown parameters. MLS and ULS captured the general crown structure of individual trees very well (Table 3), and the relationship between tree height and crown area that was sensed by the three scanning systems was very consistent (Figure 12). However, both multiple linear regression and Random Forest modeling failed to establish a reliable relationship between canopy attributes and stem DBH/AGB (R<sup>2</sup> = 0.34, RMSE = 0.075 m).

**Table 3.** Summary statistics of individual tree attributes collected by the three LiDAR systems.


**Figure 11.** Spatial distribution of tree biomass as mapped through field inventory (**a**, n = 572) and from TLS (**b**, n = 669) using DBH-based allometric equations. Agreement, omissions, and commissions between field and TLS mapping are shown in ( 93 **c**) with overlay blending.

**Figure 12.** Relationship between individual tree height and crown area, as determined from the three different LiDAR systems.

#### **4. Discussion**

Laser-based technologies have transformed our ability to measure and monitor savanna vegetation structure. The level of detail and efficiency exhibited in this study from all three sensors is very encouraging from a long-term monitoring and management perspective—providing rigorous baselines from which to assess current state and future dynamics.

#### *4.1. Efficient Monitoring of Habitat Structure*

TLS has emerged as the gold standard for 3D characterization of vegetation structure in ecosystems around the globe [17,19]. The high quality results obtained here from TLS are unsurprising, but it is worth noting that the scan position density that we employed (30 m spacing, 16 scans ha−1) was considerably lower than current recommendations in the literature which suggest aiming for 10 m spacing [37]. The TLS literature is heavily skewed towards forested ecosystems where high density scanning is critical for minimizing occlusion. In savanna woodlands and shrublands, which occupy 20% of the global terrestrial surface, the probability of occlusion is much lower and long-range scanners allow for wider scan position placement [20]. This is important in the context of habitat monitoring, as TLS is often considered too time intensive for broad extent and repeat coverage sampling.

Nonetheless, the time required for our TLS survey, and subsequent processing, was still substantially longer than that required for MLS and ULS captures (Table 1). The efficiency of MLS and ULS in terms of time are clear, but do they provide sufficient structural information for ecosystem monitoring? From a canopy height model perspective there is minimal difference among the three scanning approaches and systems. The TLS CHM could not be matched by MLS or ULS at the 0.25 m resolution, but was indistinguishable from MLS at 2 m resolution and ULS at 4 m resolution. For many ecological applications 2–4 m is more than sufficient for providing insight into how canopy growth is tracking

over time, and the extended geographic coverage may often outweigh the benefits of finely detailed local structural representation.

The minor under-estimation of woody canopy cover by ULS was consistent across the subplots of our study site. Similar findings have recently emerged from dry sclerophyll forests in Tasmania, where ULS underestimated canopy cover by a few percent in comparison to TLS at the plot scale [38]. Here we show that this under-estimation is a consistent trend across multiple height classes, indicating that the source of this bias is not only under-storey, but also small canopy features such as the edges of branches for example. MLS on the other-hand consistently over-estimated canopy cover in comparison to TLS. This is contrary to the few published studies exploring MLS characterization of woody canopy—which have mostly been conducted in forested ecosystems where tall dense canopy restricts access of the laser beam through to the canopy surface [21,39,40]. In the open savanna woodland of our study, ground-based TLS and MLS can penetrate canopies well with good fields of view. In this instance we attribute the over-estimation of canopy cover by MLS to the nosier point cloud that it produced in comparison to TLS, and these results could likely be improved through further post-processing of the MLS point cloud.

#### *4.2. Individual Trees and Biomass Scaling*

Given that MLS and ULS were able to characterize individual tree height and crown attributes very well, albeit with minor bias, we had anticipated that estimation of aboveground biomass would have been possible from these point clouds through allometric relationships with crown attributes. However, our TLS data showed that DBH-AGB did not scale in a predictable manner with tree height, crown diameter, area, and volume (R<sup>2</sup> = 0.36, RMSE = 0.075 m). This was initially surprising in light of how well AGB scales as a function of tree height and crown diameter in many ecoregions around the globe [41], but it is an important reminder of how local-scale variability is averaged out in regional and global syntheses. Savanna eucalypts have notoriously complex crown structures, exacerbated by the actions of fire, termites and cyclone damage [42]. As shown recently in a similar savanna setting in Litchfield National Park, savanna trees exhibit a wide range of crown size variability for a given DBH [6]. In addition to confounding the estimation of DBH and AGB from crown attributes, this variability also brings into question the reliability of using DBH allometry in the first instance for biomass estimation in these systems. The growing use of TLS for building Quantitative Structure Models (QSMs) of trees and accounting for biomass volumetrically has started revealing errors associated with DBH-based allometry in a variety of ecosystem types [43,44] and needs greater exploration in different savanna settings.

#### *4.3. Limitations and Future Directions*

TLS is a mature field, with well established protocols for data acquisition and processing [17]. TLS can be deployed in many different ecosystems with high confidence placed in data capture and structural representation, with efficient open-source tools for tree segmentation and structural modeling [45]. Our results obtained from MLS and ULS are encouraging and capture key elements of the savanna woody structure thoroughly. As expected both MLS and ULS did not detect some of the finer-scale ecosystem elements as well as TLS, but are more than sufficient for addressing many ecological questions.

A challenge that arises when comparing different platforms is that the end results are influenced heavily by the acquisition and pre- and post-processing parameters. As such, meaningful generalisations are difficult to make. TLS has the advantage of being a mature technology, with a large user-base and established collection and pre- and post-processing protocols. MLS and ULS are at earlier stages of development and uptake. Furthermore, it is a rapidly expanding industry and a wide variety of sensors and platform exist, from very high-end survey-grade systems (e.g., RIEGL RiCOPTER flying a RIEGL VUX-1) to lighter weight and lower cost options (e.g., DJI M600 with a Velodyne or M8 sensor). Some MLS systems rely solely on SLAM for positioning, while others have integrated RTK-GNSS, so results and noise levels will vary. Additionally, the actual flight/walk patterns and speeds have a marked impact on point-cloud characteristics—and optimal acquisition parameters may be ecosystem specific. Fortunately, since the point-cloud outputs from MLS and ULS share many features of TLS point clouds (at lower densities) analysis of these products can make use of and further advance breakthroughs in TLS processing pipelines.

A key challenge moving forward with MLS and ULS surveying will be to optimize the acquisition and processing parameters in a robust and repeatable manner. We suggest that TLS be leveraged at the start of baseline and long-term monitoring programs to provide a benchmark against which the MLS and/or ULS acquisitions can be optimized. This will enable MLS and ULS to be applied over larger areas and at more frequent time intervals, with intermittent re-calibration against nested TLS collections in key spatial locations and temporal stages throughout the monitoring program. Furthermore, establishing the scaling uncertainty of TLS measurements through to large area ULS surveying will greatly benefit the earth observation community by assisting in the calibration and validation of spaceborne optical, LiDAR and SAR missions (e.g., GEDI, ICESAT-2, NISAR, BIOMASS). When considering the faster acquisition time and the ability to collect data over larger areas, MLS and ULS systems could prove particularly useful for expanding the geographic range and representation of current calibration and validation libraries, with known uncertainties from the nested TLS models.

Lastly, while we have primarily focused here on comparing the individual point clouds and derived products produced through TLS, MLS and ULS surveying, we recognize that much could be gained through their integration. For example, a fused point cloud derived from both TLS and ULS could provide rich vegetation detail from the oblique TLS perspective and comprehensive terrain characterization from the aerial ULS viewpoint.

#### **5. Conclusions**

From a canopy modeling perspective, the differences between the TLS, MLS, and ULS were negligible for most applications, but there was a large difference in the acquisition time, with MLS and ULS offering distinct advantages. The trade-off comes through at the individual tree and stem modeling level, where TLS can reproduce DBH estimates from the field with high accuracy and precision, and capture smaller stems than are typically measured in field surveys. Comparisons among TLS, MLS and ULS are challenging as there are a large range of parameters that can be modified across a wide range of applications. TLS surveying is relatively mature, whereas optimal sampling and processing strategies for MLS and ULS are less mature, but rapidly developing. Future research should focus on using TLS as a benchmarking tool to optimize the acquisition parameters of MLS and ULS for a given research and/or management objective. Importantly, this needs to be done across a range of different ecosystem types as the strengths and weaknesses of each sensing system, and platform, will vary with site-specific conditions.

TLS provides a holistic representation of 3D structure that cannot be obtained through traditional field inventory and provides an avenue for optimizing MLS and ULS acquisition parameters, generating correction coefficients for MLS and ULS canopy products, and exploring new allometric models linking point-cloud metrics directly to above-ground woody biomass.

**Author Contributions:** Conceptualization, S.R.L. and R.B.; Data curation, T.W. and D.A.L.; Formal analysis, S.R.L. and T.W.; Methodology, S.R.L.,T.W., D.A.L. and M.R.; Project administration, R.B.; Visualization, S.R.L.; Writing—original draft, S.R.L.; Writing—review & editing, S.R.L., T.W., D.A.L., M.R. and R.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was not externally funded.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available on request due to geographic restrictions.

**Acknowledgments:** We thank Jon Schatz and Jaylen Nicholson for conducting the field inventory, with assistance from Pierre Grandclément. Linda Luck and Stephanie Johnson provided assistance with data processing. This manuscript benefited greatly from the comments of two anonymous reviewers.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

