Non-seasonal Landsat based bare area gain detection in Botswana during 2002 to 2020 Period using Maximum Likelihood Classifier (MLC).

This paper estimates bare area gain detected using cloud free Landsat 7 (ETM+) and Landsat 8 (OLI) in Botswana. From 2002 to 2020, agricultural fields shrunk by 76.4%, while built-up increased by 49.2%, and bare areas increased from 3.32% to 7.03% (or 111.7%). There is a significant seasonal change in bare area detected reaching maximum during the dry season when there is little or no ground cover. In this study, the seasonality of bare area gain was overcome by only considering a bare area pixel to contribute to bare area gain if it exists during both the winter and summer months. The probability of bare area detection was 75.0% and probability of false detection 13.3% respectively. The 13% false detection tended to be built-up areas which had similar spectral characteristics as bare areas since most built-up areas have no ground cover. The bare area gain is driven by the high population growth rate of 3.4%. From 2001 to 2017, the population of the study area has increased by 34% and now accounts for 47% of the population of Botswana.


Introduction
The ever-increasing human population with accompanying economic drivers such as concentration of people in urban areas and excessive natural resources extraction is exerting stress on the planet leading to changes in the land cover such as vegetation cover loss (Potapov et al., 2012;Ying et al., 2017). Change in land cover is key component of global environment change and drives a wide range of ecological, hydrological, and climatic processes such as loss of forests, urbanization, and increase in agricultural activities (animal grazing and crops production) to increase food production (Foley et al., 2005;Vitousek et al., 1997). It also has a negative impact on flora and fauna especially on sensitive landscapes (Cook et al., 2012;Ringrose et al., 1996;Vittor et al., 2006).
Landsat data had extensively been used in Land Use Land Cover (LULC) studies from deforestation, fire scars, agriculture monitoring, urban sprawl, and wetlands loss (Daldegan et al., 2019;Jensen et al., 1995;Spruce et al., 2018;Zhu & Woodcock, 2014). Moderate Resolution Imaging Spectroradiometer (MODIS) and Sentinel-2 datasets had been used to monitor vegetation and land cover changes from 2005 and 2015 respectively (Jin & Sader, 2005;Schug et al., 2020;Xin et al., 2012). Although non-optical techniques are also used in LULC mapping especially where clouds are frequent (Nicolau et al., 2021), optical remote sensing has long been and is promising as a systematic and cost effective way of monitoring bare ground cover and gain (Barrett et al., 2020;Dematte et al., 2009;Hansen et al., 2014;Nwilo et al., 2020). For LULC classes detection, either full pixel or subpixel scheme is used for the image classification using statistical, regression, expert classifier, artificial neural networks, and other more sophisticated algorithms (Daldegan et al., 2019;Schug et al., 2020). Although the maximum likelihood classification assumes that the statistics for each class in each band are normally distributed can still produce acceptable results even if this assumption does not hold furthermore, it is easy to implement using full pixel scheme.
Desertification is defined as the depletion or destruction of the biological capability of the land (in arid, semi-arid and dry sub-humid areas) leading to desert-like condition caused by climate variation and human activities (UNCOD, 1977;United Nations, 1993). Desertification which results from an increase in bare areas is an issue of environmental concern in semi-arid to arid parts of the world covering 40% of the global land surface and accounting for about 15% of the human population (Verón et al., 2006).
The first estimates of global bare ground gain study representing 2000 to 2012 period using satellite imagery concluded that bare ground gain represent a very small fraction of the overall global land cover (93,896±9317 km 2 ) (Ying et al., 2017;Hansen et al., 2013;Schneider et al., 2010). The huge uncertainty in this estimation (almost 10%) is because land use extent and change is difficult to quantify. Further, most of this gain occurs in the East Asia and the Pacific (45%) with China alone accounting for (78%) this gains (Ying et al., 2017). Human induced bare ground gain accounted for 95% of the total of which 39% is accounted for by commercial and residential development (builtup), but in Sub-Sahara Africa (account for 3% of global bare ground gain), the highest proportion of the bare ground is due to resources extraction followed by built-up. It has been estimated that 45,945.08 km 2 of vegetation cover was lost in northern Nigeria including the Lake Chad area (at a rate of 0.96% and 0.28%) between 1984 and 2016 (Nwilo et al., 2020). There is lack of accurate assessment of the status, change, and trend of desertification, this hinders developing global actions to prevent and eradicate the problem especially in one of the most seriously affected countries like China (Yang et al., 2005). Ying et al (2017) determined that human-induced bare ground gain accounted for 95% of the total and for Sub-Sahara Africa, this is made up of the following components: 1. commercial and residential development (37%), resource extraction/mining (32%), infrastructure development (4%), transitional (15%), Natural bare ground (12%) and greenhouses (0%).
The study presented here detects bare areas by aggregating bare areas detected during a wet season (good vegetation cover) and dry season (poor vegetation cover) in one year and comparing the bare areas change temporally (2002 to 2020). A true bare ground pixel is selected only if the pixel being considered is bare during both winter and summer months.
Studies have shown that accuracy assessment results from homogeneous land cover classes are generally higher than those from heterogeneous classes (Gessner et al., 2013); Latifovic & Olthof, 2004;Friedl et al., 2010). In the arid to semi-arid areas such as the Kalahari, bare area is actually 30% herbaceous and 70% bare (Gessner et al., 2013). This spectral mix causes confusion in LULC classification i) rural areas where bare areas, built-up, crop fields and grazing pasturelands classes are not distinct hence tend to generate a spectral mixture, ii) during the winter months, the agricultural fields and grazing lands are bare (no crops, leafless trees and shrubs, and no grass cover which is normally annual iii) urban built-up land use has a lot of open bare areas (unpaved); hence these tend to overestimate bare land class and also lowers the overall accuracy of LULC classification. Because of these inconsistencies in estimating bare areas and hence the rate of desertification at local scale, we examine the aggregating bare areas detected during two extreme conditions (winter/summer) in a year (2020) and compare to a previous period (2002)

Study Area
The study area is in Botswana (Figure 1), a relatively flat landlocked country in Southern Africa of size 581, 730 km 2 with a population of about 2.3 million people (Statistics Botswana, 2018). The country faces reoccurring droughts and a threat of desertification resulting from communal grazing which accounts for 71% of the country's land (Darkoh, 1998). The greater Gaborone area is the area south of Botswana which comprise of the capital city Gaborone and several major towns and villages accounting for 38% of the total population of Botswana (Central Statistics Office, 2000;Statistics Botswana, 2018;Tsheko, 2006).
The annual climate ranges from months of dry temperate weather during winter (April-September) humid subtropical weather interspersed with drier periods of hot weather during summer (October-March). Annual rainfall, brought by winds from the Indian Ocean, averages 500 mm, for the study area (Jonsson, 2004). In 2006, FAO estimated that bare land/built up (bare -land covered by bare soil with less herbal vegetation, such as clearing, mining and rocky area while built up -land covered build up area such as settlement and urban, and industrial area) accounted for 1.5% of Botswana land cover (FAO, 2020).
The general vegetation of the country is savanna grassland with yellow or light brown grass cover (turning green after rains) and woody plants. The savanna ranges from acacia shrub savanna in the southwest through acacia thornbush and tree savanna "parkland" into denser woodland and eventually forest as one moves north and east (Jonsson, 2004).

Methodology
Four (4) cloud free Landsat images, two (2) Landsat 8 scenes (Path 172/Row 077) acquired 31 st January and 7 th June 2020, and two (2) Landsat 7 scenes (Path 172/Row 077) acquired 21 st January and 30 th June 2002 were used in this study. Preprocessing was done with the objective to convert the digital numbers (DNs) into scaled surface reflectance values (Du et al., 2002;Ihlen & Zanter, 2019;Teillet et al., 2001). The image radiometric calibration was implemented in ENVI software Radiometric Correction Module © Exelis Visual Information Solutions using parameters from the Landsat metadata files (Chavez, 1988;Schaepman-Strub et al., 2006). The Methodology is summarized in the flow chart Figure 2 below.

Features Extraction
Six (6) classes were defined for the imagery, Agricultural fields, trees, bare, built-up, shrubs, and water. The mean reflectance values of following nine (9) and eight (8) features were calculated and extracted for each of the four classes above for both the Landsat 8 and Landsat 7 respectively as shown in Table 1 below.

Maximum likelihood classification
The images features were used in the classifications using the supervised Maximum Likelihood Algorithm ( years i and j (Wu et al., 2017). In this study only the bare area class was analysed. Where i=2002, j=2020 and b=bare, then x i =x i,b and y j =y j,b . The change detection (bare areas) was calculated as ∆ xy = x ^y for i=j and ∆ xy = x ⊕ y for i≠j.

Accuracy assessment
For validation, a total of 150 random pixels were selected for each of the four images used and high-resolution SPOT images, google map, Botswana Department of Surveys and Mapping orthophotos and field visits were used for ground truthing. Zhu & Woodcock (2014) noted that the high spatial resolution images at the temporal frequency as Landsat data, can be helpful in determining land cover change at longer time intervals. The field visits were important to resolve some validation points as LULC is highly dynamic especially in and around urban areas and dried water pans or ponds. To assess the accuracy of the detected bare areas from above process, twentyfive (25) random points were again selected for validation using the high-resolution data and field visit.

Error Matrices
Once classification of the data sets was achieved, the following matrices (overall accuracy OA, Kappa Coefficient KC, Producer Accuracy PA, and User Accuracy UA) were calculated for all the tests. For the bare areas detection an error matrix was used to estimate classification accuracy and further, probability of detection (PoD) and probability of false detection (PoFD) were also calculated.
PoD is calculated as hits divided by hits and misses while PoFD is false negatives divided by correct negatives plus false alarms.

Results and Discussion
There is complexity of the built-up spectral signature which is a spectral mixture (rural built-up, built-up, bare, industrial built-up), this resulted in built-up being confused with agriculture fields and bare land. The built-up producer accuracy was 86.7% while the user's accuracy was only 44.8% respectively. Both the producer accuracy and user's accuracy for the bare areas land class was 77.8% respectively.
The classified images for 2002 show significant influence of seasonal changes on all the LCLU classes including the water bodies (shrunk surface area due to usage and evaporation -since there is no rain fall during the dry season). One of the significant changes is the increase in bare areas especially around built-up areas and areas of agricultural activity (due to barren fields after harvesting and livestock overgrazing since livestock is 100% fed on pasture with minimal feedlots) as shown in   Botswana (Central Statistics Office, 2005, 2011. These bare areas increase significantly during the dry months and the built-up signature also increases due to the dry conditions and human activities in settlements (agricultural actives, firewood collection and forest fires).  Figure 10 shows that from 2002 to 2020, water bodies shrunk by 14.5%, agricultural fields shrunk by 76.4%, bare areas increased by 111.7%, shrubs shrunk by 1.6%, built-up increased by 49.2% and trees shrunk by 13.9%. As population increased in the greater Gaborone area, agricultural fields are converted into residential and industrial areas, it is estimated that the population growth for the study area ranges from 2.4% to 3.4% (Central Statistics Office, 2005, 2011. Bare areas have significantly increased due to land cleared for built-up, firewood collection, fallow and unused agriculture fields and overgrazing by livestock. Although the shrubs land has shrunk slightly since 2002, the natural shrub land has largely been replaced by acacia thornbush (Jonsson, 2004). The classification overall accuracy and kappa coefficient was 82.3% and 0.79 for 2002 and 81.6% and 0.78 for 2020 as shown in Table 2 and Table 3. In both cases the water LULC had a user and producer accuracies of 100% with the built-up having the lowest classification accuracy, although the bare land was satisfactorily classified, the built-up was mostly misclassified as bare land and field visits confirmed that indeed most built-up especially in rural areas had a large component of bare area pixels hence the spectral signature was closer to that of bare area. The bare area user accuracy was 87.0% and 77.8% for 2002 and 2020, respectively.   Table 3 Accuracy assessment 2020 Figure 9 and Figure 10 show that in all the six agricultural districts, bare area has significantly increased during the period 2002 to 2020. The accuracy of detected areas has an overall classification accuracy of 82.6% with a kappa coefficient of 0.72. The probability of detection was 75.0% and probability of false detection was 13.3% respectively, this means that a detected bare area has a 75% chance of being an actual bare area and there is a 13% chance that a detected bare area is not an actual bare area on the ground. This 13% tended to be built-up areas which have a lot of bare pixels. There is a lot of seasonal changes in bare areas due to the long and dry season during which there is very little or no ground cover. In this study, this effect was normalized by only considering a bare pixel as a true bare area if the pixel being considered is bare during both winter and summer months. For other studies such as soil erosion, it might be advisable to consider bare ground during the winter months since most soil is lost during the first rains after the long dry winter.

Conclusions
Agricultural fields shrunk by 76.4%, while built-up increased by 49.2%, and bare areas increased from 3.32% to 7.03% (111.7%) from 2002 to 2020. From 2001 to 2017, the population of the study area increased by 34% and now accounts for 47% of the population of Botswana. There seasonal changes in bare areas reaching maximum during the long and dry season when there is very little or no ground cover. The effect of bare area seasonality was normalized by only considering a bare ground pixel as being true only if the pixel is bare area during both winter and summer months.
Detecting bare area during the dry winter months represents temporary or seasonal bare area which tends to overestimate bare ground gain. The bare area increase is attributed to the pressures in land use due to the population growth rate of 3.4% in the study area. Bare ground gain should be detected using both dry and wet months data to avoid overestimation of the gain. Bare ground detected in the dry winter months should be used in soil erosion studies since most soil is lost during the first rains after the long and dry winter.

Recommendations
To avoid the misclassification of bare areas with bare pixels in built-up areas, it suggested that for further research work be carried out using sub-pixel classification methods.

Acknowledgement
The author is grateful to USGS for Landsat images, Botswana Department of Surveys and Mapping for the orthophotos, Botswana Ministry of Agricultural Developments and Food Security of the agricultural districts shape file, Google Maps and the University of Agriculture and Natural Resources for the field work resources. The author is grateful to the anonymous reviewers for their constructive comments and suggestions.