Impact of Soil Reflectance Variation Correction on Woody Cover Estimation in Kruger National Park Using MODIS Data

Time-series of imagery acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) has previously been used to estimate woody and herbaceous vegetation cover in savannas. However, this is challenging due to the mixture of woody and herbaceous plant functional types with specific contributions to the phenological signal and variations in soil background reflectance signatures observed from satellite. These factors cause variations in the accuracy and precision of woody cover estimates from different modelling approaches and datasets. Here, woody cover is estimated over Kruger National Park (KNP) from the MODIS 16-day composite time-series data using dry season NDVI/SAVI images and applying NDVIsoil determination methods. The woody cover estimates when NDVIsoil was ignored had R2 = 0.40, p < 0.01, slope = 1.01, RMSE (root mean square error) = 15.26% and R2 = 0.32, p < 0.03, slope = 0.79, RMSE = 16.39% for NDVIpixel and SAVIpixel, respectively, when compared to field plot data of plant functional type fractional cover. The woody cover estimated from the soil determination methods had a slope closer to 1 for both NDVI and SAVI but also a slightly higher RMSE. For a soil-invariant method, RMSE = 19.04% and RMSE = 17.34% were observed for NDVI and SAVI respectively, while for a soil-variant method, RMSE = 18.28% and RMSE = 19.17% were found for NDVI and SAVI. The woody cover estimated from all models had a high correlation and significant relationship with LiDAR/SAR based estimates and a woody cover map produced by Bucini. Woody cover maps are required for vegetation succession monitoring, grazing impact assessment, climate change mitigation and adaptation research and dynamic vegetation model validation.


Introduction
Several vegetation indices such as Leaf Area Index (LAI), normalised difference vegetation index (NDVI) and others, as well as biophysical parameters such as canopy height and aboveground biomass (AGB) have previously been used for mapping woody plant abundance [1][2][3].This is because woody fractional cover is needed as an input to many ecological models for the assessment of fire, deforestation, degradation, urban expansion and hydrological dynamics [4][5][6][7][8][9][10][11].Changes in woody cover can have profound effects and unforeseen consequences for ecosystems functioning [12,13].Information on woody fractional cover provides immense benefit to ecological modelling and helps in understanding ecosystem function in savannas [14,15], but spatially explicit information on woody fractional cover is difficult to obtain [7, 16,17].A proper understanding of vegetation structure and phenological characteristics is a key requirement for measuring woody cover [12].
Satellite remote sensing provides an indirect measurement of vegetation indices (e.g., NDVI) and biophysical parameters (e.g., LAI, fAPAR, AGB and fractional vegetation cover [FVC]) [18][19][20].Studies on woody cover estimation have been published previously [7, [21][22][23], but the accuracy and precision of woody fractional cover estimates vary with modelling approach, used datasets and ecosystems [1,21,24].Although estimates of land cover types using remote sensing data have an associated error and uncertainty of unknown magnitude, the estimate of woody fractional cover is very challenging, particularly in savannas.This is because savanna vegetation is not discrete but rather a continuum of a varying mixture of tree and grass plant functional types (PFTs), which can show clumping or patchiness at certain scales.
The existing vegetation continuous fields (VCF) product from MODIS does not accurately measure woody fractional cover, especially for open forests [10,21].White et al. [25] assessed MODIS VCF products in the Western USA, a region spanning semi-arid deserts, sparse dry woodlands, and cool mesic upland forests by using two independent ground-based tree cover databases.They reported an underestimation of low tree cover and overestimation of dense canopies in the MODIS product.An overall root mean squared error (RMSE) of 24% and 31% was found when comparing the MODIS retrievals with ground-based Forest Inventory and Analysis (FIA) data (1176 field plots) and Southwest Regional GAP (SWReGAP) datasets (2778 field plots).RMSE indicated a more positive values at > 10% cover than 15% for FIA and 12% for SWReGAP.At canopy cover >60%, the error is high (49% for FIA and 44% for SWReGAP).
Ibrahim et al. [26] estimated the fractional cover of PFTs in Kruger National Park (KNP) from MODIS NDVI time series using harmonic analysis.In that study, harmonic analysis was used to decompose the time series signal into amplitude, cycles, and phase.The field plot estimates of tree cover showed a significant correlation with the amplitude (r = −0.59,p = 0.001), phase of the first harmonic term (r = −0.73,p = 0.0001) and the number of cycles of the second harmonic term (r = 0. 56, p = 0.002).The tree cover estimated from the phase of the strongest harmonic term showed a strong linear relationship with field-measured tree cover (R 2 = 0.55, p < 0.01, slope = 0.93, RMSE = 13.26%).Ibrahim et al. [26] emphasized the importance of the phase of the significant harmonic terms for tree cover estimation.In constrast, here the MODIS NDVI data are not decomposed using harmonic analysis, but instead the impact of using a soil database for accounting for the variation of soil reflectance in the NDVI calculation for woody cover estimation in KNP is assessed.We chose KNP because it is dominated by savannah ecosystems (e.g., Skukuza thickets, open trees, dense trees, and bush savanna).The plot data we collected from the KNP span different vegetation types, geological conditions, and soil types and cover a gradient of woody/herbaceous mixtures that is very distinct in terms of structure, type, density, and distribution [26][27][28][29][30].
Many studies have applied soil-variant and invariant approaches (with respect to spatial variability of soil backgrounds in the NDVI) to estimate fractional vegetation cover [31][32][33].Zeng et al. [34] combined the International Geosphere-Biosphere Program (IGBP) land cover classification with 1 km NOAA AVHRR NDVI data and employed the fifth percentile of the histogram of the maximum NDVI for the barren or sparsely vegetated category as the NDVI soil to calculate global FVC.Zeng's method assumed that NDVI soil is invariant.Wu et al. [35] combined the Harmonized World Soil Database (HWSD) and annual minimum NDVI to calculate NDVI soil for different soil types and then estimated global FVC.Wu's et al. [35] method considers spatial variability of soil background.The HWSD provides global soil types with a spatial resolution of 1 km.Ding et al. [36] investigated the influence of variations of NDVI soil derived from FVC estimation using these two approaches proposed by Wu et al. [35] and Zeng's et al. [34].However, Ding's [36] methods used 564 reflectance spectra of soils collected in Northeast China.Validation results indicated that this approach that considered the spatial variability of soil background yields better estimates of FVC than using a soil-invariant NDVI soil value.The accuracy increases from RMSE = 7% to 10% [36].The soil-variant method is more robust when ground measurements of soil reflectance are available than using the soil database to determine the NDVI soil , but it is unclear whether ignoring NDVI soil in the absence of field spectral measurements yields higher errors than considering it using a global soil database like HWSD.This is especially true for small areas with less soil variation and low NDVI (e.g., savanna) as the uncertainty in determining the NDVI soil is higher in areas of low vegetation [36].
In this study, we investigate different approaches for accounting for the spatial variation of soil spectral signatures in estimating woody cover in African savanna using MODIS data.Specifically, the study (1) evaluates the accuracy of woody cover estimation using different regression models without accounting for soil background; (2) investigates whether accounting for NDVI soil obtained from a soil database reduces the woody cover prediction error; and (3) compares the estimated woody cover maps with other existing satellite-derived woody cover products.

The Study Area
Kruger National Park (KNP) is South Africa's largest nature reserve covering approximately 2,000,000 hectares.It extends 380 km from north to south and 60 km from east to west (Figure 1).Its elevation ranges from 260 m to 839 m above mean sea level.The mean annual rainfall ranges from 440 mm in the north to 750 mm in the south.Figure 1 shows the mean annual rainfall for the year 2002 to 2015 for the three weather stations in the area.
Remote Sens. 2019, 11, 898 3 of 25 considered the spatial variability of soil background yields better estimates of FVC than using a soilinvariant NDVIsoil value.The accuracy increases from RMSE = 7% to 10% [36].The soil-variant method is more robust when ground measurements of soil reflectance are available than using the soil database to determine the NDVIsoil, but it is unclear whether ignoring NDVIsoil in the absence of field spectral measurements yields higher errors than considering it using a global soil database like HWSD.This is especially true for small areas with less soil variation and low NDVI (e.g., savanna) as the uncertainty in determining the NDVIsoil is higher in areas of low vegetation [36].
In this study, we investigate different approaches for accounting for the spatial variation of soil spectral signatures in estimating woody cover in African savanna using MODIS data.Specifically, the study (1) evaluates the accuracy of woody cover estimation using different regression models without accounting for soil background; (2) investigates whether accounting for NDVIsoil obtained from a soil database reduces the woody cover prediction error; and (3) compares the estimated woody cover maps with other existing satellite-derived woody cover products.

The Study Area
Kruger National Park (KNP) is South Africa's largest nature reserve covering approximately 2,000,000 hectares.It extends 380 km from north to south and 60 km from east to west (Figure 1).Its elevation ranges from 260 m to 839 m above mean sea level.The mean annual rainfall ranges from 440 mm in the north to 750 mm in the south.Figure 1 shows the mean annual rainfall for the year 2002 to 2015 for the three weather stations in the area.The area is divided into two main geological zones.The western part is situated on granite and the east on basaltic bedrocks [26][27][28][29][30]. Geological bedrock influences soil formation processes, and indirectly plant species and abundance distribution, and ecological processes [30].KNP has diverse vegetation types (more than 1900 plant species) including woody and herbaceous plant functional types [27].Figure 2 shows the soil types in KNP and the main soil types in the field plot areas based on the soil database by the International Union of Soil Sciences (IUSS), the United Nations Environment Programme, the FAO, and the International Soil Reference and Information Centre (ISRIC) [37].The area is divided into two main geological zones.The western part is situated on granite and the east on basaltic bedrocks [26][27][28][29][30]. Geological bedrock influences soil formation processes, and indirectly plant species and abundance distribution, and ecological processes [30].KNP has diverse vegetation types (more than 1900 plant species) including woody and herbaceous plant functional types [27].Figure 2 shows the soil types in KNP and the main soil types in the field plot areas based on the soil database by the International Union of Soil Sciences (IUSS), the United Nations Environment Programme, the FAO, and the International Soil Reference and Information Centre (ISRIC) [37].

Field Data
Percent woody, herbaceous and bare soil cover were estimated for 25 field plots in the KNP during a field campaign carried out in March 2015 at the end of growing season.An ocular method proposed by Law et al. [38] was adopted for the estimation of percent cover for four structural vegetation types and bare soil within the plots (trees > 6 m, tree and shrubs 1-6 m, forbs and grasses) located along the road from Skukuza to Tshokwane (Figure 1).Three additional plots were added based on visual interpretation of Google Earth images to incorporate areas with denser tree cover than what was sampled in the field.The Google Earth imagery was acquired in May 2015.However, The park has been classified into landscapes (e.g., Skukuza thickets, open trees, dense trees, and bush savanna) based on geomorphology, climate, soil and vegetation pattern [27][28][29].The major woody species in the southern part include Combretum apiculatum, Acacia nigrescens, Spirostachys africana, Combretum hereroense, Sclerocarya birrea, Terminalia sericea, Combretum zeyheri, while the drier northern part is dominated by Colophospermum mopane (mopane) savanna.

Field Data
Percent woody, herbaceous and bare soil cover were estimated for 25 field plots in the KNP during a field campaign carried out in March 2015 at the end of growing season.An ocular method proposed by Law et al. [38] was adopted for the estimation of percent cover for four structural vegetation types and bare soil within the plots (trees > 6 m, tree and shrubs 1-6 m, forbs and grasses) located along the road from Skukuza to Tshokwane (Figure 1).Three additional plots were added based on visual interpretation of Google Earth images to incorporate areas with denser tree cover than what was sampled in the field.The Google Earth imagery was acquired in May 2015.However, despite the two months difference between the field data and the Google Earth imagery, we assumed that woody phenology does not vary much within this period.Previous studies found that woody phenology is more stable over time than herbaceous phenology [28].The field survey has considered the MODIS satellite pixel size (250 m × 250 m) by a sampling plot almost equal to the size (200 m × 200 m) of the pixel.Our observations in the field helped us to understand that vegetation cover does not vary that much over that spatial distance.The uncertainty is considered minimal.The 28 plots were placed in different vegetation types, geological conditions and soil types.Some of the tree species found in these plots include Combretum zeyheri/apiculatum (e.g., plot 14, 10, 4), Dichrostachys cinerea, Sclerocarya birrea, Terminalia sericea (e.g., plot 28), Acacia gradicortuna (e.g., plots 17, 25, and 19), Combretum hereroense (e.g., plot 19) and Albizia harveyi (e.g., plot 9).A detailed description of the field data can be found in Ibrahim et al. [26].

MODIS NDVI Data
MODIS NDVI data (MOD13Q1) were obtained from the National Aeronautics and Space Administration (NASA) via https://modis.ornl.gov/cgi-bin/MODIS/global/subset.pl.MOD13Q1 is a gridded level 3 product provided at 250 m spatial resolution every 16 days produced from atmospherically corrected bi-directional surface reflectance factors (BRFs) and masked for water, clouds, and cloud shadows [39,40].A previous study developed a method based on MODIS Ross-Thick and Li-Sparse kernels to estimates of BRDF effects in NOAA-AVHRR NDVI time series.The results indicated that, in most cases, the uncorrected NDVI time series do not reflect actual seasonal and interannual variation in vegetation greenness.It was found out that the techniques reduce BRDF effects in AVHRR NDVI observations by about 50% to 85% [41].MODIS NDVI has been used widely for retrieving vegetation composition such as vegetation structure and annual net primary productivity (ANPP) dynamics in grassland-shrub land areas [42], tree cover change [2], tree-grass separation/green-up dates [28], for monitoring vegetation and land use dynamics [43,44] and for the analysis of trends to assess the CO 2 fertilization effect on vegetation [45].NDVI has also been used to examine the relationship between vegetation productivity and rainfall distribution along environmental gradients [43,44].Chamaille-Jammes and Fritz [44] investigated the relationships between precipitation and primary production along a precipitation gradient.The variability of both precipitation and primary production was measured as interannual coefficient of variations (CVs), which decreased with increasing mean annual precipitation (MAP) (respectively F(1, 31) = 21,88, p, 0.0001 and F(1, 31) = 12.28, p 0.0014).The CV of annual NDVI was positively correlated with the CV of precipitation (F (1, 31) = 35.480,p, 0.0001).The study affirmed the finding of study that the sensitivity of NDVI to precipitation decreases with increasing MAP.NDVI is used here as the proxy of vegetation cover as numerous studies have identified a strong relationship between the NDVI and NPP [46][47][48].Prince and Goward [46] designed a Global Production Efficiency Model based on the production efficiency concept.The model relies on variables (e.g., NDVI, temperature, etc.) that can be remotely sensed at global scale.It is mechanistic and does not require a correlation between NDVI and primary productivity.Moreno-de las Heras et al. [42] found the strong relationships between field net primary production (ANPP) and the annual integrals (per growing cycle) of herbaceous (R 2 = 0.67, p < 0.001) and shrub (R 2 = 0.65, p < 0.001) NDVI components.Zhu and Southworth [47] predicted NPP from the GIMMS3g NDVI data and correlated it with annual MODIS NPP stratified by savannah type.Highly significant linear correlations were found for tree savanna, bush savanna, and grassland savanna with correlation coefficients of 0.77, 0.74, and 0.80, respectively.

MODIS SAVI Data
The Soil Adjusted Vegetation Index (SAVI) was developed specifically to minimize soil brightness effects on vegetation indices derived from red and near-infrared (NIR) wavelengths [49].In this study, NDVI and SAVI vegetation indices were derived from MODIS images for the period 12 July to 29 August 2014 (before the start of season) and 9 May to 26 June 2015 after the end of season.This season is best for modelling woody cover using passive optical data in savannah [50,51] because the remote sensing data for the beginning of the dry season provide the best discrimination of green woody canopies from the senesced grass layer [50,51].The SAVI index was calculated as follows: where the NIR is the near-infrared band, R is the red band and L stands for a soil correction factor (ranges from 0 to 1).L = 0.5 was used in this study as it has been recommended for savanna ecosystems by previous study [52].

MODIS Vegetation Continuous Field (VCF)
MODIS VCF is an Earth observation product describing percent tree cover, percent non-tree vegetation (mostly herbaceous) and percent bare surface and is provided globally at 250 m resolution by NASA, from the Land Processes Distribution Active Archive Centre (LP DAAC) available at http://e4ftl01.cr.usgs.gov/MOLT/.It is called MOD44B as a standard MODIS product [53,54].The collection 5 of this product (version 0051) was used in this study, being the most recent version of this dataset at the time of writing the manuscript.Specifically, MODIS VCF for the years 2008 and 2014 were used in this study for validation and comparison purposes.

Validation Datasets
Two woody cover map of KNP were used as independent validation dataset.A 2010 woody cover map was produced using 14 dual-polarized (HV, HH) 12.5 m L-band ALOS PALSAR images trained with a random forest algorithm and 25,000 ha of airborne LiDAR data [55].The L-band Synthetic Aperture Radar (SAR) data was shown to produce higher accuracy woody cover products compared to C-and X-band data in southern African savannas [56].The LiDAR data were acquired in April 2008 (end of wet season) when woody plants were leaf-on, and the SAR images in July-August 2008 (dry season, leaf-off) to avoid soil moisture effects on the radar signal [56].This was shown to be the best season to model woody cover in the region with SAR data [6,57].Woody plants of at least 1 m canopy height were included.For details of the LiDAR and SAR datasets, see [56].Validation of the SAR-map with independent LiDAR data produced an R2 = 0.8 and RMSE = 7.7% [56].The 12.5 m LiDAR/SAR product was resampled to the 250 m MODIS resolution.In addition, the Bucini's woody cover map was provided by Scientific Services (GIS unit) of SANParks [50,51].The map was modelled through multiple regression techniques using a combination of optical Landsat ETM+ data, JERS-1 SAR scenes (L-band, HH polarization), and field plots of woody cover.The best predictive model was selected based on the Akaike information criterion [50,51], and the 90 m pixel size map represented woody cover conditions for the years around 2000.

Woody Cover Estimation
The data sets, approaches and the processing procedure implemented for the estimation of woody cover are summarized in Figure 3.The reference plot data from the field were used to extract the MODIS data (NDVI pixel and SAVI pixel ).The NDVI pixel and SAVI pixel and invariant and variant NDVI and SAVI generated from the original data were used for model calibration using regression models.Woody cover was then extrapolated using the regression equations.The woody cover estimates were validated using the half of field data on woody cover.We validated the LiDAR/SAR, Bucini, and MODIS VCF products using our field data.We compared our estimates of woody cover with the LiDAR/SAR, Bucini, and MODIS VCF products (Figure 2).All statistical analyses were carried out in R. Detailed explanation on methods is given in the following sections.

Woody Cover Estimate Using NDVIpixel and SAVIpixel
The mean NDVI and SAVI over the dry period were calculated and taken as woody (NDVIpixel and SAVIpixel).In this ecosystem, woody species have two growing cycles at the time when the herbaceous layer is dormant.Grass usually dries up before the woody species lose their leaves in autumn so that two short periods with dry grass and a green woody canopy exist (before and after the wet seasons).These periods occur before and after the wet season or in autumn and spring (April-May, and October-November, respectively).The wet season starts from October and ends in April.This is considered useful in capturing the phenology of woody plants.Some woody species fully green up before the first significant rains (e.g., Sclerocarya birrea, Acacia nigrescens).Woody species in KNP usually take eight weeks to reach full leaf flush from the woody leaf-out onset.However, some woody species such as Combretum apiculatum are usually late in their leaf flush but take shorter periods to develop full leaf cover than early leaf flushers [28].To reduce an overlap of woody and herbaceous phenology occurring probably due to a delayed start or end of season, the images chosen were for the months of July and August (12 July to 29 August 2014) before the start of season and May and June (09 May to 26 June 2015) after the end of season.
The assumption we make is that the influence of bare soil on the NDVI is minimal over a small area.The soil correction factor usually applied to derive SAVI reduces the influence of bare soil.Although both vegetation indices are sensitive to fractional vegetation cover, they are also sensitive to the soil background [4,36].While the SAVI approach is subject to some uncertainty in woody cover estimation, the effects of soil reflectance due the nature of soil type characteristics (e.g., such as soil

Woody Cover Estimate Using NDVI pixel and SAVI pixel
The mean NDVI and SAVI over the dry period were calculated and taken as woody (NDVI pixel and SAVI pixel ).In this ecosystem, woody species have two growing cycles at the time when the herbaceous layer is dormant.Grass usually dries up before the woody species lose their leaves in autumn so that two short periods with dry grass and a green woody canopy exist (before and after the wet seasons).These periods occur before and after the wet season or in autumn and spring (April-May, and October-November, respectively).The wet season starts from October and ends in April.This is considered useful in capturing the phenology of woody plants.Some woody species fully green up before the first significant rains (e.g., Sclerocarya birrea, Acacia nigrescens).Woody species in KNP usually take eight weeks to reach full leaf flush from the woody leaf-out onset.However, some woody species such as Combretum apiculatum are usually late in their leaf flush but take shorter periods to develop full leaf cover than early leaf flushers [28].To reduce an overlap of woody and herbaceous phenology occurring probably due to a delayed start or end of season, the images chosen were for the months of July and August (12 July to 29 August 2014) before the start of season and May and June (09 May to 26 June 2015) after the end of season.
The assumption we make is that the influence of bare soil on the NDVI is minimal over a small area.The soil correction factor usually applied to derive SAVI reduces the influence of bare soil.Although both vegetation indices are sensitive to fractional vegetation cover, they are also sensitive to the soil background [4,36].While the SAVI approach is subject to some uncertainty in woody cover estimation, the effects of soil reflectance due the nature of soil type characteristics (e.g., such as soil brightness, moisture) is likely to be negligible because of the small spatial scale being considered [58].In addition, the soil background reflectance values are lower than the canopy reflectance due to high albedo in the tropics [33].Jiang et al. [51] found out that the nonlinearity of NDVI over partially vegetated surfaces is more prominent with darker soil backgrounds and shadow [59].
The aim of this study is to establish the fundamental relationships between NDVI pixel, SAVI pixel and field data on woody cover to develop a calibration technique to assess the extent to which such relationships are able to estimate woody cover.Both linear and nonlinear regression analyses methods have been used for data calibration as the relationship between NDVI and measurements of canopy structures vary with vegetation types and seasonality [60].

Woody Cover Estimation Using NDVI soil Determination Methods
In this method, it is first assumed that each pixel consists of three fractional covers: woody cover (T), bare soil cover (S), and grass cover (G).In savannas, during the dry season, most of the non-woody fraction is occupied by bare soil or dry grass while, in the wet season, green grass fractional cover makes up most of the contribution [28,61].In the dry season, the grass layer becomes non-photosynthetic and dries up.The green grass layer (fraction of photosynthetically active vegetation) decreases from 85% to 8% while the fraction of non-photosynthetic vegetation of the same layer increases from 7% to 79% between the beginning and end of dry season [62].Furthermore, an investigation using field reflectance measurements of bare soil, grass and woody indicated that dry grass had the lowest NDVI [62].Hence, in this study, the non-photosynthetic grass layer was merged with bare soil.NDVI pixel was attributed to only woody and bare soil as expressed in Equation ( 2).The same was applied to SAVI pixel : where f_t is tree percent cover, f_s is fractional cover of bare soil, f_g is grass percent cover and NDVI pixel is the mean dry season NDVI from the MODIS data.This means that the influence of bare soil in the vegetation index does not usually allow spectral signals of vegetation to vary.To indicate the spatial variability of PFTs, the contribution of bare soil should be accounted for.Different techniques of vegetation fractional cover estimation have previously been proposed, which account for bare soil contribution.Some of these techniques are invariant to soil types and characteristics [31][32][33].These methods rely upon the assumption that pixels with FVC = 1 and 0 exist in an image.These are described by NDVI veg and NDVI soil for maximum vegetation and bare soil, respectively.Hence, FVC is calculated as [63]: Gutman and Ignatov [31] used low spatial resolution data (0.15 × 0.15) and estimated NDVI veg as 0.52 ± 0.03 and NDVI soil as 0.05 ± 0.03.Similarly, Sobrino and Raissouni [32] provided a threshold NDVI veg of 0.5 and NDVI soil of 0.02 [35].In this study, the methods by Zeng et al. [34] and Wu et al. [35] were adopted with slight modification due our study site being relatively small, the spatial resolution, and the fact that we lacked soil spectral reflectance field measurements [31].The first method is invariant to soil characteristics for determining NDVI soil .Zeng et al. [34] determined NDVI soil by utilizing the percentile of vegetation types using the IGBP land cover classification.They used the fifth percentile of the histogram of the maximum NDVI for the barren or sparsely vegetated category as the NDVI soil , which was 0.05, to estimate global FVC.Note, however, that only woody cover is estimated as opposed to Zeng et al. [34] whose aim was to assess statistically most likely FVC using spectra of soil collected from different datasets.Therefore, woody fractional cover is estimated as: The same was applied for SAVI as follows: The procedure requires that the histogram for each land cover is computed.Considering the size of the study area, the histogram for the whole image was computed.The minimum and maximum values of NDVI are 0.12 and 0.65.Since the maximum NDVI for this image is 0.65, a lower NDVI soil is suggested for barren and sparsely vegetated areas.In this case, NDVI veg and NDVI soil are approximated as 0.02 and 0.7 for bare soil and maximum vegetation, respectively.The 0.7 is the maximum vegetation for the whole KNP.The 0.02 is the threshold for NDVI soil.The SAVI soil is thresholded at 0.05 and for maximum vegetation at 1.32 (the maximum for the whole KNP).However, the contribution of the NDVI soil for each pixel has been determined from the pure pixels based on its fractional cover of bare soil.The estimation of NDVI soil can be performed with considerable accuracy if available soil reflectance data from in situ measurements exist for the major types of soil in a given study area [4,35,36].Unfortunately, it is challenging to acquire this information [4].Consequently, many previous studies have relied on the Harmonized World Soil Database and IGBP land cover to assign NDVI soil for each vegetation type, especially at regional scale [34,35].
Wu et al. [35] used the Harmonized World Soil Database (HWSD) Version 1.21 which was produced by the International Institute for Applied Systems Analysis (IIASA) and the Food and Agriculture Organization of the United Nations (FAO) to determine NDVI soil for each soil type.The HWSD does not cover all the major soil types in KNP due to missing data from the map especially for the region of our field data.Hence, a global soil and terrain database at a scale of 1:1 million developed by the International Union of Soil Sciences (IUSS), the United Nations Environment Programme, the FAO, and the International Soil Reference and Information Centre (ISRIC) [37] was used here.Based on Wu et al. [35], the NDVI soil for the three types of soil in our plot locations (which include Regosols, Luvisols, and Nitisols) were thresholded at 0.21, 0.24, and 0.32, respectively.This approach slightly differs from Wu et al. [35].First, a linear method is applied as opposed to Wu et al. [35].Given the small size of the study area, NDVI soil is considered for each soil type and for each plot, while a single value for maximum NDVI (NDVI veg ) was considered for all locations.For the NDVI, Regosols, Luvisols, and Nitisols have been thresholded at 0.015, 0.02, and 0.02, respectively, while, for SAVI, the thresholds were 0.04, 0.05, and 0.06.The three types of soils are described based on the World Reference for Soil Resources (2014) according to the International Soil Classification System of the Food and Agriculture Organization (FAO) of the United Nations [64].Regosols are characterised as the very weakly developed mineral soils in unconsolidated materials that do not have a mollic or umbric horizon.They are generally fine-grained material and are particularly common in semi/arid areas [64].Luvisols have a higher clay content in the subsoil than in the topsoil, as a result of pedogenetic processes (especially clay migration) leading to an argic subsoil horizon [64].Nitisols are deep, well-drained, red tropical soils with diffuse horizon boundaries and a subsurface horizon with at least 30 percent clay.Nitisols are some of the most productive red tropical soils [64].

Regression Analyses
In this study, different types of regression models were used to assess the relationship between the percent woody cover from the field campaign in 2015 and the independent variables (NDVI or SAVI vegetation index).Only 50% of the field data on woody cover (14 plots) was used for model calibration while holding the remaining 50% back for validation.The power of our study design was assessed by using the Power and Sample Size Calculation software provided by Dupont et al. [65].In this study, with reference to sample size (14 plots), we found statistical power following Dupont et al. [65] methods for sample size and power calculations.Our sample size results in a statistical power of 0.98.The type of regression analysis depends on the nature of the phenology metric being used.The procedure assumes that the NDVI-fractional vegetation cover relationship (or SAVI) is a function of vegetation type, the influence of understory and bare soil [60,[66][67][68].The regression model applied for each phenology metrics is explained below: (1) Since the assumption to use the mean of the images of chosen MODIS data (NDVI pixel and SAVI pixel ) does not preclude the presence of bare soil, different regression models were tested for estimating woody cover.Although it has previously been reported that the relationship between vegetation indices (especially NDVI) and percent cover depends largely on vegetation type [60], or may even have a strong linear relationship with sparse vegetation [66,67], it is not well known how the relationship with woody cover would be in KNP.Therefore, simple linear, polynomial and logarithmic models were tested for both vegetation indices to find the best fit for percent woody cover estimation.
(2) Only a simple linear regression was applied to NDVI soil and SAVI soil determination methods.

Assessment of Model Performance
The assessment of model performance for fractional woody cover from the MODIS NDVI/SAVI time series data uses the remaining field plot data (14 plots) not used for calibration.The LiDAR/SAR-based woody cover map and the MODIS VCF datasets were first compared to the field observations to quantify their accuracies.The validation of MODIS VCF with field data uses the MODIS VCF data for the year 2014 since the field campaign was in 2015.To assess model performance for woody cover estimated in this study, the coefficient of determination (R 2 ) of the regression model was used to measure the strength of the relationship between the predicted and the observed values.The predicted data for each model are taken as the independent variable and the observed as the dependent [69].In addition, RMSE was used to determine the goodness-of-fit.All MODIS NDVI/SAVI woody cover estimates were validated with field data are from year 2014/15.

Comparison of Woody Cover Estimates with LiDAR/SAR and Bucini Woody Cover Maps
The Pearson correlation coefficient is the measure of the bidirectional linear correlation between two variables and was used to assess whether woody cover estimates from this study as well as the MODIS VCF product are related to the LiDAR/SAR and Bucini woody cover maps.The significance of the relationship was also assessed at alpha = 0.05.

NDVI pixel and SAVI pixel for Woody Cover Estimate
Data on woody and herbaceous cover, and NDVIpixel and SAVIpixel used for woody cover estimation are shown in Table 1. Figure 4 shows the relationship between NDVI pixel and SAVI pixel and percent woody cover surveyed during the 2015 field campaign.The NDVI pixel versus the percent tree cover (Figure 4a,c,e) had moderate coefficient of determination with R 2 between 0.53-0.58and p < 0.01.The relationships for linear and polynomial regressions yielded the strongest fits (R 2 = 0.56, p < 0.01 and R 2 = 0.58, p < 0.01.Although the difference between linear and nonlinear regression (Figure 4e: R 2 = 0.53, p < 0.01) is not important, the nonlinearity of the NDVI increases with increasing species composition.In the dry season, the woody species are active while other PFTs (e.g., grasses) are dry making NDVI more sensitive to vegetation.The accuracy of SAVI pixel is higher (Figure 4b,d,f) than that of NDVI pixel with R 2 = 0.56 to 0.67 (p < 0.01).The saturation of NDVI pixel occurs at a higher percent cover than that of SAVI pixel .This is expected since SAVI applied correction factors, which minimise the effect of the soil background.Overall, the existing relationship between the NDVI pixel , SAVI pixel with percent woody cover implies that they can both be used to estimate percent woody cover.Table 2 shows the field plots data for the percent woody, herbaceous, bare soil and the type of soil for each calibration plot.Table 2 also indicates the fraction of the NDVI (NDVIZeng and NDVIWu) and SAVI (SAVIZeng and SAVIWu) estimated using the two soil background effect correction methods.

NDVI and SAVI for Woody Cover Estimation with NDVI soil Determination Using a Modified
Procedure by Zeng et al. [34] and Wu et al. [35] Table 2 shows the field plots data for the percent woody, herbaceous, bare soil and the type of soil for each calibration plot.Table 2 also indicates the fraction of the NDVI (NDVI Zeng and NDVI Wu ) and SAVI (SAVI Zeng and SAVI Wu ) estimated using the two soil background effect correction methods.Figure 5a-d shows the calibration results for NDVI and SAVI.NDVI estimates for both methods show an increased accuracy much better than when NDVI soil is not removed.The invariant method for which an NDVI soil threshold of 0.02 was used has a strong relationship with percent woody cover: R 2 = 0.67, p < 0.01 (Figure 5a) while the other approach which considered the Harmonized World Soil Database to determine the NDVI soil for each soil type in the plot locations also had a strong relationship with percent cover (R 2 = 0.67, p < 0.01, Figure 5c).There is a slight difference between the two vegetation indices as only three types of soil were found in the plot locations based on the soil database.
The invariant soil determining method for SAVI Zeng which threshold NDVI soil at 0.05, is not very effective as its accuracy (R 2 = 0.50, p < 0.01) is slightly lower than the initial relationship for which the soil contribution is unaccounted.This means that the invariant method applied here may be less accurate in inferring woody cover compared to other approaches though validation results might show otherwise, while, when considering the soil type in determining the NDVI soil , a strong relationship is observed between the SAVI Wu and percent woody cover (R 2 = 0.80, p < 0.01).In this case, the NDVI soil 0.04 was threshold for Regosols while Luvisols and Nitisols at 0.05 and 0.06, respectively.Overall, all vegetation indices have shown a good relationship with the percent woody cover.

MODIS NDVIpixel and SAVIpixel Woody Cover Maps
Table 3 shows field validation plots and NDVIpixel, SAVIpixel as well their corresponding woody cover estimated from simple linear, polynomial and logarithmic regression equations.Figure 7a,

MODIS NDVI pixel and SAVI pixel Woody Cover Maps
Table 3 shows field validation plots and NDVI pixel, SAVI pixel as well their corresponding woody cover estimated from simple linear, polynomial and logarithmic regression equations.Figure 7a,c presents an accuracy assessment of woody cover from the NDVI pixel and SAVI pixel (mean of dry season images for 2014/2015) using the field data from 2015.The estimated woody cover using linear regression has an R 2 = 0.40, p < 0.01, slope = 1.01 for NDVI pixel , RMSE = 15.26% and R 2 = 0.32, p < 0.03, slope = 0.79, RMSE = 16.39% for SAVI pixel .The level of accuracy for NDVI pixel and SAVI pixel with polynomial regression are not far from the simple linear regression (NDVI pixel : R 2 = 0.40, p < 0.01, slope = 0.89, RMSE = 15.21%;SAVI pixel : R 2 = 0.32, p < 0.03, slope = 0.78, RMSE = 15.39%).The logarithmic model is slightly less accurate for both vegetation indices (NDVI pixel : R 2 = 0.40, p < 0.01, slope = 0.79, RMSE = 15.44%;SAVI pixel : R 2 = 0.32, p < 0.03, slope = 0.82, RMSE = 16.51%).These results suggest that both NDVI pixel and SAVI pixel are sensitive to percent woody cover.However, with reference to the scatterplots (Figure 7a,c), the nonlinearity of the two relationships is evident.Table 4 shows field validation plots and their corresponding woody cover estimated using two 3.3.2.NDVI and SAVI Woody Cover Estimation from the NDVI soil Determination Methods by Zeng et al. [34] and Wu et al. [35] Table 4 shows field validation plots and their corresponding woody cover estimated using two approaches which account for NDVI soil and SAVI soil in the estimation.Figure 8a shows the validation of woody cover estimates from a modified procedure of vegetation fractional cover estimation by Zeng et al. [34] for NDVI (Figure 8a) (R 2 = 0.40, p < 0.01, slope = 1.06;RMSE = 19.04%)as well as for SAVI (Figure 8a) (R 2 = 0.32, p < 0.3, slope = 1.06;RMSE = 17.34%).The woody cover estimated for both vegetation indices using Zeng's procedure indicated that the approach can be used to infer woody fractional cover using dry season images even though the accuracy of the estimated woody cover were slightly lower than when NDVI soil were unaccounted for.There is an overestimation of woody cover in the lower percent cover where the contribution of soil is higher demonstrating the implication of an invariant NDVI soil removal approach.Figure 8b shows the validation of woody cover estimated from the modified procedure of vegetation fractional cover estimates by Wu et al. [35] (NDVI: R 2 = 0.40, p < 0.01, slope = 0.98; RMSE = 18.28%,SAVI: R 2 = 0.32, p < 0.02, slope = 0.88; RMSE = 19.17%).The accuracy of this approach is slightly higher than for the previous method.The difference between the two methods is more obvious in the RMSE and slope demonstrating the importance of the NDVI soil determination method that considers soil type characteristics.5 presents a correlation coefficient and RMSE of NDVI, SAVI, and MODIS VCF tree cover estimates with LiDAR/SAR and Bucini using 14 validation plots from the field campaign in 2015.All vegetation indices have a significant relationship with previous tree cover estimates except the polynomial model in NDVI pixel with LiDAR/SAR.The linear model had the highest correlation for both vegetation indices (NDVI pixel : r = 0.52, p = 0.05 with LiDAR/SAR and r = 0.63, p = 0.014 with Bucini; SAVI pixel : r = 0.53, p = 0.05 with LiDAR/SAR and r = 0.59, p = 0.02 with Bucini).However, the logarithmic model recorded the lowest RMSE in both NDVI pixel (RMSE = 15.99) and SAVI pixel (RMSE = 14.93) compared to linear (NDVI pixel : RMSE = 16.46 and SAVI pixel : RMSE = 15.15) and polynomial models (NDVI pixel ; RMSE = 17.14 and SAVI pixel : RMSE = 15.98).The correlation between MODIS VCF with the previous tree cover estimates is however not significant r = 0.39, p = 0.16, RMSE = 23.36 with LiDAR/SAR and r = 0.40, p = 0.17, RMSE = 40.83with Bucini) and had higher RMSE.The previous products on woody cover have been validated by the providers and were found relatively accurate.The LiDAR/SAR and Bucini woody cover map have been found the most accurate using our field data (Figure 6).Despite the differences between the times of field campaign data that were used for this study (2015), the LiDAR/SAR and Bucini [48] woody cover maps are consistent with our field measurements (Figure 6).Although LiDAR/SAR woody cover has an advantage since LiDAR has the ability to measure vegetation in three dimensions [29,70], another important consideration is acquiring the SAR images in July-August 2008 (dry season, leaf-off) to avoid soil moisture effects on the radar signal [6].The Bucini woody cover map was produced from the synergy between optical (Landsat ETM+ and JER-S) and SAR data.The accuracy of this product might result due to consideration of phenology of woody species (dry season images for the optical dataset to maximize discrimination of woody vegetation) as well as for accounting the effects of climate, soil characteristics, topography, fire frequency and herbivory in a regression analysis to estimate woody cover.
Despite reported strengths of MODIS VCF datasets observed in many studies [9,[70][71][72], the accuracy of the products is lower in savannas, particularly when certain statistical observations are put into consideration.The accuracy of MODIS VCF has the highest RMSE (28.6%) of all datasets compared here, probably due its model calibration, which only considered cover of trees taller than 5 m (Figure 6).This simply means that there is an underestimation of tree or woody cover from the MODIS VCF in this region where many trees are lower than 5 m.The underestimation of woody cover from the MODIS VCF compared to in situ data observed here is similar to the study by Brandt et al. [3], whose estimate of woody cover in the Sahel was nine times higher than the MODIS VCF.Consequently, low accuracy of MODIS VCF as a proxy for overall woody cover has been reported in the scientific literature [7,25,73].
The use of a large number of phenology metrics acquired in different periods regardless of vegetation dynamics [72,74], the presence of bad pixels (haze, cloud cover and shadow), the training datasets (regression tree usually require large samples), and limitations inherent in the MODIS viewing geometry (the effects is more with the individual bands than the vegetation index-NDVI) may be responsible for the limitations of the MODIS VCF in savannas.None of the MODIS VCF (2014) pixel values was of bad quality within the field plots used for this study.Consideration of the seasonal vegetation dynamics is useful for global scale mapping of woody cover in savannas using space observations.This is due to large differences in vegetation phenology during the wet and dry seasons [75] and consequent limitations such as cloud cover and sensor viewing geometry, which may affect the interannual and seasonal variation of woody/herbaceous phenology [41,76].

The NDVI pixel and SAVI pixel
The composite dry season images used for vegetation indices to estimate woody cover in KNP were accurate (Figure 7a-c).The validation of woody cover estimated from the NDVI showed a lower error (RMSE = 15.26%,15.21%, 15.44% for the linear, polynomial and logarithmic model respectively) for all models compared to SAVI pixel (RMSE = 16.39%,16.39%, 16.51% for linear, polynomial and logarithmic model, respectively) (Figure 7a-c), illustrating the strong dependence of NDVI on woody canopy structure during the dry season in KNP.This also demonstrated the presence of a photosynthetically active woody layer in the dry season as previously observed [50].The relationship between the NDVI and PFT types depends on the nature of the ecosystem in question [33,41,60].
The relationship between NDVI pixel / SAVI pixel and percent woody cover is approximately linear (Figure 7a).Most of the plots used for woody cover in this study are within the granitic zone, where soil types differ from the basaltic zone.The results discussed in this section agree with previous studies [60,66,67].A recent study found a linear relationship between FVC and the EVI as well as the SAVI vegetation index.In the same study, the relationship between the NDVI and FVC was nonlinear due to saturation effects at high vegetation fractions, the presence of shadow as well as the influence of soil background [77].
There is strong correlation between SAVI pixel / NDVI pixel and Bucini / LiDAR/SAR woody cover (Table 5).At the lower percent woody cover, the points are clustered.This demonstrated the influence of radiative transfer from the surface on canopy reflectance especially where there is mixed woody, herbaceous and bare soil fractions [70,[78][79][80].Therefore, in a sparse vegetation cover, radiative interactions cause soil to be prominent and canopy visible reflectance will contain a strong backscatter component.The soil reflectance effect is less prominent in the NIR canopy reflectance, since multiple scattering of NIR radiation by vegetative components dominate.Hence, the influences of soil reflectance reduces with decreasing canopy gaps [36,81].

Woody Cover Estimated Using Two NDVI soil and SAVI soil Determination Methods
The assessment of woody cover maps from the soil determination methods indicated a moderate linear relationship between the predicted and observed percent woody cover (Figure 8a,c).Although the slope of the regression line for both vegetation indices were higher (slope ranges from 0.88 to 1) (Figure 8a,c) compared to woody cover estimates without removal of the soil signal (slope ranges from 0.78 to 1) (Figure 8a), the RMSE is still high with soil determination methods.The RMSE for woody cover estimates without soil signal removal ranges between 15.21-16.51%.In contrast, the RMSE for the invariant method is 19.04% and 17.34% for the NDVI and SAVI vegetation indices estimates, respectively, while the RMSE for the variant method is 18.28% and 19.17% for the NDVI and SAVI estimates, respectively.Therefore, RMSE increases to more than 3% for the NDVI and about 5% for SAVI.The RMSE for SAVI vegetation index increases to about 2% with an invariant method and 4% with variant method.The SAVI vegetation index was found to be less sensitive to soil signal removal than the NDVI.SAVI is one of the vegetation indices specifically developed to reduce soil backgrounds effects.Although soil colour is useful for differentiating soil reflectance [80], soil moisture was considered an important factor in influencing vegetation indices [82].From the results in Figure 8, uncertainty in the estimates of percent woody cover at lower cover is high in all methods.The second approach of the soil determination methods overestimated low percent woody cover values.This can be explained by the sensitivity of the canopy NDVI or SAVI to soil background or the result of changing canopy structure, which might decrease in the NIR reflectance and increasing visible reflectance consequently leading to reduced NDVI.In addition, the sensitivity of vegetation indices to soil backgrounds was found to be greatest in the intermediate level of vegetation cover [36,83].It is usual to overestimate percent woody cover whenever soil contribution is underestimated and vice versa [4].Overall, despite the uncertainty in the estimation of percent woody cover with the soil determination approaches used in this study, heterogeneity is present in the woody cover estimates derived from these approaches.

The Uncertainties and Sources of Errors and Proposed Improvements
While our results demonstrated the potential of MODIS data to estimate percent woody cover from vegetation indices, woody cover estimated in this study has some limitations and remaining uncertainties that must be considered: (1) Phenology Phenology of PFTs in savannas is usually influenced by many environmental factors [84] Specifically, woody phenology is influenced mainly by temperature and day length [85], or precipitation and disturbance in certain conditions [86].For these reasons, the estimates of percent woody cover from passive optical remote sensing are less accurate compared to high-resolution datasets.While high resolution data such as LiDAR (active optical data) can determine woody canopy cover by measuring its 3D structure, the estimates from the passive optical imaging datasets mostly rely on the green canopy cover within a pixel [3].
The grass layer (fraction of photosynthetic vegetation) in savannas may have changed from 85% to 8% and the fraction of non-photosynthetic vegetation of the same layer may increase from 7% to 79% between the beginning and end of dry season [62].Although careful attention was given to defining the dry season, woody/herbaceous cover separation can be affected by certain grass species that are supported by soil moisture and temperature [28].In savannas, certain environmental factors favor grass growth and influence its phenology, productivity and biomass allocation [87,88].Temperature and soil moisture are generally not limiting factors of grass productivity in KNP [28].This might contribute to the overestimation of woody cover, especially in a highly mixed woody/herbaceous area.In addition, the time differences between the LiDAR/SAR, Bucini, and field data campaign may affect the accuracy assessment of the two products from MODIS compared here.Although the MODIS vegetation indices (e.g., NDVI) are less sensitive to illumination than individual bands, as previously reported in the literature, the estimates of woody cover in this study may well contain remaining uncertainties despite being specific to a particular season [41,76].The temporal granularity of the MODIS NDVI products of 16 days limits the temporal precision of the detection of changes in greenness.Denser time-series data would be preferable.
(2) The ground data (field plot data on percent woody/herbaceous cover) The calibration data used in our models may not be the representative of all species over the KNP landscape.The field method for woody cover estimation is a visual approach which may also constrain the accuracy of our models due to remaining uncertainties in the field data collection.However, the results presented in this study demonstrated that percent woody cover can be estimated from vegetation indices in savannas, and that a single regression model based on our field data indicated RMSE between 15-21%.However, while the accuracy assessment indicates lower errors, percent woody cover of <40%, where spectral signatures are probably dominated by understorey (dry grass) and soil, is not highly correlated with the field measurements.This might arise due to the presence of soil and dry grass, underestimation of percent woody cover in the field campaign or changes in woody phenology due to fire over the period [89].For both calibration and validation plots, there are only a few plots with percent woody cover >40%.This means that, if the regression models are to be developed and applied consistently over the large area, it is important that they are established on a much larger sample than as presented in the current study.This may reduce uncertainty and increase model accuracy.
(3) NDVI soil estimation The use of in situ measurements of soil reflectance remains a crucial step for an effective determination of NDVI soil in a pixel to estimate woody cover fraction [4,82,90], especially in savannas where the availability of vegetation indices at MODIS resolution of 250 m is essential for capturing vegetation and bare soil.One of the greatest challenges for woody cover estimation is the lack of soil reflectance since soil reflectance varies with soil type and characteristics (e.g., soil moisture).
In this study, as demonstrated by the validation results using different models, the influence of NDVI soil is minimal.Though the estimates of woody cover from the linear regression using soil determining methods had a slope closer to 1 (Figure 8a,b) for both vegetation indices, the NDVI pixel and SAVI pixel had the lowest RMSE (Figure 7a).This might be due to smaller size of the field plot data as well as difficulties in matching the soil type characteristics for both validation and calibration plots.Woody cover estimation with these approaches might introduce errors in a situation where only woody cover is being estimated.Challenges in woody/herbaceous or soil separation remain critical to model accuracy.However, the methods would have been more accurate if a larger environment is considered as some of these approaches are insensitive to a particular land cover type [34].The influence of spectral response pattern of both vegetation and soil can have strong temporal and spatial effects [58].The spatial effects may be negligible if a small area is considered.The temporal effects for soil [91] and vegetation are important as the species keep changing throughout the growing period, coupled with sensor limitations [33,41,92].
Smallman et al. [88] evaluated the role of repeated woody biomass estimates in constraining the dynamics of the major ecosystem carbon pools.They highlighted the challenges with dead organic carbon stocks and soil using the Harmonised World Soil Database (HWSD) to account for bare soil.In their estimates of carbon stocks, the in situ soil carbon observations have a lower uncertainty than those based on the HWSD.They stressed that the impact of the HWSD prior is reduced due to lack of a robust assessment of the uncertainty associated with the database and the lack of information on the time for which the priors are representative, necessitating a conservative use of the database.In situ measurements of soil reflectance (if available) would be more useful regardless of spatial scale [36,82,91].

Conclusions
A remote sensing-based model of woody cover retrieval in African savanna was developed from vegetation index metrics based on NDVI and SAVI derived from MODIS data and field data from 28 sites in KNP.The models were developed on the understanding that during the dry season only woody species are photosynthetically active.A strong linear relationship was found between the phenology and woody cover observations from a field campaign in 2015.The accuracy of the estimated woody cover had R 2 = 0.40, p < 0.01, slope = 1.01,RMSE = 15.26% and R 2 = 0.32, p < 0.03, slope = 0.79, RMSE = 16.39% for NDVI pixel and SAVI pixel , respectively.The percent woody cover estimated from the soil determination methods had an improved slope for both NDVI and SAVI but a slightly higher RMSE.Although it was not the primary objective of this study, it turns out that the LiDAR/SAR estimate is more accurate for woody cover estimates (R 2 = 0.45, p < 0.001, Slope = 0.5, RMSE = 15.90%)than other products tested in this study.The maps of woody cover will be useful in understanding woody/grass interactions in wooded savannas.Future work will have to ascertain the transferability of the method to savanna sites worldwide.

Figure 1 .
Figure 1.Location of the study area of Kruger National Park (KNP) (Southern part) in South Africa and its main river courses, indicating the locations of the sample plots of the field data collection in 2015.Blue circles indicate the field plot locations.The red points indicate the mean annual rainfall for the year 2002 to 2015 (14 years) of three stations in KNP.

Figure 1 .
Figure 1.Location of the study area of Kruger National Park (KNP) (Southern part) in South Africa and its main river courses, indicating the locations of the sample plots of the field data collection in 2015.Blue circles indicate the field plot locations.The red points indicate the mean annual rainfall for the year 2002 to 2015 (14 years) of three stations in KNP.

Figure 2 .
Figure 2. Soil types in KNP, the main soil types in the field plot area, field plot data (the blue and pink points indicate calibration and validation plots, respectively.

Figure 2 .
Figure 2. Soil types in KNP, the main soil types in the field plot area, field plot data (the blue and pink points indicate calibration and validation plots, respectively.

Figure 3 .
Figure 3. Flowchart showing data sets and processing procedure for the estimates of woody cover in KNP.

Figure 3 .
Figure 3. Flowchart showing data sets and processing procedure for the estimates of woody cover in KNP.

Figure 4 .
Figure 4.The relationship between NDVIpixel, SAVIpixel and field percent cover estimates with regression analyses, simple linear (a and b), polynomial (c and d) and nonlinear regression (e and f).

3. 1 . 2 .
NDVI and SAVI for Woody Cover Estimation with NDVIsoil Determination Using a Modified Procedure by Zeng et al.[34] and Wu et al.[35]

Figure 4 .
Figure 4.The relationship between NDVIpixel, SAVIpixel and field percent cover estimates with regression analyses, simple linear (a,b), polynomial (c,d) and nonlinear regression (e,f).

Figure 5 .
Figure 5. Calibration of NDVI and SAVI for woody estimate considering two methods for correcting soil variation background effects using a modified procedure by Zeng et al. [34] and Wu et al. [35]; NDVIZeng (a) and NDVIWu (c), NDVIZeng (b) and NDVIWu (d).

Figure 5 .
Figure 5. Calibration of NDVI and SAVI for woody estimate considering two methods for correcting soil variation background effects using a modified procedure by Zeng et al. [34] and Wu et al. [35]; NDVI Zeng (a) and NDVI Wu (c), NDVI Zeng (b) and NDVI Wu (d).

Figure 6 .
Figure 6.Validation of woody cover estimated from MODIS VCF, LiDAR/SAR, Bucini woody cover estimates (2001) observed woody cover from the field plot data 2015.

Figure 6 .
Figure 6.Validation of woody cover estimated from MODIS VCF, LiDAR/SAR, Bucini woody cover estimates (2001) observed woody cover from the field plot data 2015.

Figure 7 .
Figure 7. Validation of woody cover estimates derived with NDVIpixel and SAVIpixel using regression analyses, (a) with simple linear; (b) polynomial; (c) logarithmic regressions.

Figure 7 .
Figure 7. Validation of woody cover estimates derived with NDVI pixel and SAVI pixel using regression analyses, (a) with simple linear; (b) polynomial; (c) logarithmic regressions.

Figure 8 .
Figure 8. Validation of NDVI and SAVI woody cover estimates from the (a) modified procedure by Zeng et al. [34] and the (b) modified procedure by Wu et al. [35].3.3.3.Comparison of Estimated Woody Cover with LiDAR/SAR-Derived and Bucini's Woody Cover Estimates Using Pearson Correlation Analysis and RMSE Table 5 presents a correlation coefficient and RMSE of NDVI, SAVI, and MODIS VCF tree cover estimates with LiDAR/SAR and Bucini using 14 validation plots from the field campaign in 2015.All

Figure 8 .
Figure 8. Validation of NDVI and SAVI woody cover estimates from the (a) modified procedure by Zeng et al. [34] and the (b) modified procedure by Wu et al. [35].

Table 1 .
Woody, grass, bare soil, NDVI pixel and SAVI pixel for woody cover estimation.

Table 2 .
[35]mates of percent cover using NDVI and SAVI with soil determining methods using a modified procedure by Zeng et al.[34]and Wu et al.[35]

Table 3 .
Validation of MODIS NDVI pixel and SAVI pixel woody cover estimates.

Table 4 .
Validation of MODIS NDVI woody cover maps estimated using two methods of soil signal removal.

Table 4 .
Validation of MODIS NDVI woody cover maps estimated using two methods of soil signal removal.
3.3.3.Comparison of Estimated Woody Cover with LiDAR/SAR-Derived and Bucini's Woody Cover Estimates Using Pearson Correlation Analysis and RMSE

Table 5 .
Correlation and RMSE of estimated woody cover with LiDAR/SAR and Bucini estimates.