Multi-Temporal Landsat Images and Ancillary Data for Land Use/Cover Change (LULCC) Detection in the Southwest of Burkina Faso, West Africa

Abstract

Accurate quantification of land use/cover change (LULCC) is important for efficient environmental management, especially in regions that are extremely affected by climate variability and continuous population growth such as West Africa. In this context, accurate LULC classification and statistically sound change area estimates are essential for a better understanding of LULCC processes. This study aimed at comparing mono-temporal and multi-temporal LULC classifications as well as their combination with ancillary data and to determine LULCC across the heterogeneous landscape of southwest Burkina Faso using accurate classification results. Landsat data (1999, 2006 and 2011) and ancillary data served as input features for the random forest classifier algorithm. Five LULC classes were identified: woodland, mixed vegetation, bare surface, water and agricultural area. A reference database was established using different sources including high-resolution images, aerial photo and field data. LULCC and LULC classification accuracies, area and area uncertainty were computed based on the method of adjusted error matrices. The results revealed that multi-temporal classification significantly outperformed those solely based on mono-temporal data in the study area. However, combining mono-temporal imagery and ancillary data for LULC classification had the same accuracy level as multi-temporal classification which is an indication that this combination is an efficient alternative to multi-temporal classification in the study region, where cloud free images are rare. The LULCC map obtained had an overall accuracy of 92%. Natural vegetation loss was estimated to be 17.9% ± 2.5% between 1999 and 2011. The study area experienced an increase in agricultural area and bare surface at the expense of woodland and mixed vegetation, which attests to the ongoing deforestation. These results can serve as means of regional and global land cover products validation, as they provide a new validated data set with uncertainty estimates in heterogeneous ecosystems prone to classification errors.

1. Introduction

Remote sensing plays an important role in the management of the earth’s surface by providing spatio-temporal information on land use/cover (e.g., water, forest, bare area, and cropland). Change in land use/cover (LULC) affects the global environment [1], its biodiversity [2], local, regional and global climate [3,4], and also can accelerate, among others land degradation, which reduces ecosystem services and functions [5]. Monitoring land use/cover change (LULCC) is, therefore, relevant for sustainable landscape and environmental management. For instance, in regions like West Africa, where LULC is known to change rapidly [6], regular map updates could lead to better estimation of deforestation and land degradation rates, which are key components of the UN REDD (Reducing Emissions from Deforestation and Forest Degradation) program and UNCCD (United Nation Convention to Combat Desertification) strategy “Zero Net Land Degradation”, respectively.

Remote sensing provides an effective means for mapping LULC, especially over large areas [7]. Different types of classifications have been applied for LULC mapping [8,9,10,11]. Some of these classifications are based on mono-temporal image (e.g., [12,13]), multi-temporal images (e.g., [10,14]), or in combination with ancillary data (e.g., [11,15,16]) such as biophysical variables (e.g., elevation, slope and soil types). Mono-temporal classification is widely used in literature (e.g., [6,12,17,18]); this analysis focuses on a single date image for LULC mapping. The processing of single date image is faster as compared to multi-temporal classification. In multi-temporal classification, bands from more than one date, season, or year are combined and classified [19]. This approach involves the use of more images, which increases processing time. Multi-temporal classification has been used to classify the LULC of various landscapes (e.g., [10,14,20]).

A comparison of mono-temporal and multi-temporal classifications in different areas revealed contradictory results. For instance, in two states in the United States of America (USA), Lunetta and Balogh [10] used single date Landsat image against two dates of imagery and found that the latter increased wetlands identification accuracy from 69% to 88%. Key et al. [14] also noticed a positive effect of multi-temporal classification for discriminating individual tree species in a temperate hardwood forest (West Virginia University forest). Contrary to the aforementioned studies, Langley et al. [19] compared single date and multi-temporal satellite image classifications in semi-arid grassland and found out that single date classification produced more satisfactory results than multi-temporal classification. The integration of ancillary data with mono-temporal image has been found in some places suitable for improving LULC classification accuracy (e.g., [11,15]).

Although mono-temporal and multi-temporal classifications as well as the addition of ancillary data have been compared in elsewhere, such investigations remain rare in West Africa, especially in the Sudanese savannah zone. In addition, determining an efficient LULC mapping method is crucial for performing accurate post-classification change analysis and to monitor environmental degradation in this region, which is particularly exposed to effects of climate change and population growth. Furthermore, an accurate LULCC map could be an efficient decision-making support tool for mitigating carbon emission and for proper environmental change management, such as forest cover restoration through reforestation campaigns.

LULCC assessment in West Africa is challenging, as the region lacks atmospherically undisturbed and cloud free satellite images, which is a limiting factor for LULC mapping. Beside, lack of ground truth information for past years limits assessments of historical LULC (e.g., [21]) as well as validation of LULCC map (e.g., [6,22,23]).

The present study aimed at comparing mono-temporal and multi-temporal LULC classifications as well as their combination with ancillary data and to determine LULCC across the heterogeneous landscape of southwest Burkina Faso using accurate classification results.

2. Materials and Method

2.1. Study Area

The study area is located in the Black Volta basin of Burkina Faso with an area of about 5120 km² (Figure 1) and belongs to the South-Sudanese climatic zone that is characterized by two main seasons. The rainy season extends from May to October, and the dry season occurs from November to April. The average monthly temperature ranges from 26 °C to 32 °C, and the average annual rainfall in the period 1981–2012 was estimated to be 862.87 mm based on the data collected from the national direction of meteorology of Burkina Faso. Rainfall is marked by high inter-annual variability, and the vegetation is Sudanese savannah of which dominant woody species are Gardenia sp, Combretum micranthum, Parkia biglobosa, Vitellaria paradoxa, Bombax costatum, Berlinia grandifotia, among others.

The main livelihood activity of the local population is agriculture, which is characterized by low inputs (e.g., fertilizers) with slash and burn being the frequently applied farming practice [6]. According to Reenberg et al. [24], the typical agricultural practice is small-scale subsistence farming on rainfed fields, while the main crop productions are cotton and cereals (e.g., sorghum, maize and millet). Rice is cultivated as dry-rice in topographic depression or as irrigated wet-rice close to lakes and dammed-up rivers [13].

Figure 1. Situation of the study area.

Figure 1. Situation of the study area.

2.2. Data and Pre-Processing

2.2.1. Landsat Images

Landsat TM data (30 m × 30 m resolution, scene 196/052) were used for LULC classification. In total, ten monthly images were downloaded from United States Geological Survey (USGS) website for the years 1999, 2006 and 2011 (Table 1) using a criterion of cloud cover less than 10%. Five Landsat bands were considered in this study: blue (0.45–0.52 µm), green (0.52–0.60 µm), red (0.63–0.69 µm), near infrared (0.76–0.90 µm) and middle infrared (1.55–1.75 µm).

The Landsat data were calibrated and corrected to remove atmospheric influences. Firstly, the digital numbers (DN) values of the TM images were calibrated into radiance based on information from the meta data files provided by USGS (e.g., sun elevation, acquisition date). Afterwards, the radiance data were converted to surface reflectance using the ENVI 5 FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) module, which has shown to be efficient for satellite images atmospheric correction (e.g., [25]). This module requires information such as sensor altitude, initial visibility, atmospheric model, aerosol model, flight date, pixel size, and scene center location [26]. Some of these parameters were obtained from meta data files, while others were already included into the module such as atmospheric and aerosol models which were set to tropical and rural, respectively, during the atmospheric correction process to match with the study area environment.

Image-to-image co-registration was performed in order to ensure good alignment of pixels in the respective images. A root mean square error of less than one pixel was achieved for all the co-registrations. The images were already georeferenced to the Universal Transverse Mercator (UTM) projection WGS84 zone 30 north.

Table 1. Landsat TM images of 2011, 2006 and 1999.

2011	2006	1999
3 March	31 October	20 October
6 May	16 November	14 February
7 June	18 December	-
9 July	-	-
29 October	-	-

Table 1. Landsat TM images of 2011, 2006 and 1999.

2011	2006	1999
3 March	31 October	20 October
6 May	16 November	14 February
7 June	18 December	-
9 July	-	-
29 October	-	-

2.2.2. Environmental Ancillary Data

Ancillary data were also included in the LULC classification process (Table 2): slope, elevation, soil types and geomorphology. Slope and elevation were derived from version 2 of Digital Elevation Model (DEM) of ASTER, which has a spatial resolution of 30 m × 30 m. The DEM was georeferenced to UTM WGS84 zone 30 north. Vector layers (shapefiles) for geomorphology and soil types were obtained from the National Soil Office of Burkina Faso (BUNASOL) at a scale of 1/500,000. They highlighted different geomorphological units and soil types found in the study area. The shapefiles were projected to UTM WGS84 zone 30 north, then rasterized and resampled to 30 m × 30 m using ArcGIS 10.1.

Table 2. Ancillary data used in LULC classification.

Ancillary data	Source	Description	resolution
Elevation	ASTER	Height	30 m × 30 m
Slope	ASTER	degree	30 m × 30 m
Soil	BUNASOL-BF	Soil types	30 m × 30 m
Geomorphology	BUNASOL-BF	geomorphological units	30 m × 30 m

Table 2. Ancillary data used in LULC classification.

Ancillary data	Source	Description	resolution
Elevation	ASTER	Height	30 m × 30 m
Slope	ASTER	degree	30 m × 30 m
Soil	BUNASOL-BF	Soil types	30 m × 30 m
Geomorphology	BUNASOL-BF	geomorphological units	30 m × 30 m

2.2.3. Reference Data Sources

High-resolution images (Table 3) and field data were used to collect LULC reference samples to train and validate the Landsat images classifications. The high-resolution images were acquired as close as possible to the acquisition dates of the Landsat TM images. RapidEye image (from April 2011) was collected from the RapidEye Science Archive Team (RESA) of the German Aerospace Center (DLR) at level 3A (i.e., orthorectified with a spatial resolution of 5 m × 5 m and georeferenced to UTM projection). In addition, Quickbird image (2.4 m × 2.4 m) was acquired from October 2012, and scanned aerial photo (June 1999) was obtained from the National Geographic Institute of Burkina Faso (IGB) with a resolution of 2.3 m × 2.3 m. High-resolution images of Google Earth (from October 2006 and November 2007) were also used in this research. Apart from the RapidEye image that was already corrected, the other spatial data were geometrically adjusted (co-registration) to the Landsat images and georeferenced to UTM WGS84 zone 30 north.

A field campaign was conducted in October 2013 to collect LULC ground truth samples with a handheld Global Positioning Systems (GPS) using the projection system UTM WGS84 zone 30 north. LULC areas that remained stable since 2011 were sampled based on local population knowledge. Five broad LULC classes were identified (Table 4) using a modified LULC classification scheme of the FAO [27]: Woodland, bare surface, agricultural area, water and mixed vegetation. Homogeneous areas of 30 × 30 m were surveyed for each LULC (to match Landsat pixel), and the coordinate of the center recorded. In total, 150 samples were recorded across the study area.

Table 3. High-resolution images used in the analysis.

Images	Date	Resolution	Extent	% of study area covered
RapidEye	April 2011	5 m × 5 m	625 km2	12.20
Quickbird	October 2012	2.4 m × 2.4 m	25.7 km2	0.50
Aerial photo	June 1999	2.3 m × 2.3 m	188 km2	3.70
Google Earth image 1	November 2007	2.4 m × 2.4 m	306.9 km2	6.00
Google Earth image 2	October 2006	2.4 m × 2.4 m	309.8 km2	6.05

Table 3. High-resolution images used in the analysis.

Images	Date	Resolution	Extent	% of study area covered
RapidEye	April 2011	5 m × 5 m	625 km2	12.20
Quickbird	October 2012	2.4 m × 2.4 m	25.7 km2	0.50
Aerial photo	June 1999	2.3 m × 2.3 m	188 km2	3.70
Google Earth image 1	November 2007	2.4 m × 2.4 m	306.9 km2	6.00
Google Earth image 2	October 2006	2.4 m × 2.4 m	309.8 km2	6.05

Table 4. Adapted LULC classification scheme, modified from FAO [27].

Non-Modified	Modified
Level 1	Level 2	Adopted LULC classes
Vegetated	Woodland	Woodland
	Mixture of grasses, shrubs and trees	Mixed vegetation
	Cultivated area	Agricultural area
Non-vegetated	Bare land	Bare surface
	Built up
	Tarred road
	Rock
	Rivers	Water
	Artificial water bodies
	Lakes

Table 4. Adapted LULC classification scheme, modified from FAO [27].

Non-Modified	Modified
Level 1	Level 2	Adopted LULC classes
Vegetated	Woodland	Woodland
	Mixture of grasses, shrubs and trees	Mixed vegetation
	Cultivated area	Agricultural area
Non-vegetated	Bare land	Bare surface
	Built up
	Tarred road
	Rock
	Rivers	Water
	Artificial water bodies
	Lakes

2.3. LULC Classification

2.3.1. Image Combinations

In order to assess the classification of mono-temporal and multi-temporal satellite data as well as their combination with ancillary data, four spatial data set combinations were produced using Landsat images and ancillary data (Table 5). These images combinations are: (1) mono-temporal image; (2) mono-temporal image plus ancillary data; (3) multi-temporal images; and (4) multi-temporal images plus ancillary data. For the combination of mono-temporal image with ancillary data (Table 5, column 2), only the image that produced the highest overall accuracy amongst all mono-temporal classifications (Table 5, column 1) was selected.

Table 5. Different approaches of image combinations used for the year 2011.

	Mono-Temporal Image	Mono-Temporal Image Plus Ancillary Data	Multi-Temporal Images	Multi-Temporal Images Plus Ancillary Data
Landsat bands	Five images (October, July, June, May and March)	Image which achieved the highest accuracy in mono-temporal classification	Five images (October, July, June, May and March)	Five images (October, July, June, May and March)
Ancillary data		Elevation		Elevation
		Slope		Slope
		Geomorphology		Geomorphology
		Soil types		Soil types

Table 5. Different approaches of image combinations used for the year 2011.

	Mono-Temporal Image	Mono-Temporal Image Plus Ancillary Data	Multi-Temporal Images	Multi-Temporal Images Plus Ancillary Data
Landsat bands	Five images (October, July, June, May and March)	Image which achieved the highest accuracy in mono-temporal classification	Five images (October, July, June, May and March)	Five images (October, July, June, May and March)
Ancillary data		Elevation		Elevation
		Slope		Slope
		Geomorphology		Geomorphology
		Soil types		Soil types

2.3.2. Classification Algorithm: Random Forest Classification Algorithm

Supervised classifications of the Landsat images were performed using Breiman’s [28] non-parametric Random Forest (RF) classifier. RF is a machine learning algorithms [29] that can incorporate diverse sources of data (e.g., biophysical and remotely sensed data) [30], and it was found more accurate than other classifiers such as maximum likelihood classification (e.g., [31,32]), support vector machine (e.g., [31]) and neural networks (e.g., [33]).

RF is an ensemble of classification trees in which each tree contributes with a single vote for the assignment of the most frequent class to the input vector (X), ${\widehat{C}}_{rf}^B=majority vote {\left\{{\widehat{C}}_b\left(X\right)\right\}}_1^B$ [28,34], where ${\widehat{C}}_b\left(X\right)$ is the class prediction of the bth random forest tree, and B the total number of trees.

In fact, RF collects different subsets of training data to grow the trees (Ntree) using bagging [35] which is a technique used for training data creation by randomly resampling the original dataset with replacement (i.e., with no deletion of the data selected from the input sample for generating the next subset) [34,36]. A second random sampling is operated by RF to select a subset of predictive variables (Mtry) for the division of every node, which reduces the generalization error [29,37]. This randomness included in RF process decreases the correlation between trees in the forest and increases accuracy [38]. RF also computes the contribution of each variable to the classification using Mean Decrease Gini (MDG) and Mean Decrease Accuracy (MDA). However, Nicodemus [39] noticed that MDG is sensitive to within-predictors correlation and differences in category frequencies, while MDA is robust to these data characteristics.

Field data (only for 2011) and visual interpretation of the high-resolution images provide the possibility to collect a first set of LULC truth points (for each years) to train the RF classifier. Polygons of homogeneous pixels were drawn around each truth point for each LULC class and saved as vector layer of training areas. Landsat pixels that overlap the training areas were then used to train the RF classifier using R statistical software. The number of trees built was set to 800, and the number of selected variables at each node split was, by default, the square root of the total number of variables in each image combination.

2.4. LULCC Mapping: Post Classification Change Detection

Pixel-based post-classification change detection was performed using the classified maps of 1999 and 2011. Post-classification change detection consists in comparing individual LULC maps from two different periods in order to determine change areas. It emphasizes on “from to” change detection technique, which provides a change matrix and enables tracking the trajectory of each pixel between the two time steps of observation. In the present study, the post-classification change detection was performed in ArcGIS 10.1 where the LULC maps of 1999 and 2011 were combined into one raster file using the function “combine”. This enabled the detection of 25 LULC conversions or LULCC classes using the attribute table of the produced raster. However, four generalized categories of LULCC classes were considered: stable natural vegetation (SNV), natural vegetation loss (NVL), stable non-natural vegetation (SNNV), and other change (OC). Based on a new field created in the attribute table of the combined raster file, the four LULCC classes were assigned to each pixel following the rules shown in Table 6 below.

Table 6. LULC conversions included in the four LULCC classes between 1999 and 2011.

Name	Stable Natural Vegetation	Natural Vegetation Loss	Stable non-Natural Vegetation	Other Change
Change classes	Stable woodland	Woodland to other LULC unless mixed vegetation	Stable agricultural area	Agricultural area to all other LULC
	Stable mixed vegetation	Mixed vegetation to other LULC unless woodland	Stable bare surface	Bare surface to all other LULC
	Woodland to mixed vegetation
	Mixed vegetation to woodland
			Stable water area	Water to all other LULC

Table 6. LULC conversions included in the four LULCC classes between 1999 and 2011.

Name	Stable Natural Vegetation	Natural Vegetation Loss	Stable non-Natural Vegetation	Other Change
Change classes	Stable woodland	Woodland to other LULC unless mixed vegetation	Stable agricultural area	Agricultural area to all other LULC
	Stable mixed vegetation	Mixed vegetation to other LULC unless woodland	Stable bare surface	Bare surface to all other LULC
	Woodland to mixed vegetation
	Mixed vegetation to woodland
			Stable water area	Water to all other LULC

2.5. Accuracy and Area Assessment

2.5.1. Sampling Design

The sampling design is the protocol for selecting the subset of spatial units that will form the basis of the accuracy assessment [40,41]. In this study, pixels were used as spatial units (30 m × 30 m). Stratified random sampling was applied for collecting a second set of reference points to test RF classifications and validate the LULCC map; this is because the statistics computed in this study were based on the adjusted error matrix that requests a probability sampling for collecting reference data [42]. Stratified random sampling is a probability sampling design and a key element of a statistically rigorous assessment [41]. It enables statistical inference for computing estimates with confidence intervals.

Here, the different LULC and LULCC classes were considered as strata. Stratification is conducted when the strata are of interest for reporting results (e.g., accuracy and LULC class area), and to improve the precision of the accuracy and area estimates [40]. The strata sample allocation can be done either by equal sample size (same sample size for each class) or by proportional sample size (sample proportional to the spatial extent of each class in the map). The first method favors user’s accuracy against overall and producer’s accuracies [43], while, in the second method, the standard errors of estimating producer’s and overall accuracies become smaller as compared to equal allocation [40]. In order to take advantage of both sample allocation methods, the recommendation of Olofsson et al. [40] to increase the sample size of the rarer classes was followed in this study.

2.5.2. Response Design

The response design is the protocol for determining the reference LULC classification of a sampling unit [41]. The high-resolution images were used as reference data (Table 3).

For each year, the portions of the classified Landsat data that overlapped the high-resolution images were extracted and converted to polygon vector layer (in ArcGIS 10.1) from which the RF training areas were excluded. Each class was isolated as an individual vector layer based on which a set of randomly selected pixels were generated for each stratum (class). Each pixel received a LULC label by visual interpretation of the high-resolution images. The labels of each reference pixels on the maps and the high-resolution images were then compared to produce error matrices. Altogether, 371, 328 and 345 pixels were selected as samples for the years 2011, 2006 and 1999, respectively. The samples allocated to each stratum are shown in Table 7.

For the LULCC of the period 1999–2011, the overlapping area between LULCC map, aerial photo of 1999 and high-resolution images of 2011 were considered for reference pixels collection following the same procedure as above. However, Landsat images (1999 and 2011) were included especially for LULCC classes that were underrepresented in the high-resolution images. In all, 300 LULCC reference pixels were selected based on the classes established in Section 2.4. The numbers of samples are given in Table 8. Finally, the labels of the reference pixels on the change map and also on the reference images of 1999 and 2011 were compared to build an error matrix.

Table 7. Sample size allocated to each LULC class in 1999, 2006 and 2011.

	Sample Allocated
LULC classes (strata)	2011	2006	1999
Water	42	40	42
Woodland	100	97	111
Bare surfaces	59	45	45
Mixed vegetation	75	78	89
Agric. Area	95	68	58
Total column	371	328	345

Table 7. Sample size allocated to each LULC class in 1999, 2006 and 2011.

	Sample Allocated
LULC classes (strata)	2011	2006	1999
Water	42	40	42
Woodland	100	97	111
Bare surfaces	59	45	45
Mixed vegetation	75	78	89
Agric. Area	95	68	58
Total column	371	328	345

Table 8. Sample size allocated to each LULCC class.

LULCC Classes (strata)	Sample Allocated
Stable natural vegetation	125
Natural vegetation loss	56
Stable non-natural vegetation	42
Other change	77
Total	300

Table 8. Sample size allocated to each LULCC class.

LULCC Classes (strata)	Sample Allocated
Stable natural vegetation	125
Natural vegetation loss	56
Stable non-natural vegetation	42
Other change	77
Total	300

2.5.3. Analysis

Accuracy assessment

In order to determine the accuracy of a classified map with q categories, an error matrix is constructed [44]. Following the suggestions of Olofsson et al. [42] for assessing the accuracy of LULCC maps based on stratified random sampling strategies and pixel-based classifications, the error matrices of both LULCC and LULC classification were adjusted by the area of each category on the maps, and error matrices based on area proportions (${\widehat{P}}_{ij})$were produced:

$${\widehat{P}}_{ij}=W_i\frac{n_{ij}}{n_{i+}}$$

(1)

where $W_i$ is the proportion of area of category $i$ in the map, $n_{ij}$is the number of samples mapped as $i$ and belonging to category $j$ in the reference data, $n_{i+}$ is the number of samples mapped as category $i$ in the map. In the adjusted error matrix (Table 9), each cell element ${\widehat{P}}_{ij}$ indicates the probability that a randomly selected area is classified under category $i$ in the image and under category $j$ in the reference data [45].

Table 9. Adjusted error matrix of estimated area proportions.

			Reference
		1	2	…	j	…	q	Total
Map	1	${\widehat{P}}_{11}$	${\widehat{P}}_{12}$	…	${\widehat{P}}_{1j}$		${\widehat{P}}_{1q}$	${\widehat{P}}_{1+}$
	2	${\widehat{P}}_{21}$	${\widehat{P}}_{22}$	…		…		${\widehat{P}}_{2+}$
	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
	i	${\widehat{P}}_{i1}$	${\widehat{P}}_{i2}$	…	${\widehat{P}}_{ij}$	…	${\widehat{P}}_{iq}$	${\widehat{P}}_{i+}$
	⋮	⋮	⋮	⋮		⋮	⋮	⋮
	q	${\widehat{P}}_{q1}$	${\widehat{P}}_{q2}$	…	${\widehat{P}}_{qj}$	…	${\widehat{P}}_{qq}$	${\widehat{P}}_{q+}$
	Total	${\widehat{P}}_{+1}$	${\widehat{P}}_{+2}$	…	${\widehat{P}}_{+j}$	…	${\widehat{P}}_{+q}$	1

Table 9. Adjusted error matrix of estimated area proportions.

			Reference
		1	2	…	j	…	q	Total
Map	1	${\widehat{P}}_{11}$	${\widehat{P}}_{12}$	…	${\widehat{P}}_{1j}$		${\widehat{P}}_{1q}$	${\widehat{P}}_{1+}$
	2	${\widehat{P}}_{21}$	${\widehat{P}}_{22}$	…		…		${\widehat{P}}_{2+}$
	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
	i	${\widehat{P}}_{i1}$	${\widehat{P}}_{i2}$	…	${\widehat{P}}_{ij}$	…	${\widehat{P}}_{iq}$	${\widehat{P}}_{i+}$
	⋮	⋮	⋮	⋮		⋮	⋮	⋮
	q	${\widehat{P}}_{q1}$	${\widehat{P}}_{q2}$	…	${\widehat{P}}_{qj}$	…	${\widehat{P}}_{qq}$	${\widehat{P}}_{q+}$
	Total	${\widehat{P}}_{+1}$	${\widehat{P}}_{+2}$	…	${\widehat{P}}_{+j}$	…	${\widehat{P}}_{+q}$	1

Based on this new error matrix, accuracies (overall, user’s and producer’s) were calculated. The overall accuracy ($\widehat{O})$ indicates the overall proportion of area correctly classified. It is the sum of ${\widehat{p}}_{ii}$ of the adjusted error matrix diagonal.

$$\widehat{O}={\displaystyle \sum }_{i=1}^q{\widehat{P}}_{ii}$$

(2)

User’s accuracy (${\widehat{U}}_i$) of class i is the proportion of the area mapped as class i that has reference class i, and producer’s accuracy (${\widehat{P}}_j$) of class j is the proportion of the area of reference class j that is mapped as class j. Both accuracies were computed according to Equations (3) and (4), respectively.

$${\widehat{U}}_i=\frac{{\widehat{P}}_{ii}}{{\widehat{P}}_{i+}}$$

(3)

$${\widehat{P}}_j=\frac{{\widehat{P}}_{jj}}{{\widehat{P}}_{+j}}$$

(4)

Classification types performance

McNemar’s test was used to evaluate the performance of multi-temporal and mono-temporal classification as well as their integration with ancillary data. It is a non-parametric test that is more precise and sensitive than the Kappa z-test [46]. McNemar’s test evaluates the difference between paired proportions. Here, the test is based on a pair of confusion matrices of correctly and wrongly classified reference samples. It produces a chi-square (${\chi }^2$) statistics that is computed by Equation (5). The difference between two classification approaches is significant when the p-value is less than 0.05.

$${\chi }^2=\frac{{(f_{12}-f_{21})}^2}{\left(f_{12}+f_{21}\right)}$$

(5)

where $f_{12}$ indicates the number of cases that are wrongly classified by approach 2, but correctly classified by approach 1, and $f_{21}$ is the number of cases that are correctly classified by approach 2, but wrongly classified by approach 1.

Area estimates and uncertainty

In accordance with the systematic approach suggested by Olofsson et al. [42], the adjusted error matrix was used to compute an area estimator based on the proportion of the area of category j. Equation (6) gives the area of category j (${\widehat{A}}_j$).

$${\widehat{A}}_j= A_{tot}\times {\widehat{P}}_{+j}$$

(6)

where $A_{tot}$ is the total area, and

$${\widehat{P}}_{+j}={\displaystyle {\sum }_{i=1}^q}{W}_i\frac{n_{ij}}{n_{i+}}$$

(7)

This area estimator is an error-adjusted estimator of area that includes the area of map omission error of category j and removes the area of map commission error [42]. Its standard error $S\left({\widehat{A}}_j\right)$was computed as follows.

$$S\left({\widehat{A}}_j\right)=A_{tot} \times S\left({\widehat{P}}_{+j}\right)$$

(8)

where the standard error for the stratified estimator of proportion of area$S\left({\widehat{P}}_{+j}\right)$ is computed as:

$$S\left({\widehat{P}}_{+j}\right)=\sqrt{{\displaystyle \sum }_{i=1}^q\frac{W_i{\widehat{P}}_{ik-}{\widehat{P}}_{ik}^2}{n_{i+}-1}}$$

(9)

The use of ${\widehat{P}}_{+j}$ (estimated from the reference samples) instead of ${\widehat{P}}_{i+}$(map areas) is because it allows the assessment of uncertainty of the area estimates in the form of sampling variability that can be computed as confidence interval. For ${\widehat{A}}_j$the approximate 95% confidence interval (CI) was derived as:

$$CI={\widehat{A}}_j\pm z\times S\left({\widehat{A}}_j\right)$$

(10)

where z corresponds to the percentile from the curve of the standard normal distribution, and for 95% confidence, z = 1.96.

3. Results

3.1. Suitable Period for Mono-Temporal LULC Classification

In order to detect the most adequate period to classify the observed LULC, the year 2011 was selected because it had images from dry season (March), early rainy season (May and June), mid rainy season (July) and late rainy season (October), which cover the entire vegetation period. The results of the mono-temporal RF classifications are shown in Figure 2, and it can be observed that there is an increase in overall, average user’s and average producer’s accuracies from March to October. Table 10 presents the pixels count error matrices obtained from the five mono-temporal classifications. Major confusions were recorded, for instance, between agricultural area and natural vegetation, agricultural area and bare surface and among natural vegetation types. The use of October image for classification outperformed all other mono-temporal attempts, because it reduced confusion between natural vegetation classes (Woodland and Mixed vegetation) and between agricultural area and natural vegetation, among others. The matrices clarify that at the end of the dry period (March) the class agriculture area seems to spectrally resemble other classes, which explains the aforementioned low accuracies obtained for that period. The main finding is that the LULC classification of the late rainy season images (e.g., October) performed better than those applied to images of the mid rainy, early rainy and dry seasons. Therefore, the mono-temporal image of October served as benchmark for assessing the effect of multi-temporal images and ancillary data on LULC classification.

Figure 2. Overall, user’s and producer’s accuracies of the mono-temporal classifications of 2011 derived from adjusted error matrix. User’s and producer’s accuracies are averages of the class-wise assessments.

Figure 2. Overall, user’s and producer’s accuracies of the mono-temporal classifications of 2011 derived from adjusted error matrix. User’s and producer’s accuracies are averages of the class-wise assessments.

3.2. LULC Classification Accuracies According to Images Combinations

The results presented in Table 11 show that mono-temporal classification (of the best performing data set, October) had the lowest accuracy among the four combinations of images. In contrast to the optimal mono-temporal approach, multi-temporal classification improved the overall accuracy from 88% to 94%, the average user’s from 87% to 93%, and the average producer’s from 86% to 91%. The incorporation of ancillary data has also positively influenced the LULC discrimination. Indeed, their integration with mono-temporal image (October) yielded an overall accuracy comparable to the use of multi-temporal data (94%). There was also an improvement in average user’s accuracy (87% to 93%) and average producer’s accuracy (86% to 91%). However, the best performance in terms of the accuracy assessment was recorded by the combination of multi-temporal images and ancillary data, which produced the highest overall, average user’s and average producer’s accuracies with values of 95%, 95% and 92%, respectively. This superiority of multi-temporal classification and the addition of ancillary data over mono-temporal classification is supported by the results of McNemar’s test (Table 12). The test showed significant difference (p < 0.05) in LULC classifications produced by multi-temporal, mono-temporal plus ancillary data and multi-temporal plus ancillary data classifications over mono-temporal data. However, no significance difference was found between multi-temporal, mono-temporal plus ancillary data and multi-temporal plus ancillary data classifications.

The positive impact of multi-temporal images and ancillary data was also perceptible on each LULC class. The user’s and producer’s accuracies of each class, according to Table 13 and Table 14, indicate that the discrimination of the individual LULC classes was improved with the use of multi-temporal images and the addition of ancillary data.

Table 10. Error matrices of mono-temporal LULC classifications of 2011, the values in all the tables refer to the numbers of pixels. 1: Water; 2: Woodland; 3: Bare surface; 4: Mixed vegetation; 5: Agricultural area.

March
	1	2	3	4	5	Total
1	15	0	2	0	25	42
2	0	79	0	5	16	100
3	0	4	35	7	11	57
4	1	16	12	34	14	77
5	0	0	0	5	90	95
Total	16	99	49	51	156	371
May
	1	2	3	4	5	Total
1	13	0	16	0	13	42
2	3	91	0	3	3	100
3	0	2	44	5	8	59
4	0	20	4	39	12	75
5	1	4	2	5	83	95
Total	17	117	66	52	119	371
June
	1	2	3	4	5	Total
1	37	0	0	0	5	42
2	4	91	2	1	2	100
3	0	5	46	3	5	59
4	0	20	9	36	10	75
5	0	3	1	5	86	95
Total	41	119	58	45	108	371
July
	1	2	3	4	5	Total
1	42	0	0	0	0	42
2	6	87	1	3	3	100
3	0	6	46	2	5	59
4	0	9	8	43	15	75
5	0	2	1	0	92	95
Total	48	104	56	48	115	371
October
	1	2	3	4	5	Total
1	40	2	0	0	0	42
2	0	90	0	6	4	100
3	0	0	48	6	5	59
4	0	5	4	57	9	75
5	0	0	0	5	90	95
Total	40	97	52	74	108	371

Table 10. Error matrices of mono-temporal LULC classifications of 2011, the values in all the tables refer to the numbers of pixels. 1: Water; 2: Woodland; 3: Bare surface; 4: Mixed vegetation; 5: Agricultural area.

March
	1	2	3	4	5	Total
1	15	0	2	0	25	42
2	0	79	0	5	16	100
3	0	4	35	7	11	57
4	1	16	12	34	14	77
5	0	0	0	5	90	95
Total	16	99	49	51	156	371
May
	1	2	3	4	5	Total
1	13	0	16	0	13	42
2	3	91	0	3	3	100
3	0	2	44	5	8	59
4	0	20	4	39	12	75
5	1	4	2	5	83	95
Total	17	117	66	52	119	371
June
	1	2	3	4	5	Total
1	37	0	0	0	5	42
2	4	91	2	1	2	100
3	0	5	46	3	5	59
4	0	20	9	36	10	75
5	0	3	1	5	86	95
Total	41	119	58	45	108	371
July
	1	2	3	4	5	Total
1	42	0	0	0	0	42
2	6	87	1	3	3	100
3	0	6	46	2	5	59
4	0	9	8	43	15	75
5	0	2	1	0	92	95
Total	48	104	56	48	115	371
October
	1	2	3	4	5	Total
1	40	2	0	0	0	42
2	0	90	0	6	4	100
3	0	0	48	6	5	59
4	0	5	4	57	9	75
5	0	0	0	5	90	95
Total	40	97	52	74	108	371

Table 11. Classification accuracy in 2011 (%) according to images combinations.

Images Combinations	Overall Accuracy	Av. User’s acc.	Av. Producer’s Acc.
Mono-temporal	88	87	86
Mono-temporal plus ancillary	94	93	91
Multi-temporal	94	93	91
Multi-temporal plus ancillary data	95	95	92

Table 11. Classification accuracy in 2011 (%) according to images combinations.

Images Combinations	Overall Accuracy	Av. User’s acc.	Av. Producer’s Acc.
Mono-temporal	88	87	86
Mono-temporal plus ancillary	94	93	91
Multi-temporal	94	93	91
Multi-temporal plus ancillary data	95	95	92

Table 12. McNemar’s test results between image combinations.

	Mono-Temporal
	$f_{12}$	$f_{21}$	Chi-square	p-value
Multi-temporal	4	22	12.5	0.0004
Mono-temporal plus ancillary data	6	23	10	0.001
Multi-temporal plus ancillary data	1	19	16.2	0.00005
	Mono-temporal plus ancillary
	$f_{12}$	$f_{21}$	Chi-square	p-value
Multi-temporal	6	7	0.08	0.8
Multi-temporal plus ancillary data	10	5	1.7	0.2
Multi-temporal plus ancillary data	Multi-temporal
	$f_{12}$	$f_{21}$	Chi-square	p-value
	8	3	2.3	0.1

Table 12. McNemar’s test results between image combinations.

	Mono-Temporal
	$f_{12}$	$f_{21}$	Chi-square	p-value
Multi-temporal	4	22	12.5	0.0004
Mono-temporal plus ancillary data	6	23	10	0.001
Multi-temporal plus ancillary data	1	19	16.2	0.00005
	Mono-temporal plus ancillary
	$f_{12}$	$f_{21}$	Chi-square	p-value
Multi-temporal	6	7	0.08	0.8
Multi-temporal plus ancillary data	10	5	1.7	0.2
Multi-temporal plus ancillary data	Multi-temporal
	$f_{12}$	$f_{21}$	Chi-square	p-value
	8	3	2.3	0.1

3.3. Contribution of Remotely Sensed Bands and Ancillary Data to LULC Classification

The contribution of the remotely sensed bands and ancillary data to LULC classification in the study area is shown in Figure 3, which highlights the mean decrease accuracy (MDA) score of each variable. This was computed by RF algorithm based on the classification of mono-temporal image plus ancillary data. In general, near infrared, middle infrared and elevation were the three most important variables for LULC classification. They are followed by geomorphology, red, green and blue bands, whereas slope and soil types were the least important variables. Among the four environmental data used in this research, elevation and geomorphology have more contributed to LULC discrimination compared to soil types and slope.

Table 13. LULC user’s accuracy (%) per image combination in 2011.

Classification Approach	Water	Woodland	Bare Surface	Mixed Vegetation	Agricultural Area
Mono-temporal	95	90	81	76	94
Mono-temporal plus ancillary	95	98	90	88	95
Multi-temporal	95	96	90	89	97
Multi-temporal plus ancillary	98	99	93	88	97

Table 13. LULC user’s accuracy (%) per image combination in 2011.

Classification Approach	Water	Woodland	Bare Surface	Mixed Vegetation	Agricultural Area
Mono-temporal	95	90	81	76	94
Mono-temporal plus ancillary	95	98	90	88	95
Multi-temporal	95	96	90	89	97
Multi-temporal plus ancillary	98	99	93	88	97

Table 14. LULC producer’s accuracy (%) per image combination in 2011.

Classification Approach	Water	Woodland	Bare Surface	Mixed Vegetation	Agricultural. Area
Mono-temporal	100	96	66	81	85
Mono-temporal plus ancillary	100	100	69	92	93
Multi-temporal	100	100	71	92	93
Multi-temporal plus ancillary	100	98	69	94	97

Table 14. LULC producer’s accuracy (%) per image combination in 2011.

Classification Approach	Water	Woodland	Bare Surface	Mixed Vegetation	Agricultural. Area
Mono-temporal	100	96	66	81	85
Mono-temporal plus ancillary	100	100	69	92	93
Multi-temporal	100	100	71	92	93
Multi-temporal plus ancillary	100	98	69	94	97

Figure 3. Remotely sensed bands and ancillary data contributions to LULC classification based on mean decrease accuracy (MDA) score of RF mono-temporal image plus ancillary data classification.

Figure 3. Remotely sensed bands and ancillary data contributions to LULC classification based on mean decrease accuracy (MDA) score of RF mono-temporal image plus ancillary data classification.

3.4. The Dynamics of LULC in the Study Area over the Years 1999, 2006 and 2011

In accordance with the results of the 2011 classifications, the combination of multi-temporal images and ancillary data was used to classify the images of 1999 and 2006. The monthly Landsat images of 1999 and 2006, specified in Table 1, were combined with ancillary data (Table 2), and the classification yielded 94% (93%) as overall accuracy for the LULC map of 1999 (2006). The distribution of each LULC area in 1999, 2006 and 2011 is presented in Table 15, Table 16 and Table 17, respectively. The mapped area is slightly different from the estimated area for all LULC classes. In all cases, it could be found within the confidence interval (95%) around the estimated area indicating a high reliability of the produced maps [42].

Table 15. Proportion of LULC types in 1999.

LULC	Mapped Area %	Estimated Area %	Confidence Interval %
Water	1.7	1.7	± 0.1
Woodland	41.5	40.1	± 2.0
Bare surface	1.8	2.4	± 1.0
Mixed vegetation	33.2	33.2	± 2.3
Agricultural area	21.8	22.6	± 2.3
Total	100	100

Table 15. Proportion of LULC types in 1999.

LULC	Mapped Area %	Estimated Area %	Confidence Interval %
Water	1.7	1.7	± 0.1
Woodland	41.5	40.1	± 2.0
Bare surface	1.8	2.4	± 1.0
Mixed vegetation	33.2	33.2	± 2.3
Agricultural area	21.8	22.6	± 2.3
Total	100	100

Table 16. Proportion of LULC types in 2006.

LULC	Mapped Area %	Estimated Area %	Confidence Interval %
Water	0.7	0.7	± 0.04
Woodland	39	38.1	± 1.3
Bare surface	1.8	2.9	± 1.3
Mixed vegetation	32.5	32.3	± 1.7
Agricultural area	26	26	±2 .6
Total	100	100

Table 16. Proportion of LULC types in 2006.

LULC	Mapped Area %	Estimated Area %	Confidence Interval %
Water	0.7	0.7	± 0.04
Woodland	39	38.1	± 1.3
Bare surface	1.8	2.9	± 1.3
Mixed vegetation	32.5	32.3	± 1.7
Agricultural area	26	26	±2 .6
Total	100	100

Table 17. Proportion of LULC types in 2011.

LULC	Mapped Area%	Estimated Area %	Confidence Interval %
Water	0.2	0.2	± 0.01
Woodland	35.3	36.4	± 2.1
Bare surface	2.9	3.3	± 1.4
Mixed vegetation	31.6	30.4	± 1.3
Agricultural area	30	29.7	± 1.6
Total	100	100

Table 17. Proportion of LULC types in 2011.

LULC	Mapped Area%	Estimated Area %	Confidence Interval %
Water	0.2	0.2	± 0.01
Woodland	35.3	36.4	± 2.1
Bare surface	2.9	3.3	± 1.4
Mixed vegetation	31.6	30.4	± 1.3
Agricultural area	30	29.7	± 1.6
Total	100	100

The results revealed that agricultural area and bare surface have increased in the study area at the expense of woodland and mixed vegetation, which decreased over the years. For instance, the proportion of agricultural area increased from 22.6% ± 2.3% in 1999 to 26% ± 2.6% in 2006 and further to 29.7% ± 1.6% in 2011. These dynamics are well captured in Figure 4, which highlights the expansion of agricultural area (in yellow) across the study area. Contrary, the areas of woodland and mixed vegetation, as observed in 1999, decreased in 2006 and 2011 (Table 15, Table 16 and Table 17). A decreasing trend was also observed for areas covered by water.

Figure 4. LULC spatial distribution in the study area in 1999, 2006 and 2011.

Figure 4. LULC spatial distribution in the study area in 1999, 2006 and 2011.

3.5. LULCC in the Study Area between 1999 and 2011

The adjusted error matrix of the LULCC map derived between 1999 and 2011 is given by Table 18. The accuracy assessment returned an overall accuracy of 92%, while the user’s accuracy ranged from 86% to 95% and the producer’s accuracy from 84% to 95%.

Table 18. Adjusted error matrix of LULCC map between 1999 and 2011. SNV: stable natural vegetation; NVL: natural vegetation loss; SNNV: stable non-natural vegetation; OC: other change.

	SNV	NVL	SNNV	OC	Total	User’s (%)
SNV	0.547	0.014	0	0.014	0.575	95
NVL	0.016	0.151	0.003	0.003	0.173	87
SNNV	0.004	0.01	0.129	0	0.143	90
OC	0.007	0.004	0.004	0.094	0.109	86
Total	0.574	0.179	0.136	0.111	1
Producer’s (%)	95	84	95	85		92

Table 18. Adjusted error matrix of LULCC map between 1999 and 2011. SNV: stable natural vegetation; NVL: natural vegetation loss; SNNV: stable non-natural vegetation; OC: other change.

	SNV	NVL	SNNV	OC	Total	User’s (%)
SNV	0.547	0.014	0	0.014	0.575	95
NVL	0.016	0.151	0.003	0.003	0.173	87
SNNV	0.004	0.01	0.129	0	0.143	90
OC	0.007	0.004	0.004	0.094	0.109	86
Total	0.574	0.179	0.136	0.111	1
Producer’s (%)	95	84	95	85		92

Table 19 presents the estimated and the mapped LULCC areas. In this case also, the map area of each LULCC class fell within the computed confidence interval (95%). Between 1999 and 2011, 57.4% ± 2.7% of the study area was dominated by stable natural vegetation. Area under loss of natural vegetation represented 17.9% ± 2.5%. Stable non-natural vegetation and other change covered 13.6% ± 1.5% and 11.1% ± 1.9% of the study area, respectively. The distribution of the LULCC classes in the study area is shown in Figure 5 where most of the losses of natural vegetation are located in areas dominated by agriculture.

Table 19. LULCC area between 1999 and 2011.

	Map Area	Estimated Area
LULCC classes	%	%	Conf. interval (%)
Stable natural vegetation	57.5	57.4	± 2.7
Natural vegetation loss	17.2	17.9	± 2.5
Stable non-natural vegetation	14.3	13.6	± 1.5
Other change	11	11.1	± 1.9
Total	100	100

Table 19. LULCC area between 1999 and 2011.

	Map Area	Estimated Area
LULCC classes	%	%	Conf. interval (%)
Stable natural vegetation	57.5	57.4	± 2.7
Natural vegetation loss	17.2	17.9	± 2.5
Stable non-natural vegetation	14.3	13.6	± 1.5
Other change	11	11.1	± 1.9
Total	100	100

Figure 5. Distribution of LULCC in the study area between 1999 and 2011.

Figure 5. Distribution of LULCC in the study area between 1999 and 2011.

The main LULC conversions observed between 1999 and 2011 (Table 20) consisted of changes from woodland to agricultural area (occurred on 8.6% of the study area), woodland to mixed vegetation (8.4%), mixed vegetation to woodland (7.2%) and mixed vegetation to agricultural area (6.8%). Stable LULC are dominated by woodland, mixed vegetation and agricultural area, which covered 24.1%, 17.8% and 13.8% of the area, respectively. Although there was an important loss of natural vegetation (17.9% ± 2.5%) in the study area, the transfer matrix highlighted that also regrowth of natural vegetation was recorded. As example, 2.7%, 4.5% and 0.8% of the study area were converted from agricultural area to woodland, agricultural area to mixed vegetation and bare surface to mixed vegetation, respectively. Those conversions could be cropland left under fallow, and potentially the result of reforestation and afforestation campaigns that are frequent in this area.

Table 20. LULC transfer matrix between 1999 and 2011 expressed as percentage of area.

		2011
		Water	Woodland	Bare Surfaces	Mixed Veg.	Agric. Area	Area 1999
1999	Water	0.2	1.1	0	0.1	0.3	1.7
	Woodland	0	24.1	0.4	8.4	8.6	41.5
	Bare surfaces	0	0.2	0.3	0.8	0.5	1.8
	Mixed veg.	0.01	7.2	1.4	17.8	6.8	33.2
	Agric. Area	0	2.7	0.8	4.5	13.8	21.8
	Area 2011	0.2	35.3	2.9	31.6	30	100

Table 20. LULC transfer matrix between 1999 and 2011 expressed as percentage of area.

		2011
		Water	Woodland	Bare Surfaces	Mixed Veg.	Agric. Area	Area 1999
1999	Water	0.2	1.1	0	0.1	0.3	1.7
	Woodland	0	24.1	0.4	8.4	8.6	41.5
	Bare surfaces	0	0.2	0.3	0.8	0.5	1.8
	Mixed veg.	0.01	7.2	1.4	17.8	6.8	33.2
	Agric. Area	0	2.7	0.8	4.5	13.8	21.8
	Area 2011	0.2	35.3	2.9	31.6	30	100

4. Discussion

4.1. LULC Classification

The use of late rainy season image (e.g., October image) was found more suitable than single image from other seasons (mid rainy, early rainy and dry seasons) for LULC classification in the Sudanese savannah of southwest BF. This is due the fact that late rainy season image provides better spectral information to separate confusing LULC classes such as agricultural area and natural vegetation. Medicinal, sacred, and fruit trees maintained in croplands creating intermixed areas (Figure 6) can be seen as one major reason for this spectral confusion between these two classes in the study region.

The result obtained for the Sudanese savannah are in agreement with the study by Lunetta and Balogh [10] and Key et al. [14], who noticed an improvement in LULC classification accuracy when they used multi-temporal images as compared to mono-temporal. However, findings in this study are in contrast with those reported in Langley et al. [19]. This is because the study by Langley et al. [19] focused on areas dominated by grassland located in the Jornada del Muerto plain of southern New Mexico, while the present work investigated an area dominated by woodland and mixed vegetation with many patterns of agricultural areas. These findings therefore highlight that differences in environment and LULC types might play a critical role.

Multi-temporal classification enables the analysis of images acquired at different phenological stages, which adds useful information for the classification, and in turn permits a better class discrimination. It has to be noted that this research did not assess different multi-temporal image compositions, like various combinations of months, which could be helpful for the choice of adequate temporal images [47,48]. However, the mono-temporal results strongly indicate that for multi-temporal mapping, the inclusion of late rainy season images is of relevance for accurate LULC mapping in southwest Burkina Faso. The remaining challenge is that frequent cloud cover in the region may reduce the availability of data in such important temporal windows suitable for classification [49].

In the face of frequently disturbed atmospheric conditions and poor availability of cloud-free satellite data in the study region, it is of high interest that ancillary data can significantly improve LULC classification to almost the same accuracy level as multi-temporal classifications. In the study area for instance, natural vegetation distribution (e.g., Woodland) is influenced by topography [13]. This explains the usefulness of including environmental variables into LULC discrimination especially elevation and geomorphology units that added most inputs to the classification after the remotely sensed infrared bands. These results confirm similar conclusions of previous studies, which utilized the environmental settings as add-on for an improved classification. Sesnie et al. [15] found the overall accuracy increasing from 82.4% to 87.4% with the addition of terrain variables (e.g., elevation and slope) in the classification of nine land cover types in Nicaraguan tropical dry forest. For monitoring of land-cover change in San Diego County, Rogan et al. [11] noted that ancillary variables (elevation, fire history, and Slope) contributed 15% to the overall accuracy of land-cover change mapping.

The most accurate classification was achieved when multi-temporal images and ancillary data were combined. The accuracy of this approach significantly differed from mono-temporal classification (as shown by the McNemar tests at the 95% confidence level), but not from classification results using a combination of specific mono-temporal image with ancillary data. This again shows that combining mono-temporal satellite data with information about environmental settings can be assessed to be a suitable alternative even in the heterogeneous savannah of southwest Burkina Faso.

The ancillary data showed positive impact on the classification results, but it needs to be stated that geomorphology and soil data, were produced at a scale of 1:500,000, which is at a coarse spatial resolution in comparison to 30 m Landsat pixels. The use of higher spatial resolution data matching the satellite data pixel could possibly improve the accuracy of LULC classifications. In general, LULC/LULCC maps presented in this study can be seen as a valuable input for training and validating global or regional LULC maps, which in turn will be used for regional assessments e.g. of the carbon budget [50]. Local maps are of particular value within the very complex landscape of West Africa, which extremely challenges the accurateness of existing regional to global maps [51]. Assessing the accuracy of regional maps requires an estimation of error propagation and hence documented quality information about the local input maps. The maps presented here, have such quality attributes even though it must be stated that validation was carried out on only 28.5% of the study area. For applying the validation throughout the study area, the availability of more high-resolution data sets would be desirable. With the increasing number of high-resolution data sources (e.g., QuickBird, RapidEye and Ikonos), future studies will potentially consider larger portions of their study area for accuracy assessment, which may give more confidence to the results. However, for instance validation activities of other LULC mapping studies relying on field observations can usually cover only those parts of the study areas which can be accessed either by car or by bike (e.g., [13]), which also questions the spatial distribution of the samples and hence the overall validity of the accuracy assessment. However, a validated LULC map based on some portions of the study area remains more reliable than a non-validated LULC map as observed in previous studies conducted in West Africa (e.g., [21]).

Figure 6. Harvested crop intermixed with trees in the study area.

Figure 6. Harvested crop intermixed with trees in the study area.

4.2. LULCC in the Study Area

Availability of adequate reference data provided a unique baseline for conducting LULCC map accuracy assessment in southwestern Burkina Faso. The observed changes of LULC in the study area were characterized by the increase in agricultural area and bare surface at the expense of natural vegetation (mixed vegetation and woodland). These dynamics attest to the ongoing deforestation in the southwest of Burkina Faso. Human activities are the main drivers of these changes in LULC, because in this region, where farming is characterized by low inputs such as fertilizers, farms are usually expanded and scattered across the landscape to increase yield [52]. Such a practice combined with population growth could explain the conversion of relatively huge areas of natural vegetation to agricultural area, as was observed in the period from 1999 to 2011. This conversion was also noticed elsewhere in the Sudanese savannah (e.g., [17,18,22]). Besides, in this region, fallowing practice has been reduced [53] leading to excessive cultivation of the same land over many years resulting in soil fertility loss, which favors the augmentation of bare surface. Furthermore, expansion of human settlement across the study area, as in Diebougou, Dano, Fara and Dissin, could also be a contributory factor to increasing bare surface at the expense of natural vegetation. Apart from agriculture and urban growth, other anthropogenic activities, such as wood harvesting for charcoal production, bushfire and mining, contribute to deforestation in the study area (Figure 7).

The loss of natural vegetation could have also been exacerbated by climate variability, which acts as a catalyst to the anthropogenic pressure. Indeed, in West Africa, IPCC [54] reported that climate change has increased rainfall variability and the recurrence of extreme events such as droughts. Figure 8 shows the temporal variability of rainfall (1981–2012) in the study area where the red bars indicates deficit of water and the blue represents the surplus of water. Three out of four stations (Fara, Diebougou and Dissin) highlighted that the years 1999 and 2006 were more humid than 2011. In addition, the year 1999 was characterized by floods in the southwest of Burkina Faso [55]. This further buttress the observed increase in bare surface in 2011 compared to 2006 and 1999. In fact, in this region, a decrease in rainfall leads to an increase in barren land area [12], mainly due to the close relationship between rainfall variability and vegetation dynamics [56].

Figure 7. (a) Trees cut for new cropland; and (b) fuel wood collected in an area devastated by bushfire.

Figure 7. (a) Trees cut for new cropland; and (b) fuel wood collected in an area devastated by bushfire.

Figure 8. Rainfall variability in the study area (1981–2012) expressed as standardized precipitation index (SPI) distribution (Rainfall data collected from the national direction of meteorology of Burkina Faso).

Figure 8. Rainfall variability in the study area (1981–2012) expressed as standardized precipitation index (SPI) distribution (Rainfall data collected from the national direction of meteorology of Burkina Faso).

5. Conclusions

This study, carried out in the Sudanese savannah of southwest Burkina Faso, compared mono-temporal and multi-temporal classification as well as their combination with ancillary data to determine the best classification method in order to produce an accurate post-classification LULCC map. In the study region, late rainy season images, used as input to a RF classifier algorithm, produced higher overall accuracy (88%) than images from other periods. Multi-temporal classification significantly improved mono-temporal classification of LULC in terms of overall and class-wise accuracies. It was found that the overall accuracy increased from 88% to 94%, the average user’s from 87% to 93%, and the average producer’s from 86% to 91%. The inclusion of ancillary data was found beneficial for multi-temporal classifications and also enhances mono-temporal classification accuracies to the level of multi-temporal classifications (94%). Among the ancillary data (elevation, geomorphology, slope and soil types), elevation and geomorphology were found to be the most contributors to LULC classification in the study area. The results suggest that, at least for the heterogeneous Sudanese savannah of southwest Burkina Faso, the inclusion of ancillary data reduces the data requirements for accurate LULC/LULCC maps, where atmospheric conditions are the limiting factors for the availability of multi-temporal data. However, at least a cloud-free satellite data set from late rainy season period should be available for successful mapping. Of course, this method needs to be tested in other regions not only within the Sudanese savannah for a more generalized assessment of its performance and spatial transferability to LULC/LULCC maps in West Africa and beyond.

The study also found unsustainable LULC dynamics within the savannah of southwest Burkina Faso, which is characterized by increasing agricultural area and bare surface at the detriment of woodland and mixed vegetation. For instance, between 1999 and 2011, agricultural area increased from 22.6% ± 2.3% to 29.7% ± 1.6%, and contrary, woodland and mixed vegetation decreased from 40.1% ± 2.0% to 36.4% ± 2.1% and from 33.2% ± 2.3% to 30.4% ± 1.3%, respectively. The LULCC map produced an overall accuracy of 92% and the loss of natural vegetation was estimated to be 17.9% ± 2.5% in the period 1999–2011. Human activities are supposed to be the main drivers of the observed changes in LULC, but rainfall variability most probably also contributed to the observed LULCC. The loss of natural vegetation increases carbon emission in the atmosphere and could accelerate the process of land degradation in the study area. Of course, bi-temporal comparisons cannot attribute all observed alterations to land use change but can be seen as indicators of ongoing developments. Trend analysis of vegetation cover such as those obtained from moderate resolution time series (e.g., Moderate-Resolution Imaging Spectroradiometer, MODIS) can substantiate the hints of the high-resolution change map and allows for better assessment of intra-annual variability of vegetation cover, e.g., due to climate variability within the change analysis. However, the LULCC maps presented in this study can successfully support monitoring of ongoing reforestation activities. It can also be used to locate areas where adequate measures, for more efficient and sustainable land use management, can be introduced to safeguard the natural environment and livelihood of local population in the southwest Burkina Faso.