Optimisation of Savannah Land Cover Characterisation with Optical and SAR Data

: Accurately mapping savannah land cover at the regional scale can provide useful input to policy decision making efforts regarding, for example, bush control or overgrazing, as well as to global carbon emissions models. Recent attempts have employed Earth observation data, either from optical or radar sensors, and most commonly from the dry season when the spectral difference between woody vegetation, crops and grasses is maximised. By far the most common practice has been the use of Landsat optical bands, but some studies have also used vegetation indices or SAR data. However, conﬂicting reports with regards to the effectiveness of the different approaches have emerged, leaving the respective land cover mapping community with unclear methodological pathways to follow. We address this issue by employing Landsat and Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar (ALOS PALSAR) data to assess the accuracy of mapping the main savannah land cover types of woody vegetation, grassland, cropland and non-vegetated land. The study area is in southern Africa, covering approximately 44,000 km 2 . We test the performance of 15 different models comprised of combinations of optical and radar data from the dry and wet seasons. Our results show that a number of models perform well and very similarly. The highest overall accuracy is achieved by the model that incorporates both optical and synthetic-aperture radar (SAR) data from both dry and wet seasons with an overall accuracy of 91.1% ( ± 1.7%): this is almost a 10% improvement from using only the dry season Landsat data (81.7 ± 2.3%). The SAR-only models were capable of mapping woody cover effectively, achieving similar or lower omission and commission errors than the optical models, but other classes were detected with lower accuracies. Our main conclusion is that the combination of metrics from different sensors and seasons improves results and should be the preferred methodological pathway for accurate savannah land cover mapping, especially now with the availability of Sentinel-1 and Sentinel-2 data. Our ﬁndings can provide much needed assistance to land cover monitoring efforts to savannahs in general, and in particular to southern African savannahs, where a number of land cover change processes have been related with the observed land degradation in the region.


Introduction
Savannahs are important ecosystems that are found on almost half of the African continent and a fifth of the Earth's surface [1]. They consist of different densities of grasses and woody vegetation and are host to a number of ongoing land system science debates, such as equilibrium dynamics [2], the role of anthropogenic processes in land degradation [3] and the contribution of drylands to the • Do Landsat-based spectral metrics or L-band PALSAR data map the main savannah land cover types more efficiently? • Does the integration of Landsat and PALSAR data improve mapping? • Does single or multi-seasonal data produce the most accurate land cover models?

Study Area and Description of Main Land Cover Types
The area, comprising of 44,467 km 2 , was chosen in a north-south transect in a way that would allow it to be large enough to include a representative variation of the main land cover types found in the region. The choice was also dictated by the availability of the radar data. About 70% of the study area is in the Northwest Province of South Africa. It also extends to the Northern Cape Province (~10%), the Free State Province (~10%), and to Botswana (~10%; Figure 1). Temperatures range from 17 • C to 31 • C in the summer and from 3 • C to 21 • C in the winter. The annual rainfall is~400 mm, with nearly all of it falling during the summer wet months, between October and April [36]. Around 80% of the area falls within the Savannah Biome (Bushveld vegetation). The remainder falls within the Grassland Biome, which contains a variety of grasses typical of arid regions [37]. Within the South African part of the study area, 16 different vegetation types are found, mostly belonging to the thornveld, bushveld or savannah grassland categories [37]. In addition, a large part of the study area is arable land (e.g., the Vaalharts Irrigation Scheme [38]), contributing to a significant portion of the country's maize, groundnuts, sunflowers, dry beans and grain sorghum [39]. The area is also a vegetable and citrus fruit producer [40].
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 18 in the region. The choice was also dictated by the availability of the radar data. About 70% of the study area is in the Northwest Province of South Africa. It also extends to the Northern Cape Province (~10%), the Free State Province (~10%), and to Botswana (~10%; Figure 1). Temperatures range from 17 °C to 31 °C in the summer and from 3 °C to 21 °C in the winter. The annual rainfall is ~400 mm, with nearly all of it falling during the summer wet months, between October and April [36]. Around 80% of the area falls within the Savannah Biome (Bushveld vegetation). The remainder falls within the Grassland Biome, which contains a variety of grasses typical of arid regions [37]. Within the South African part of the study area, 16 different vegetation types are found, mostly belonging to the thornveld, bushveld or savannah grassland categories [37]. In addition, a large part of the study area is arable land (e.g., the Vaalharts Irrigation Scheme [38]), contributing to a significant portion of the country's maize, groundnuts, sunflowers, dry beans and grain sorghum [39]. The area is also a vegetable and citrus fruit producer [40]. We followed the Land Cover Classification System (LCCS) nomenclature developed by the Food and Agriculture Organisation of the United Nations (FAO, UN [41]) to identify the main land cover types of the study area. This is the same system employed by the South African National Geospatial Information (NGI) [42] to develop the South African national product of 2009-2013 [43], which we introduce in our analysis as a comparison dataset (in the absence of in-situ ground truth data). The LCCS class labelling approach describes the four main land cover types of the study area as follows (CSIR 2010): i Woody vegetation: trees, shrubs and bushes consisting of areas with either tree cover densities of >15% and the remaining woody cover, consisting of shrubs and bushes, no more than twice that of the tree cover, or of areas dominated by shrubs and bushes with a woody cover of >15% and a tree cover less than twice that of the remaining woody cover. Shrubs and bushes may vary in height up to 3 m. Trees are considered as woody vegetation with a crown elevation of >1.5 m above ground and a total height of >3 m. We followed the Land Cover Classification System (LCCS) nomenclature developed by the Food and Agriculture Organisation of the United Nations (FAO, UN [41]) to identify the main land cover types of the study area. This is the same system employed by the South African National Geospatial Information (NGI) [42] to develop the South African national product of 2009-2013 [43], which we introduce in our analysis as a comparison dataset (in the absence of in-situ ground truth data). The LCCS class labelling approach describes the four main land cover types of the study area as follows (CSIR 2010): i Woody vegetation: trees, shrubs and bushes consisting of areas with either tree cover densities of >15% and the remaining woody cover, consisting of shrubs and bushes, no more than twice that of the tree cover, or of areas dominated by shrubs and bushes with a woody cover of >15% and a tree cover less than twice that of the remaining woody cover. Shrubs and bushes may vary in height up to 3 m. Trees are considered as woody vegetation with a crown elevation of >1.5 m above ground and a total height of >3 m. ii Grassland: areas dominated by grasses with >4% vegetation cover. Areas may contain up to 15% woody cover. Grasslands may be almost bare during the dry season and during drought episodes. iii iii Cropland: areas where natural vegetation has been removed or modified and replaced by other types of vegetation cover of anthropogenic origin. This vegetation is artificial and requires human activities to maintain it in the long term. All vegetation that is planted or cultivated with an intent to harvest is included in this class, including cultivated herbaceous graminoids such as maize, sugar cane and cereals. It also includes small fields with subsistence crop farming, all herbaceous non-graminoids such as cotton, sunflower, potatoes, etc., pulses and orchards. The fields may be fallow at certain times during the year. iv iv Non-vegetated (or bare): this class consists of artificial cover as a result of human activities as well as areas with <4% vegetative cover, including bare rock, sands and unconsolidated bare soil such as animal feed lots, visible erosion, fine rock and soil fragments, as well as areas with bare soil resulting from unfavourable conditions (e.g., low rainfall).
Examples of these can be seen in the photographs of Figure 2.
Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 18 ii Grassland: areas dominated by grasses with >4% vegetation cover. Areas may contain up to 15% woody cover. Grasslands may be almost bare during the dry season and during drought episodes. iii Cropland: areas where natural vegetation has been removed or modified and replaced by other types of vegetation cover of anthropogenic origin. This vegetation is artificial and requires human activities to maintain it in the long term. All vegetation that is planted or cultivated with an intent to harvest is included in this class, including cultivated herbaceous graminoids such as maize, sugar cane and cereals. It also includes small fields with subsistence crop farming, all herbaceous non-graminoids such as cotton, sunflower, potatoes, etc., pulses and orchards. The fields may be fallow at certain times during the year. iv Non-vegetated (or bare): this class consists of artificial cover as a result of human activities as well as areas with <4% vegetative cover, including bare rock, sands and unconsolidated bare soil such as animal feed lots, visible erosion, fine rock and soil fragments, as well as areas with bare soil resulting from unfavourable conditions (e.g., low rainfall).
Examples of these can be seen in the photographs of Figure 2.

Aerial photographs
The sampling was carried out using colour aerial photography of 0.5 m pixels, freely-available from the NGI [42]. The NGI aims to capture 40% of the country every three years and the remaining areas every five years. We used the ESRI ArcGIS 10.5© environment for collecting the samples as the software makes the NGI aerial photographs available as a mosaicked product for the entire country

Aerial photographs
The sampling was carried out using colour aerial photography of 0.5 m pixels, freely-available from the NGI [42]. The NGI aims to capture 40% of the country every three years and the remaining areas every five years. We used the ESRI ArcGIS 10.5© environment for collecting the samples as the software makes the NGI aerial photographs available as a mosaicked product for the entire country [44]. Our study area is covered by photographs taken mainly in 2010 and, to a lesser extent, in 2008 (<5%) and in 2013 (<10%).

Landsat
The choice of Landsat data was driven by the need to coincide with the availability of the reference data (which cover the study area with aerial photographs taken between 2008 and 2013), and with the PALSAR data which were only available for the year 2008. All Climate Data Record (CDR) images with less than 80% cloud cover were obtained from the USGS EROS Data Center for the six WRS-2 tiles covering the study area ( Figure 1; Table 1). The images were preprocessed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) procedure, and clouds and cloud shadows were removed using F-mask [45,46]. Finally, the Normalised Difference Vegetation Index (NDVI [47]) was calculated. From the resulting 7-band image stacks (the six 30 m-pixel bands, plus the NDVI), spectral-temporal variability metrics were calculated, i.e., statistical values were derived from all co-located observations offering an effective method for generating meaningful land cover predictors [32,48]. We calculated five statistics for each of the seven bands: the mean, the standard deviation and three percentiles (25th, 50th and 75th), totaling 35 layers per pixel per season. To generate seasonal composites, metrics were calculated over two non-overlapping periods: the dry and the wet season. Moreover, in order to have sufficient data to allow for the calculation of the metrics, and to compensate for cloud contamination and data availability issues, imagery from five consecutive seasons were considered (Table 1). To choose the specific temporal windows for the dry and the wet seasons we consulted with the daily and dekadal rainfall estimates of the University of Reading TAMSAT group (version 3: https://www.tamsat.org.uk/data/rfe/; [49]). According to the TAMSAT data, in 2009 and 2010, the rains arrived after late October, and the green-up therefore started in November. Equally, in some years (2008-2010), the wet season lasted until mid-May, delaying the start of the dry season (browning). We therefore selected the period from 2 June to 1 October for the dry season, and from the 21 November to 31 March for the wet season. Figure 3 shows the Landsat 5 and 7 observations available per pixel for the study area for both the dry and wet seasons. It also includes an example, from 2008, of the data scarcity for the region, even for the dry season when cloud contamination was not an issue: large parts of the study area, especially in the south, are covered by less than three observations, which renders the calculation of metrics irrelevant. On average, about 30 non-cloudy observations per pixel were found within each 5-year wet and dry period (Table S1).

PALSAR
ALOS PALSAR-1 was a JAXA-managed L-band radar mission, operating from 2006 to 2011 at a wavelength of 23.6 cm. Data were collected across four polarizations: fine-beam single (FBS: HH or HV; spatial resolution = 10 m); fine-beam dual (FBD: HH + HV or VV + VH; spatial resolution = 20 m), or quad full polarimetry (HH + HV + VH + VV; spatial resolution = 30 m). ALOS operated a global image acquisition strategy; however, due to seasonal requirements for different regions, not all areas were imaged consistently.
Radar images were downloaded from the Alaska Satellite Facility (https://www.asf.alaska.edu/). For the dry season, the only data available within the 5-year window (that the Landsat metrics were calculated for), was the year 2008. Dual polarisation FBD images were acquired (HH, HV). For the wet season, only FBS (HH) bands were available for 2006-2007 and 2007-2008, and the latter was chosen. To ensure consistency across the study area, images were selected with the smallest possible time difference: a maximum of 17 days for both seasons (Table 2). Following download, the raw images were first converted to backscatter using the standard conversion equation: where σ 0 is backscatter, DN is the raw digital number and CF is a calibration factor (−83). Secondly, a Lee filter with a 7 × 7 window size was applied to minimise speckle effects [50], using Erdas IMAGINE 2013© [51]. Pixels were aggregated to match the Landsat resolution (30 m) using the nearest neighbour resampling. Finally, in addition to the original backscatter, we also calculated the dry season ratio image of HH/HV as it has been proved to be useful to land cover classification [52][53][54]. To ensure consistency across the study area, images were selected with the smallest possible time difference: a maximum of 17 days for both seasons (Table 2). Following download, the raw images were first converted to backscatter using the standard conversion equation: where σ 0 is backscatter, DN is the raw digital number and CF is a calibration factor (−83). Secondly, a Lee filter with a 7 × 7 window size was applied to minimise speckle effects [50], using Erdas IMAGINE 2013© [51]. Pixels were aggregated to match the Landsat resolution (30 m) using the nearest neighbour resampling. Finally, in addition to the original backscatter, we also calculated the dry season ratio image of HH/HV as it has been proved to be useful to land cover classification [52][53][54]. To increase the utility of the SAR data, we calculated a series of Gray-Level Co-Occurrence Matrix (GLCM) texture variables. GLCMs are a series of localised texture metrics that quantify the statistical properties of a layer over a moving window [55]. We calculated seven GLCM layers (mean, variance, homogeneity, contrast, dissimilarity, entropy, and second moment [56]). These statistics were calculated over both 3 × 3 and 9 × 9 windows, resulting in 15 layers per SAR backscatter (one backscatter + seven 3 × 3 GLCM layers + seven 9 × 9 GLCM layers).

Classification and Validation
Our analysis focused on identifying the optimal combination of seasonal and sensor information for savannah land cover mapping. Moving toward this aim, we developed a series of classification models consisting of combinations of the datasets described in the previous sections. This resulted in a total of 15 independent models (Table 3). Table 3. The 15 independent models tested in this study and the number of parameters in each. Land cover classifications were carried out using the machine learning algorithm of Random Forests which combines decision trees with bootstrapping and aggregation [57]. Random Forests has been proven to be a successful classifier for remote sensing imagery, due to their effective handling of correlated predictors and reduced tendency toward over-fitting [58]. Models were constructed using the "randomForest" package version 4.6 within the R statistical environment [59]. Training data were acquired through random sampling, with 2400 points (600 per class) representative of a 30 m Landsat pixel selected.

Model # Parameters Included Number of Parameters
The classified land cover maps were validated using the NGI aerial photographs and a two-stage stratified random sampling procedure following best-practice guidelines [60]. First, an initial sample of 100 points per class was collated (400 samples in total). The user's accuracy was then calculated, which, combined with the mapped cover proportions from the all-variable model (Model #9 in Table 3), allowed for the estimation of an appropriate stratified random samples size, with a target standard error of 0.01 [60][61][62]. This resulted in a final validation set of 1082 points, which were assigned a land cover class based on the visual interpretation of the NGI imagery. Finally, the model accuracy and class area estimates were adjusted using probability equations and best practice guidelines [60].

Results
The model that scored the highest overall accuracy of 91.1 ± 1.7% was the one that incorporated Landsat and PALSAR metrics from both the dry and wet seasons (Model #9, Table 3, Figure 4). Its accuracy statistics are shown in Table 4, and the resulting map can be seen in Figure 5a. This model included 130 parameters and the ten most important ones, according to the mean decrease in the Gini index [63], can be seen in Figure S1.
Another five models (#1, 2, 5, 6, and 10) also scored very high overall accuracies (from 90.4 ± 1.8% to 91 ± 1.7%), i.e., within the confidence interval of the best-performing model. The six top performing models were able to map the main savannah types with similar per class accuracies ( Figure 6; Table S2).
Remote Sens. 2018, 10, x FOR PEER REVIEW  8 of 18 3), allowed for the estimation of an appropriate stratified random samples size, with a target standard error of 0.01 [60][61][62]. This resulted in a final validation set of 1082 points, which were assigned a land cover class based on the visual interpretation of the NGI imagery. Finally, the model accuracy and class area estimates were adjusted using probability equations and best practice guidelines [60].

Results
The model that scored the highest overall accuracy of 91.1 ± 1.7% was the one that incorporated Landsat and PALSAR metrics from both the dry and wet seasons (Model #9, Table 3, Figure 4). Its accuracy statistics are shown in Table 4, and the resulting map can be seen in Figure 5a. This model included 130 parameters and the ten most important ones, according to the mean decrease in the Gini index [63], can be seen in Figure S1.
Another five models (#1, 2, 5, 6, and 10) also scored very high overall accuracies (from 90.4 ± 1.8% to 91 ± 1.7%), i.e., within the confidence interval of the best-performing model. The six top performing models were able to map the main savannah types with similar per class accuracies ( Figure 6; Table S2).
According to the "All Landsat All SAR" model (#9), woody vegetation and grassland were the most prevalent land cover classes, with 20,833 km 2 (95% CI: 165 km 2 , i.e., 51 ± 0.4% of the study area) and 13,485 km 2 (95% CI: 569 km 2 , i.e., 33 ± 1.5%), respectively. Cropland was the third largest (5930 ± 501 km 2 or 15 ± 1.2% of the study area), while non-vegetated land occupied 3338 km 2 (95% CI: 710 km 2 , i.e., 8 ± 1.7% of the study area).      Table 3); (b) error-adjusted area estimates and 95% confidence interval margins for the four land cover classes estimated from the best performing model; and (c) the 2010 high resolution (SPOT 5 and colour infrared aerial photography-based) land cover data of the South African National Geospatial Information (NGI) mapping agency [43]. Locations A, B and C are the example subsets that appear in Figure 7.
The most reliably mapped class was woody cover, with user's and producer's accuracies of 98.1% and 93.7%, respectively (Table 4). All classes demonstrated user's and producer's accuracies of over ~80%, with the exception of the user's accuracy for the non-vegetated land (67.3%). The main source of error came from the commission of this land cover type, which can often occur when attempting to map classes that cover smaller areas [32,64]. Moreover, there is a degree of confusion between the grassland and cropland classes, which resulted in the relatively large confidence interval for cropland, even larger for grassland. This is also anticipated as a number of crops cultivated in the region (e.g., maize) have similar spectral properties to certain grasses [65]. In general, well-irrigated (pivot) fields were well mapped ( Figure 7C), while smaller scale agriculture was less so ( Figure 7A).  Table 3); (b) error-adjusted area estimates and 95% confidence interval margins for the four land cover classes estimated from the best performing model; and (c) the 2010 high resolution (SPOT 5 and colour infrared aerial photography-based) land cover data of the South African National Geospatial Information (NGI) mapping agency [43]. Locations A, B and C are the example subsets that appear in Figure 7.
The most reliably mapped class was woody cover, with user's and producer's accuracies of 98.1% and 93.7%, respectively (Table 4). All classes demonstrated user's and producer's accuracies of over 80%, with the exception of the user's accuracy for the non-vegetated land (67.3%). The main source of error came from the commission of this land cover type, which can often occur when attempting to map classes that cover smaller areas [32,64]. Moreover, there is a degree of confusion between the grassland and cropland classes, which resulted in the relatively large confidence interval for cropland, even larger for grassland. This is also anticipated as a number of crops cultivated in the region (e.g., maize) have similar spectral properties to certain grasses [65]. In general, well-irrigated (pivot) fields were well mapped ( Figure 7C), while smaller scale agriculture was less so ( Figure 7A).
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 18 Figure 6. The class-level accuracy statistics for the 15 models tested.
Our best performing model (Figure 5a) is spatially similar to the 10 × 10 m land cover map produced by the NGI for the Northwest Province in 2010 using SPOT 5 data (Figure 5c) [43]. We used the NGI map to compare our results in the absence of ground truth data, other than the very highresolution aerial photography. The two maps follow similar general patterns with regards to the extent of woody cover, cropland and grassland. If the NGI data are considered as the reference, then our land cover map manifests a general underestimation of woody vegetation and cropland and a respective overestimation of grassland and non-vegetated land. We compared the NGI land cover with our best-performing model statistically: considering each land cover class separately, the data were first aggregated into percentage cover for 1 km grid cells. Pearson's correlation tests were then applied to the corresponding pixels. This is considered a robust approach for comparing land cover as it minimises errors caused by geolocation and pixel misalignment [66]. The correlations were the highest for the cropland class (94.4%), followed by woody vegetation (74.5%) and grassland (61.7%), while a relatively low correlation was found for the non-vegetated class (54.5%).
Overall, the Landsat-only models were more accurate than the PALSAR models. The latter performed marginally worse than their optical counterparts, with overall accuracies of 81.6% (±2.3%) for the SAR model with both wet and dry season parameters (All SAR); 81.3% (±2.3%) for the dry season model (dry SAR) and 73.7% (±2.6%) for the wet season model (wet SAR), respectively ( Figure  4). However, the PALSAR imagery were very effective in identifying woody cover, achieving very low commission errors: 3.4% for the "all SAR" and 4.2% for the "dry SAR". The "wet SAR" model was not as successful, scoring a much higher woody vegetation commission error of 13.6% (Table S2). All other classes generally had higher commission and omission errors than the Landsat-based models, with the exception of the producer's accuracy of the non-vegetated areas achieved by the "all SAR" and the "dry SAR" models (Table S2, Figure 6). Most interestingly, the results also show that the combination of the optical with the L-band SAR data improved the predictions, especially when the dry season PALSAR data were added to the models (Figures 4 and 5; Table S2). Our best performing model (Figure 5a) is spatially similar to the 10 × 10 m land cover map produced by the NGI for the Northwest Province in 2010 using SPOT 5 data (Figure 5c) [43]. We used the NGI map to compare our results in the absence of ground truth data, other than the very high-resolution aerial photography. The two maps follow similar general patterns with regards to the extent of woody cover, cropland and grassland. If the NGI data are considered as the reference, then our land cover map manifests a general underestimation of woody vegetation and cropland and a respective overestimation of grassland and non-vegetated land. We compared the NGI land cover with our best-performing model statistically: considering each land cover class separately, the data were first aggregated into percentage cover for 1 km grid cells. Pearson's correlation tests were then applied to the corresponding pixels. This is considered a robust approach for comparing land cover as it minimises errors caused by geolocation and pixel misalignment [66]. The correlations were the highest for the cropland class (94.4%), followed by woody vegetation (74.5%) and grassland (61.7%), while a relatively low correlation was found for the non-vegetated class (54.5%).
Overall, the Landsat-only models were more accurate than the PALSAR models. The latter performed marginally worse than their optical counterparts, with overall accuracies of 81.6% (±2.3%) for the SAR model with both wet and dry season parameters (All SAR); 81.3% (±2.3%) for the dry season model (dry SAR) and 73.7% (±2.6%) for the wet season model (wet SAR), respectively ( Figure 4). However, the PALSAR imagery were very effective in identifying woody cover, achieving very low commission errors: 3.4% for the "all SAR" and 4.2% for the "dry SAR". The "wet SAR" model was not as successful, scoring a much higher woody vegetation commission error of 13.6% (Table S2). All other classes generally had higher commission and omission errors than the Landsat-based models, with the exception of the producer's accuracy of the non-vegetated areas achieved by the "all SAR" and the "dry SAR" models (Table S2, Figure 6). Most interestingly, the results also show that the combination of the optical with the L-band SAR data improved the predictions, especially when the dry season PALSAR data were added to the models (Figures 4 and 5; Table S2). The most commonly employed approach of using dry season optical data did not prove to be more accurate than its wet season counterpart. The overall accuracies for both models were very The most commonly employed approach of using dry season optical data did not prove to be more accurate than its wet season counterpart. The overall accuracies for both models were very similar, with the wet model actually performing slightly better: 81.8 ± 2.3% for the wet season and 81.7 ± 2.3% for the dry season, respectively (Figure 4). The wet season model also scored lower grassland and cropland commission errors, as well as lower grassland and non-vegetated omission errors, than did the dry season model ( Figure 6, Table S2). Most importantly, when the wet season data were added to the 'dry' models, results always improved for all land cover types ( Figure 6, Table S2).
Similarly, when SAR data were used in the model, the addition of the wet season data improved the model performance for most land cover types. With the exception of the producer's accuracy for the grassland and the cropland and of the user's accuracy for the latter, in all other cases the per-class accuracies increased when the wet season data were added.

Discussion
It is commonly understood that savannah landscapes are difficult to map, mostly due to the low inter-class but high intra-class spectral variability, confounded by seasonal and within-pixel variation [32,67]. Our best performing model, incorporating 130 parameters and estimated from dry and wet season optical and radar data, was able to achieve a high overall accuracy (91.1 ± 1.7%), as well as user's and producer's accuracies ranging from 80% to 98% for woody cover, grasses and crops. These three land cover types together occupy 92% of the study area in southern Africa. Difficulties were identified and high commission errors were produced when attempting to discriminate the non-vegetated areas (32.7%). This was anticipated and explained by both the low coverage of the area by such land cover types (8%), as well as the spectral confusion with the fallow land [32,64]. The latter was also noted by Eggen et al. [67] in their Landsat-based Ethiopian study, where they report very low accuracies for the barren class (54% user's and 48% producer's accuracy).
Our results, produced from multi-sensor and multi-season data, were similar to the accuracy achieved using coarser resolution data by Mishra et al. [68] who employed MODIS annual and intra-annual vegetation indices in a southern African savannah (Central Kalahari, Botswana). They presented overall accuracies of 93% when attempting to map the six main morphological vegetation types (woodland, dense shrub land, open shrub land, very open shrub land, grassland and pans/bare areas), which were very similar to ours (91.1%). The very high accuracy achieved with the MODIS data is commendable, especially since the methodology suggested by Mishra et al. [68] does not involve an onerous multi-sensor undertaking. However, when a finer resolution mapping product is required, MODIS-based approaches are not sufficient, and more complicated sensor fusion techniques need to be considered: for instance the one proposed in our study, or Landsat-MODIS fusions (e.g., STARFM [69] and ESTARFM [70]).
The use of metrics was also identified as being beneficial for land cover mapping by Müller et al. [32] in a study on another savannah environment in the Brazilian Cerrado, which was however much wetter. They used Landsat time series to map cropland, pasture and natural savannah vegetation and achieved a very similar overall accuracy (93%). They make a note that such high accuracies were possible due to the temporal depth of their study period, which was 3 years. In our study, we had to resort to using a slightly longer period (5 years) due to data availability and cloud contamination issues (Figure 3). The use of a 5-year compositing epoch may result in some land cover changes occurring within our study period. However, most change processes in savannahs (e.g., shrub encroachment, overgrazing, and overexploitation of woody shrubs and trees for fuelwood) are gradual transitions, which reduces the risk of abrupt conversions.

Landsat or PALSAR? or Both?
Radar data are known for their ability to map certain savannah vegetation properties, such as woody and canopy structure and volume [26,71]. Here we tested their fitness for mapping different savannah land cover types. When it came to the ability to identify woody cover, it was the combination of PALSAR dry and wet season metrics that produced the lowest commission errors (3.4%) and outperformed all other optical-only model combinations. This is in support of a number of studies who tested the performance of SAR data when mapping woody vegetation or the percentage of woody cover. Naidoo et al. [7] compared L-band PALSAR dry-season data with Random Forest models incorporating Landsat bands and vegetation and textural indices from different seasons, in order to map a woody canopy cover in an area within the Kruger National Park in South Africa. They found that the SAR-only models outperformed the best Landsat-only model by 9% (R 2 = 0.81 and R 2 = 0.72, respectively).
Our results, however, also showed that the Landsat-based models were generally slightly better than their SAR-based counterparts in identifying the other three targeted savannah land cover types. Walker et al. [72], in their study in the southeastern Brazilian Amazon, along with Laurin et al. [73], who compared SAR and optical classifications of 8 land cover classes in the West African region, have found that Landsat-based models performed better than PALSAR-based ones. In a more similar savannah environment in the Limpopo Province of South Africa, Higginbottom et al. [28] also tested the performance of models incorporating Landsat and PALSAR data. They mapped the fractional woody vegetation cover and used the global annual PALSAR mosaics, rather than the actual data we use here (http://www.eorc.jaxa.jp/ALOS/en/palsar_fnf/fnf_index.htm). They found that the Landsat-based models outperformed the PALSAR-based ones. However, the amount of Landsat imagery available, which makes the use of metrics more effective [32], renders the comparison with the relatively sparse PALSAR archive somewhat disproportionate: the SAR only models and the respective dry and wet season GLCM texture variables were developed from data gathered over one year only (2008; Table 2), while hundreds of Landsat images were involved in the calculations of the spectral-temporal variability metrics (Table 1). Still, the top five variables of the best-performing model, out of a total of 130, were all calculated from the PALSAR data ( Figure S1), which is a clear indication of the importance of incorporating radar data in savannah land cover mapping efforts.
We found that multi-sensor combinations performed better than optical-or radar-only models: our dry season Landsat result (Figure 7(A5,B5,C5)) was improved by~10% when the dry season SAR data were included in the model (Figure 7(A3,B3,C3)), bringing the overall accuracy up from 81% to~91%, almost the maximum achieved by any model combinations. This is in agreement with Naidoo et al. [7] who found that the combination of multi-sensor data produces the highest overall accuracies, reporting an improvement of 8% and 17% from the PALSAR-only and Landsat-only models, respectively. Laurin et al. [73] also compared PALSAR and optical (Landsat and AVNIR-2) land cover classifications in a wet tropical area and reported an improvement from the multi-sensor integration. In another study in the southeastern Brazilian Amazon, Walker et al. [72] compared PALSAR-based model results to Landsat ones and concluded that multi-sensor and multi-temporal radar and optical data is the way forward.
Given the promising results achieved by [26,74] with the RADARSAT-2 C-band when mapping savannah woody vegetation structure and by [75] who combined the Sentinel-1 C-band and Sentinel-2 data (with improved spatial, spectral and temporal resolutions compared to the Landsat data used in this study) in order to map land cover in a Mediterranean environment (overall accuracy =~94%, k = 0.928), it is anticipated that the synergy between the two Sentinels should provide considerable support in the efforts to accurately map savannah land cover.

Dry or Wet? or Both?
Using dry season optical data alone to map savannah land cover types, as in the case of most land cover mapping studies [76,77], produced a very similar overall accuracy (~82 ± 2%) and very similar omission and commission errors to their wet season counterparts. Wet season data, however, are inherently problematic due to the fact that cloud contamination and Landsat data availability is often limited in many dryland regions [78]. It is, therefore, reasonable that most analysts favour the use of dry season data since the accuracies achieved are relatively high, especially for mapping woody vegetation (with commission errors in our study of 6%). However, we also found that the combination of dry with wet season metrics improved results (an overall accuracy improvement of~4% and per class statistics improvements throughout) and it is this approach that should be preferred if enough data are available. This is in agreement with Higginbottom et al. [28] who tested fractional woody cover estimation models with metrics estimated from the dry and the wet period and reported higher accuracies when data from both seasons were used. The higher temporal resolution of the Sentinel-2 data should therefore provide the means for improved savannah land cover mapping endeavours.
The PALSAR dry season model (Model #14 in Table 3) outperformed the wet season one (Model #15 in Table 3) in all mapped land cover classes (Table S2). It is also worth noting that the dry season HV backscatter outscored all other parameters as the best-performing model according to the mean variance Gini index ( Figure S1). However, the wet season HV polarisation data were missing for our area and period of study and, therefore, we are unable to directly compare the dry and wet season SAR model performances. Nevertheless, we did find that the addition of the wet season SAR data to their dry season counterparts increased the model performance (Figures 4 and 5; Table S2); consequently, if radar data from both seasons are available (e.g., Sentinel-1 C-band), this should be the preferred methodological approach for savannah land cover mapping.
An approach that has not been considered here but could potentially improve results is not to consider dry and wet season predictors separately, but to combine them. For example, a difference between a mean wet and a mean dry vegetation index could be calculated for the optical data and a similar mean index estimated from the different seasons and polarisations of the SAR data. On a similar note, another approach could be to use a phenological analysis to derive surrogate maps of the vegetated land cover types by analyzing the temporal evolution of an index (e.g., NDVI; [79]). High temporal-resolution data are required for this and, as a consequence of this, coarser resolution data, such as MODIS, have traditionally been employed [80]. However, new approaches that extract phenology estimates using high-resolution imagery from a single season [81], or that combine high spatial-resolution data (e.g., GF-1) with phenological parameters from high temporal-resolution sensors (e.g., MODIS) to classify land cover, have recently been developed [82]. To this end, our ability to now employ Sentinel-2 data, with their high repeat frequency, offers a promising outlook.

Conclusions
Savannah ecosystems cover a fifth of the Earth's surface and are one of the least studied and most difficult biomes to map, increasingly undergoing processes of degradation. The most recent attempts to map the savannah land cover have incorporated optical or radar EO data from the dry season or from different phenological periods. In this study, we focus on a southern African area to assess how accurately woody vegetation, grassland, cropland and non-vegetated areas can be mapped in order to provide assistance in future attempts to characterise and monitor savannah land cover. Given the overall accuracy and the per class statistics produced from the 15 models tested, as well as the spatial similarity of an independent finer-scale land cover map with our results, it is safe to conclude that the main cover types of woody vegetation, grassland and crops can be accurately mapped, when optical and radar data from the dry and wet seasons are integrated. The recent developments in EO sensors and the open-access availability of higher spectral, spatial and temporal resolution data (e.g., the European Space Agency's Sentinel-1 and Sentinel-2 satellites), offer promising pathways for accurate multi-sensor and multi-season savannah land cover mapping.