A Forest Monitoring System for Tanzania

Abstract

Tropical forests provide essential ecosystem services related to human livelihoods. However, the distribution and condition of tropical forests are under significant pressure, causing shrinkage and risking biodiversity loss across the tropics. Tanzania is currently undergoing significant forest cover changes, but monitoring is limited, in part due to a lack of remote sensing knowledge, tools and methods. This study has demonstrated a comprehensive approach to creating a national-scale forest monitoring system using Earth Observation data to inform decision making, policy formulation, and combat biodiversity loss. A systematically wall-to-wall forest baseline was created for 2018 through the application of Landsat 8 imagery. The classification was developed using the extreme gradient boosting (XGBoost) machine-learning algorithm, and achieved an accuracy of 89% and identified 45.76% of the country’s area to be covered with forest. Of those forested areas, 45% was found within nationally protected areas. Utilising an innovative methodology based on a forest habitat suitability analysis, the forest baseline was classified into forest types, with an overall accuracy of 85%. Woodlands (open and closed) were found to make up 79% of Tanzania’s forests. To map changes in forest extent, an automated system for downloading and processing of the Landsat imagery was used along with the XGBoost classifiers trained to define the national forest extent, where Landsat 8 scenes were individually downloaded and processed and the identified changes summarised on an annual basis. Forest loss identified for 2019 was found to be 157,204 hectares, with an overall accuracy of 82%. These forest losses within Tanzania have already triggered ecological problems and alterations in ecosystem types and species loss. Therefore, a forest monitoring system, such as the one presented in this study, will enhance conservation programmes and support efforts to save the last remnants of Tanzania’s pristine forests.

1. Introduction

Tropical forests contain the most distinct and complex biome on Earth, with unique plant species of high economic value, and support habitat for many animal species. They provide numerous valuable ecosystem services while also aiding the mitigation of climate change [1,2,3,4]. Africa is home to some of the world’s most magnificent tropical forests. With more more than 60 million people dwelling within or near these forests, they are relied upon for many ecosystem systems, with livelihoods dependent on them for providing food, medicinal plants, fuel, fibres, and non-timber forest products. At the same time, they are also important for societal and cultural purposes [5]. Despite the importance of these forests, the spatial and temporal variation in tropical forest cover has raised the dynamic state and increased forest cover decline. However, throughout many parts of Africa, forest cover change remains poorly understood due to field-based monitoring challenges and a lack of remote sensing studies. Government policies have often failed to prevent illegal forestry activities due to the absence of defined tenure, which has also increased over-utilisation. This is also linked to the lack of national capacity for monitoring and reporting deforestation, especially in many sub-Saharan Africa countries [3] and hence remains a barrier.

Therefore, the loss of tropical forests has been of international concern [6] in various Multilateral Environmental Agreements (MEA) such as the Convention on Biological Diversity (CBD) and the United Nations Framework Convention on Climate Change (UNFCCC) [7], as it is linked to increases in erosion, run-off and flooding, CO₂ concentration, climate change, and biodiversity loss [4,8,9,10]. This demonstrates the need for more conservation of tropical forests, particularly threats from the expansion of agriculture, urban expansion and other forms of deforestation in contemporary society. Accurate information on the extent and extinction of tropical forests is necessary as the yardstick for policy and decision makers.

The application of Earth Observation data and growing computational power can help characterise tropical forests and monitoring at the landscape level to address inadequate information [11]. The provision of satellite-based data has come a long way since the 1970s, with more missions to be launched over the coming years [12] that will support the mapping of changes in tropical forest areas with broad geographic coverage. However, efforts to provide ways to operationalise Earth Observation data for tropical forest monitoring have, to date, been heavily biased toward particular regions, such as South America [11], while Africa, apart from the Congo Basin, has been severely understudied. For example, from 1995 to 2003, approximately two-thirds of studies focused on the Amazon Basin, 18% focused on central Africa, and 17% on Southeast Asia [13]. Therefore, expanding tropical forest monitoring beyond these areas is required to capitalise on Earth Observation data availability with extensive area coverage. Perhaps most importantly, attempts to combine societal constraints while contributing to an increase in scientific knowledge and provide up-to-date and appropriate maps and statistics of tropical forest conservation activities on the ground are needed.

The situation is more complicated in developing countries such as the United Republic of Tanzania (hereafter referred to as Tanzania), where many households are highly dependent on forest resources. Despite the value of forests, changes in use patterns pose a noticeable threat to forest resource sustainability on both socio-economic and ecological functioning. This has led to the increased scarcity of forest resources, further aggravated by the continuing high deforestation rate [14] and associated with a change in climatic conditions experienced in many parts of Tanzania. Therefore, the need for data and information on the state of Tanzania’s forest resources is of increasing importance. Yet, Tanzania’s forest resource status and trends are mostly unknown, with current data being fragmented and outdated [15].

The lack of institutional capacity has largely constrained data reliability on Tanzanian forest resources, with inadequate national-wide forest monitoring coverage using an Earth Observation-based system. Using freely accessible satellite data and advanced remote sensing methods can provide a cost-effective and timely approach to achieving systematic wall-to-wall information for Tanzania’s forests. This information is required to support national policy processes aimed at improving sustainable forest management—at the same time addressing the issues of Reduced Emissions from Deforestation and Forest Degradation (REDD+) and Green House Gas (GHG) as international reporting obligations. Also the 2030 Agenda for Sustainable Development Goals (SDGs) for enhancing life on land (goal 15) and its targets through combating deforestation [16].

Large-scale forest monitoring by remote sensing has been reported from global forest change studies and provides initial forest loss estimates. Hansen et al. [14] provides a global forest loss estimate from 2001 to 2019 with a spatial resolution of 30 m through the publically available Global Forest Watch dataset (version 1.7). As a global and freely available dataset, it offers extensive forest change information and the accuracy of this dataset should be assessed for Tanzania.

This study aimed to create the basis for a long-term national forest monitoring system for Tanzania. Such a system is needed to help the country bridge the information gap and knowledge concerning remote sensing in forest monitoring at a national level as a cost-effective and timely means, instead of time-consuming and costly ground surveys. Specifically, this study will focus on the following research questions: (i) What is the current (baseline) spatial extent of forest cover in Tanzania? (ii) What is the distribution of the different forest types in Tanzania? (iii) How can changes in forest extent be mapped as part of an ongoing monitoring system?

2. Methods

2.1. Study Area

Our study area is located in mainland Tanzania. The country’s territory covers 945,100 km${}^2$ (29${}^{\circ }$ and 41${}^{\circ }$ E and 1${}^{\circ }$ and 12${}^{\circ }$ S) (Figure 1). Tanzania has diverse terrain, with a combination of plains, hills and forested mountains, with the highest peak at 5895 m a.s.l. The country has different bioclimatic and topographic zones, ranging from dry regions where precipitation levels are below 400 mm to humid areas where precipitation levels reach over 2000 mm per year [17], with a maximum mean temperature range of 26.6–33.1 ${}^{°}$C and a minimum of 5.3–18.3 ${}^{°}$C [18]. This broad diversity of geographic conditions has given rise to various ecosystems and habitats, including semi-arid, humid, and tropical subhumid zones [19]. The primary natural forest types are montane, lowland, mangrove forest and woodlands (open, closed and thickets) [20] and managed timberline trees, with different localised habitat distributions [21].

2.2. Software and Data Processing

Data acquisition and processing were undertaken using the open-source Remote Sensing and GIS Software Library (RSGISLib; [22]), the KEA file format [23] and the XGBoost [24] machine learning library. Earth Observation Data Downloader (EODataDown) was used for downloading and processing the Landsat-8 data to an analysis-ready data (ARD) product [25]. For the detection of change, EODataDown plugins were implemented, automating the scene-based processing following the generation of the ARD product. The EODataDown software is written in Python and uses a PostgreSQL database for meta-data storage with dependencies on RSGISLib and Atmospheric and Radiometric Correction of Satellite Imagery (ARCSI) [26]) software. For this study, the system was deployed on the SuperComputing Wales (SCW) high-performance computing (HPC) infrastructure. The use of SCW significantly reduced the processing time by paralleling the data processing, where approximately 100 cores were used at each processing stage.

2.3. Landsat-8 Pre-Processing

The forest baseline for Tanzania was developed using multispectral data from the Landsat 8 Operational Land Imager (OLI) [27]. At the time of analysis, all the collection-1 images (2013 to 2018) for May to November with a cloud cover threshold of >80% were downloaded. The images were selected from a dry period to minimise the cloud and cloud shadow and discrepancies in reflection caused by seasonal vegetation fluxes [28]. A total of 3200 Landsat 8 OLI images were downloaded using the Google Cloud API.

Pre-processing to surface reflectance was undertaken using the ARCSI [26] software, which uses a dark object subtraction (DOS) [29] in the visible bands to retrieve an estimate of the aerosol optical depth (AOD) [26] from the image. This is then used to parameterise the 6S radiative transfer model [30] and apply the resulting correction to create standardised reflectance using the topographic and bidirectional reflection correction proposed by Shepherd and Dymond [31]. More details of this processing chain are provided in [32].

2.4. Classification Methodology

The classification followed a hierarchical approach, first delineating a binary forest extent (level 1) and then classifying the forest pixels into forest types (level 2) (Figure 2), generating the forest baseline for Tanzania.

2.5. Forest/Non-Forest Classification

2.5.1. Defining Training Data

The training polygons for the forest/non-forest classification were generated with reference to the Landsat-8 imagery, higher-resolution Google Earth imagery and field knowledge. A total of 46,176 training polygons were collected (forest, n = 22,440 and non-forest, n = 23,736). These samples were then rasterised onto each of the 3200 images and the associated image pixels were extracted, creating 435,808,135 forest samples and 1,423,875,598 non-forest samples.

Training with this number of samples would create an XGBoost classifier with many levels using very large trees that would be very slow to train and apply. Therefore, 100 training sets were derived where 500 (of the 3200) scenes were randomly selected. The associated training samples were merged and then balanced to ensure an equal number of forest and non-forest samples in each training set. The number of samples within the 100 sets varied, with a minimum of 60,874,216 and a maximum of 78,833,485; the mean was 68,118,012 and the standard deviation was 35,37,652. Each of the 100 sets was then split by taking random non-overlapping samples to create training (50%), validation (25%) and testing (25%) datasets.

2.5.2. XGBoost Classifier

The XGBoost algorithm was used in this analysis. The XGBoost algorithm utilises a gradient boosting (GBM) approach that creates scalable tree structures designed for memory efficiency with parallel processing capability. The XGBoost algorithm also includes tree pruning and regularisation to avoid overfitting or bias to provide better classifier performance [24]. Therefore, the algorithm can use the large training datasets [33] produced as part of this study.

2.5.3. Optimising the XGBoost Parameters

The XGBoost algorithm has a large number of hyperparameters which can impact the classifier performance and therefore require optimisation. A Bayesian optimisation was therefore used to optimise the parameters for each of the 100 training sets. The full training and validation datasets were not used for the hyperparameter optimisation to reduce the time required to perform the optimisation. Randomly selected subsets of the training and validation data were therefore created; 60,000 training samples of forest and non-forest were used (i.e., 120,000 in total) and 20,000 for each (i.e., 40,000 in total) for the validation dataset.

2.5.4. Training the XGBoost Classifiers

Using the full training and validation datasets and the optimised hyperparameters, the 100 forest/non-forest classifiers were consecutively trained using 40 cores per job. It took 35 days for all 100 models to be trained. Using the testing datasets, the average accuracy of the classifiers was 99%.

2.5.5. Creating the Final Forest Extent Map

To generate the forest extent classification, the 100 classifiers were applied to each of the 3200 scenes. The 100 classifications were summarised on a per-pixel basis providing a percentage for the number of times a pixel was classified as forest. The scene was then thresholded using thresholds of 30%, 50% and 80%. The thresholds were selected based on a visual inspection of the subsample of scenes and chosen to capture the probability of pixels classified as a forest at the three different levels. To merge the scene-based classifications, to create a national forest mask, a 100 km tiling was used to allow parallel processing. The percentage of times it was classified as forest was calculated for each pixel, resulting in three output images for each scene-based threshold (i.e., 30%, 50% and 80%). Those outputs were subsequently thresholded using the same 30%, 50% and 80% thresholds, creating 9 forest extent maps for Tanzania (e.g., scene threshold of 50% and national threshold of 80%). An independent accuracy assessment was used to identify the optimal forest extent map.

2.5.6. Accuracy of the Forest Extent Map

The National Forest Inventory (NFI; NAFORMA) collected by Tanzania Forest Services from 2011 to 2015 [20] and other local forest inventories for the period 2016–2018 were considered for this assessment of accuracy. Still, the temporal- and spatial-scale differences in defining the forest extent from these data were difficult. Therefore, the NFI data were not considered reliable reference data to assess the forest extent map against.

Therefore, the accuracy assessment was conducted using stratified random sampling in 9 185 km by 180 km subregions distributed throughout Tanzania, ensuring that the validation took into account the variability in forest and non-forest land cover across the country. High-resolution images in base-map layers from virtual globe web-based maps (i.e., ESRI Satellite, Bing Satellite and Google Satellite), available through QGIS, were used as reference data as these layers have been found to be adequate for validation by previous studies (e.g., [34]) and enabled national coverage of the reference points rather than biased based on high-resolution imagery availability. For each subregion, 1000 sample points for forest and non-forest classes were generated, making 2000 points for each sample area and 18,000 reference points in total.

The ClassAccuracy QGIS plugin [32] was used to efficiently verify each point with an overlay on a virtual globes web-based map between the two classes. The accuracy metrics from the error matrix were summarised as both an overall accuracy (OA), user and producer accuracy (UA and PA), allocation disagreement (AD), quantity disagreement (QD) [35,36], F1-score [37], and the Matthews Correlation Coefficient (MCC) [38]. The metrics enable the users to understand the distribution of errors in the products.

2.6. Forest Type Classification

The second step (level 2) of forest classification was to classify the forested pixels into forest types (montane, lowland, mangrove, plantation forest, closed woodland, open woodland and thicket). The classification of forest types is necessary for generating detailed forest distribution to evaluate forest ecological systems and support monitoring and management practices. To constrain this analysis, each forest type habitat suitability, previously published in John et al. [21], was used. This novel approach was selected to minimise the classification error such that a pixel was only considered for the forest types the habitat suitability analysis had identified. Therefore, it constrained the classification of forest types based on their adaptation and corresponding bioclimatic patterns, minimising misclassification.

2.6.1. Forest Types Mask

The habitat suitability extent maps were intersected to merge the individual forest type suitability maps, identifying 34 combinations. For a small number of areas, the habitat suitability result provided suitability for only a single class (e.g., open woodland). However, this would not allow the classifier to perform a classification, so a second class was added in these cases. For example, for the areas which only had suitability for open woodland, then closed woodland was added (

Table 1. A sample of predicted combined forest types suitability depicting single classes occurrences. Therefore, a second class was added in these situations to enable the classifier to perform classification.

	Predicted Combination
Value	Habitat Type Predicted	Yes = 1/No = 0	Added Class
2	Montane forest	1	Lowland forest, plantation forest
	Lowland forest	0
	Mangrove forest	0
	Plantation forest	0
	Closed woodland	0
	Open woodland	0
	Thicket	0
4	Montane forest	0
	Lowland forest	0
	Mangrove forest	0
	Plantation forest	0
	Closed woodland	0
	Open woodland	0
	Thicket	1	Open woodland, closed woodland
8	Montane forest	0
	Lowland forest	0
	Mangrove forest	0
	Plantation forest	1	Montane, lowland forest
	Closed woodland	0
	Open woodland	0
	Thicket	0
16	Montane forest	0
	Lowland forest	0
	Mangrove forest	0
	Plantation forest	0
	Closed woodland	0
	Open woodland	1	Closed woodland, lowland forest
	Thicket	0
32	Montane forest	0
	Lowland forest	0
	Mangrove forest	0
	Plantation forest	0
	Closed woodland	1	Open woodland, lowland forest
	Open woodland	0
	Thicket	0
64	Montane forest	0
	Lowland forest	1	Montane forest, closed woodland
	Mangrove forest	0
	Plantation forest	0
	Closed woodland	0
	Open woodland	0
	Thicket	0
128	Montane forest	0
	Lowland forest	0
	Mangrove forest	1	Lowland forest, closed woodland
	Plantation forest	0
	Closed woodland	0
	Open woodland	0
	Thicket	0

).

The habitat suitability analysis was undertaken at a pixel resolution of 1 km as this was the resolution of the environmental variables used for the analysis, which allowed local environmental variability to be captured [39] at a country scale, and the corresponding bioclimatic patterns enable inferring of relationships between different forest type habitats. Therefore, a nearest neighbour resampling was used to create a 30 m resolution product required for the Landsat classification. However, due to this resolution change, there were 30 m pixels that were within the forest mask but did not have a habitat suitability. For example, along the coast and other forest/non-forest boundaries present at 1 km (e.g., wetlands and lakes). A k-Nearest Neighbour (k-NN) fill was performed to fill these pixels, where an unknown pixel was filled with the mode of the k spatially nearest pixels, $k=5$ was used for this analysis.

2.6.2. Defining the Training Dataset

Similar to the forest/non-forest classification, the training polygons were defined for each forest type and extracted from all 3200 Landsat scenes. A total of 20,370 sample polygons were defined for the forest types and resulted in 249,302,636 pixel samples (

Table 2. Summarise the training dataset.

Forest Type	Sample Polygons	Pixel Samples
Montane	1272	7,638,017
Lowland	2053	19,192,926
Mangrove	1407	2,225,183
Plantation	794	2,118,903
Closed Woodland	3264	34,636,347
Open Woodland	11,070	165,439,923
Thicket	510	18,051,337

).

To define the number of samples for each of the 34 class combinations (e.g., ‘mangroves’ and ‘lowland’), the data were combined, creating 34 training sets. The training data were balanced using random sampling to ensure an equal number of samples per class, avoiding bias towards the majority class. The maximum number of samples for a class was limited to 10,000,000 to reduce processing time. Therefore, for a combination of Mangroves, Lowland Forest, and Closed Woodland, the Lowland Forest and Closed Woodland samples were subset to 2,225,183 (i.e., the number of samples of mangroves as this was the smallest of the three classes). However, if the combination were Closed Woodland and Open Woodland, then the number of samples would have been limited to 10,000,000.

2.6.3. Training of the Classifiers

A single classifier was trained for each combination where the samples were split into training (50%), validation (25%) and testing (25%) sets. A 10% sample of the training and validation datasets were randomly extracted to optimise the XGBoost hyperparameters using Bayesian optimisation. The XGBoost classifiers were then trained using the optimised hyperparameters and full training and validation datasets. Using the independent testing dataset, the average classifier accuracy was 99%.

2.6.4. Final Forest Types Map

As with the forest/non-forest classification, the classification was applied on an individual scene basis. The habitat suitability derived mask was used to define the classifier applied to each pixel. Therefore, each pixel was only considered for the classes defined by the habitat suitability analysis. To summarise the scene-based forest type classifications, creating a national map, the mode of the scene-based analysis was taken on a per-pixel basis using the same 100 km tiles to allow for parallelisation of the processing.

2.6.5. Accuracy of the Forest Types Map

To assess the forest type map, the National Forest Inventory (NFI; NAFORMA; [20]) was used where it was masked using the forest mask from this study. Any remaining points within the NFI defined as being outside of a forested area were then removed and the points for the period 2011–2015 were checked for validity using virtual globe web-based maps in QGIS (e.g., Google Earth data). The temporal differences prevented the NFI data from being used to assess the forest extent map. However, to assess forest type classification, the NFI data were masked using the forest extent defined by this study, removing the majority of the temporal change from the dataset. If a pixel was defined as forest in our 2013–2018 baseline mapping and with the NFI data (2011–2015), it was unlikely that the forest type could have changed. Forest types are also more difficult to assess by analysing high-spatial-resolution remotely sensed data, particularly the differences between classes such as open and closed woodlands. Therefore, access to a field derived dataset is strongly preferred.

A final total of 13,200 NFI field points were used to assess the accuracy of forest type classification ($n=3895$ closed woodland, $n=1708$ lowland forest, $n=57$ mangrove forest, $n=401$ montane forest, $n=6721$ open woodland, $n=216$ plantation forest and $n=202$ thicket). The number of these points related to their spatial distribution across Tanzania with 6721 points for open woodland while just 57 for mangroves.

2.7. Forest Cover Change and Monitoring

2.7.1. System Architecture

One of the requirements for identifying change was the generation of a monitoring system, rather than just a classification of forest extent change. Therefore, the system was based on the EODataDown software system (Figure 3). EODataDown can be configured to automatically download and process Landsat, Sentinel-2 and Sentinel-1 data to an analysis-ready data product. The analysis is executed in date order with the oldest image first. The EODataDown can also execute a set of user-defined plugins that perform a set of data analysis tasks (e.g., the detection of change). Each time the system is executed, the latest imagery is downloaded and analysed. Using a tool such as cron, the system can be automated to run independently at a set time interval (e.g., daily or weekly) and create an automated monitoring system. An advantage of the EODataDown system is that it allows the end-user to focus on their data analysis, while the EODataDown architecture manages the data storage and processing; creating a monitoring system.

2.7.2. Landsat 8 Imagery

EODataDown uses the ARCSI software to generate optical ARD data products, and therefore the same ARD product is used for the change analysis as the baseline classification. However, to minimise the identification of False-Positives for change, a further processing step was applied where an addition ‘clear-sky’ mask was also derived through the RSGISLib software [22] and implemented as an EODataDown plugin.

The aim of the ‘clear-sky’ mask is to identify the large continuous areas of available imagery which have a ‘clear view to the sky’, discarding small isolated regions of imagery close to clouds that are sources of error within the cloud masking and ARD generation. The generation of the ‘clear-sky’ product is a two-step process. The first step is to buffer the identified cloud and cloud shadows by 30 km. Those regions outside of the resulting mask are clumped, and only those with a size greater than 3000 pixels are selected. Those regions are then grown, so they were not within 10 km (30 Landsat pixels) of a cloud or cloud shadow pixel.

While applying the ‘clear-sky’ mask reduces errors in the following change analysis associated with omissions in the cloud masking, it also reduces the extent of data available for the change analysis, including regions of valid data.

2.7.3. Forest Change Definition

Forest change is described as the complete or partial removal of forest cover (Figure 4) that causes changes in forest structure [40]. In the context of remote sensing, this needs to consider the resolution of the imagery being used, and at least three pixels ($30\times 30$ m), an area of approximately 0.27 ha, was considered the minimum mapping unit for this analysis. Defining the minimum mapping unit minimises the number of false-positive changes due to the complex land surface conditions, especially in savanna ecosystems, which are naturally more variable. The forest change detection is designed to detect and track abrupt forest change events from anthropogenic and natural catastrophes. It augments the forest information as a novel source of forest change map products in Tanzania. A summary of the forest change detection process is presented in Figure 4.

Within savanna ecosystems, fire is a natural process. Therefore, it is normal for these ecosystems to have a patchwork of burnt areas and rarely produces long-term changes to the ecosystem. The proper use of fire in woodland savannas is essential for maintaining these ecosystems. Early burning is carried out to reduce more severe fire damage later in the fire season. Therefore, in the context of this study, these changes are not considered “real change”. To separate fire from other changes within the analysis, the Normalised Burn Ratio (NBR) index (Equation (1); [41]) and Burn Area Index (BAI; Equation (2); [42]) were used. Thresholds of NBR $>-0.02$ and BAI$<100$ were used to define the unburnt areas, and it was within these areas that the remaining change analysis was undertaken. These thresholds were identified through a visual sensitivity analysis across a range of scenes and locations throughout Tanzania.

$$NBR=\frac{NIR-SWIR}{NIR+SWIR}$$

(1)

$$BAI=\frac{1}{{(0.1+R)}^2+{(0.6+NIR)}^2}$$

(2)

where in Landsat 8: NIR = band 5, SWIR = band 7 and R = band 4

Plantations are heavily managed forests within the landscape and are not considered in terms of national forest change statistics as they are already considered ‘changed’ and under anthropogenic modification. Change regularly occurs in plantation forests with partial forest loss after harvesting. However, replanting occurs shortly after and therefore, while the land cover may have temporally changed, the land use has not. To mask the plantation forest areas, a mask was extracted from the forest type classification.

2.7.4. Scene-Based Change Detection

The change analysis was undertaken in two steps. The first identified possible change pixels within the forest extent baseline previous defined, following masking for burnt areas and plantations. A normalised difference vegetation index (NDVI) threshold of <0.35 was applied to identify potential change areas. This threshold was identified based on expert knowledge of the environment and a visual sensitivity analysis across a number of scenes and time periods. To avoid seasonal changes (i.e., loss of leaves due to phenology), particularly within the savanna’s of western Tanzania, an additional threshold was applied to the whole scene where if >30% of the forested pixels were identified as a possible change, then the whole scene was ignored. At the extent of a Landsat scene, even following the removal of cloud, the extent of change within a scene will be small (i.e., probably much less than <1%). Finally, features of less than 3 pixels are removed from the layer to reduce noise.

Using the pixels identified in the first stage, the second identifies the change through classification, using the same 100 XGBoost classifiers trained to generate the forest extent baseline. Reducing the number of pixels classified significantly reduces the processing time for a scene, as only a small percentage of the total number of pixels within the scene is being classified. The 100 classifications are merged, and a threshold of 50% was applied to identify forest and non-forest regions. The resulting non-forest regions are considered the final change features for the scene.

2.7.5. Confirming Changes and Updating the Forest Baseline

On a per-scene basis, the identification of a change is considered to be of low confidence due to omissions in the cloud masking, topographic shadowing and misclassification of the forest extent. However, changes can be confirmed through multiple observations as errors are unlikely to coincide on consecutive scenes. Through experimentation, a threshold of 5 observed changes was selected to confirm a change. The changes were summarised to provide an annually updated forest baseline (Figure 3), as required for national reporting. The date of the first observed change was used to define the time point of the change occurring. EODataDown can be implemented as a scoring system where confirmed changes are provided as automated alerts to end users.

Change occurrence summaries were therefore generated for 2018, 2019 and 2020. The 2018 change layer included all changes identified against the forest extent for the 2013–2018 period of the imagery used to define the baseline and application of the 2018 changes to the baseline was considered to create a forest baseline with a discrete date of 1 January 2019. The change alerts for 2019 were identified and applied to generate a baseline for 1 January 2020. Finally, while the Landsat imagery was processed for 2020, few changes were confirmed, with 5 observations due to the frequency of cloud-free Landsat observations, meaning that many changes are only confirmed once the following years’ data (i.e., 2021) are available. Lowering the threshold to 3 observations increased the amount of change identified for 2020, indicating that change is similar to 2019. However, there was also an increase in the number of false positives. Therefore, the focus will be on reporting and validating the changes identified for 2019.

2.7.6. Forest Change Accuracy Assessment

Assessing change detection accuracy is difficult as change is rare, with only a small percentage of pixels changing within a given period (e.g., one year). Therefore, to assess the accuracy of the changes identified, relatively small but intensively sampled plots were used. A single plot was defined as a $1\times 1$ km area. The extent was observed as the most suitable trade-off, providing sufficient area for a representative assessment but being viable to interpret accuracy points intensively. 16 sample plots were selected in areas where changes were known to have occurred and stratified across the different forest types. For each $1\times 1$ km plot, 1000 reference points were randomly generated, with a minimum spacing of 30 m between the points. The minimum distance constraint was designed to ensure a pixel was only assessed once. In total, 16,000 verification points were used for the assessment. The verification points were not stratified using the generated change layer to avoid bias in the points’ location, but the densely allocated points should enable the estimation of change omissions, a significant challenge when assessing the accuracy of a change product [43].

The reference points were visually examined using the available Landsat 8 and, for some areas, Planet 3 m datasets. The ClassAccuracy accuracy assessment tool in QGIS [32] was used to efficiently assess the points. The images were acquired in October or November between the two target years (i.e., 2018 and 2019) to minimise seasonal variation were used for the analysis. For example, the validation was undertaken by observing the earlier image from the year 2018 and was confirmed using the later image of the year 2019, whether “real” or “false” forest change.

The evaluation of forest change was assessed using precision, recall, F1-score, user and producer accuracy, and overall accuracy metrics, which match binary classification and are widely applied in remote sensing classification methods [44]. Therefore, the evaluation focused on generating forest change information based on true positive, true negative, false positive, and false negative. The final forest loss accuracy assessment was compared with the available global forest change dataset from Hansen et al. [14] version 1.7 for 2019.

3. Results

3.1. Forest/Non-Forest Classification

The binary analysis (forest/non-forest) produced nine forest extent maps for Tanzania (

Table 3. Evaluating classification models performance using accuracy assessment metrics for binary classification (forest/non-forest).

Classification Model
Single-Scene (%)	Multi-Scene (%)	OA%	AD	QD	PC	TD	MCC
30	30	68.46 ± 0.50	0.009	0.202	0.787	0.212	0.50
30	50	81.77 ± 0.50	0.029	0.123	0.847	0.152	0.67
30	80	88.17 ± 0.40	0.065	0.047	0.887	0.112	0.75
50	30	73.16 ± 0.50	0.016	0.181	0.802	0.197	0.56
50	50	88.11 ± 0.40	0.049	0.061	0.889	0.110	0.77
50	80	85.93 ± 0.40	0.039	0.080	0.879	0.120	0.72
80	30	80.84 ± 0.50	0.034	0.127	0.838	0.161	0.66
80	50	89.66 ± 0.40	0.100	0.003	0.896	0.103	0.78
80	80	81.24 ± 0.40	0.014	0.111	0.873	0.126	0.64

OA: overall accuracy; AD: allocation disagreement; QD: quantity disagreement; PC: proportional correct; TD: total disagreement; MCC: Matthews Correlation Coefficient.

). To identify which map to take forward for further analysis, the accuracy of each map was considered.

3.1.1. Accuracy Assessment and Model Selection

The reported accuracy for the 9 models exhibited a satisfactorily level of OA ranging from 68.46 ± 0.50% to 89.66 ± 0.40% (Table 3). The best three models were further evaluated to select the final model. Therefore, the final chosen model (Figure 5) had a single-scene threshold of 80% and a multi-scene threshold of 50% with an overall accuracy of 89.66 ± 0.40%, F1-score of 0.87 and MCC value of 0.78 (

Table 4. Detailed accuracy metrics for the best three classification models.

Classification Model
Single-Scene (%)	Multi-Scene (%)	Cover Type	UA%	PA%	F1	P	R	C	O
80	50	Forest	87.69 ± 0.70	87.86 ± 0.60	0.87	0.87	0.87	0.053	0.050
		Non-forest	91.11 ± 0.50	90.98 ± 0.40	0.91	0.90	0.91	0.050	0.051
30	80	Forest	79.35 ± 0.90	91.59 ± 0.60	0.85	0.91	0.79	0.080	0.032
		Non-forest	94.65 ± 0.40	86.20 ± 0.50	0.90	0.86	0.94	0.032	0.080
50	50	Forest	95.02 ± 0.40	80.43 ± 0.60	0.87	0.80	0.95	0.024	0.085
		Non-forest	83.04 ± 0.70	95.79 ± 0.40	0.88	0.95	0.83	0.085	0.024

UA: user accuracy; PA: producer accuracy; P: precision; R: recall; C: commission; O: omission.

), sufficiently separated primary forests from non-forest classes. Classification errors were generally associated with wetlands and agricultural cropping being missed classified as forest. In contrast, the under-representation of forest extent was mainly associated with open woodlands where strong seasonal patterns, related to intermix with edaphic areas, especially in semi-arid regions and disturbances often by frequent fire and image availability may have contributed to the misclassification.

3.1.2. Forest Area Estimates

Table 5. Estimated forest area from the classification results for the nine models; compared with the National Forest Inventory (NFI)—NAFORMA [20].

Classification Model
Single-Scene (%)	Multi-Scenes (%)	Cover Type	Estimated Area (km${}^2$)	Area (%)	NFI Area (km${}^2$)	Area (%)
30	30	Forest	756,686	84.87
		Non-forest	134,913	15.13
30	50	Forest	612,041	68.65
		Non-forest	279,558	31.35
30	80	Forest	321,124	36.02
		Non-forest	570,476	63.98
50	30	Forest	708,875	79.51
		Non-forest	182,725	20.49
50	50	Forest	521,238	58.46
		Non-forest	370,361	41.54
50	80	Forest	249,171	27.95
		Non-forest	642,428	72.05
80	30	Forest	614,830	68.96
		Non-forest	276,770	31.04
80	50	Forest	407,976	45.76	481,000	54.4
		Non-forest	483,624	54.24	402,000	45.6
80	80	Forest	156,134	17.51
		Non-forest	735,466	82.49

presents the forest cover extent for Tanzania, as generated from the nine classification models. The results were compared with the previous national field inventory (2011–2015; [20]; NAFORMA). As shown in Table 3, Table 4 and Table 5, the classification identified as having the highest accuracy also resulted in forest area estimates closest to the NFI estimates [20]. Table 5 also demonstrates the relationship between the forest area mapped and the single-scene and multi-scene thresholds. As seen, the differences in the area mapped between the thresholds are also quite large (Table 5). This implies that an improvement in the classification accuracy might be possible with further refinement of the threshold selection and is an area for further study. Similarly,

Table 6. Summarised forest extent using a map ranked by regions.

Rank	Region	Area (ha)	Area (%)
1	Lindi	5,526,955	13.55
2	Ruvuma	5,059,867	12.40
3	Morogoro	4,295,711	10.53
4	Katavi	3,478,155	8.53
5	Tabora	3,054,801	7.49
6	Mbeya	2,163,746	5.30
7	Kigoma	1,988,670	4.87
8	Iringa	1,803,936	4.42
9	Pwani	1,739,369	4.26
10	Singida	1,689,458	4.14
11	Njombe	1,427,548	3.50
12	Songwe	1,220,258	2.99
13	Mtwara	1,203,790	2.95
14	Kagera	1,105,736	2.71
15	Tanga	951,165	2.33
16	Manyara	861,542	2.11
17	Rukwa	768,848	1.88
18	Geita	736,440	1.81
19	Dodoma	644,469	1.58
20	Kilimanjaro	377,028	0.92
21	Arusha	265,288	0.65
22	Shinyanga	177,709	0.44
23	Mara	92,345	0.23
24	Mwanza	89,022	0.22
25	Dar Es Salaam	43,970	0.11
26	Simiyu	31,773	0.08

summarises the forest extent by region in Tanzania.

3.2. Forest Types Classification

The novel forest type classification methodology, which used the forest habitat suitability to constrain the classification, resulted in an OA of 85%. This classification result provides important information on the current status of unique forest ecosystems and patterns in Tanzania. With a complex forest landscape with varying climatic conditions, from dry savanna to moist montane forest.

3.2.1. Accuracy Assessment

The overall accuracy of the forest type classification map (Figure 6) was 85%, with F1-scores ranging from 0.77 to 0.99. A quantity disagreement of 0.02 was calculated along with an allocation disagreement of 0.11, demonstrating that error in the classification has resulted in areas being misclassified but that these areas largely cancelled each other out to provide a more accurate overall area estimation (

Table 7. Thematic accuracy measures of the forest types classification.

Forest Type	UA(%)	PA(%)	F1	P	R	C	O
Montane forest	88.53 ± 3.10	89.65 ± 2.80	0.89	0.89	0.88	0.002	0.003
Lowland forest	96.02 ± 0.90	88.55± 1.30	0.92	0.88	0.96	0.005	0.015
Mangrove forest	98.24 ± 0.34	100 ± 0.00	0.99	1	0.98	0.000	0.000
Closed woodland	72.35 ± 1.40	82.76 ± 1.10	0.77	0.83	0.72	0.067	0.049
Open woodland	89.73 ± 0.70	85.23 ± 0.60	0.87	0.85	0.89	0.057	0.068
Plantation forest	77.78 ± 5.5	90.32 ± 4.00	0.84	0.90	0.77	0.004	0.001
Thicket	89.60 ± 4.20	79.04 ± 4.00	0.84	0.79	0.89	0.001	0.003
Overall accuracy (OA)	85.22 ± 0.50
Allocation disagreement (AD)	0.11
Quantity disagreement (QD)	0.02
Proportion correct (PC)	0.86
Total disagreement (TD)	0.13

UA: user accuracy; PA= producer accuracy; P: precision; R: recall; C: commission; O: omission.

).

The majority of the classification error (15%) was found between the deciduous forest types, closed and open woodland. These are challenging classes to define a boundary (Figure 7a,b) as the class definitions are described by a variation in tree cover rather than species composition. Therefore, the spectral difference between the classes is associated with the percentage of background soil and grassland reflectance verse canopy leaf reflectance. This is also associated with a broad geographic overlap between open and closed woodland and other forest communities (Figure 7c). For example, the possibility of habitat overlap between closed and open woodland is estimated at 80% and closed woodland and lowland is approximately 40% (Figure 7c). However, this study has demonstrated the ability to reliably retrieve quantitative information for forest types within the heterogeneous landscape of Tanzania and at a national level, meeting the requirements for national forest monitoring and reporting.

3.2.2. Forest Type Area Estimates

The forest type areal estimates were compared with the National Forest Inventory (NAFORMA), with some small differences between the forest types (

Table 8. Estimated area for forest type classification as compared with the National Forest Inventory (NFI) NAFORMA assessment [20].

Forest Types	Map Area (km${}^2$)	Area(%)	NFI Area (km${}^2$)	Area(%)
Montane forest	9716	2.35	9953	2.03
Lowland forest	60,670	14.65	16,565	3.38
Mangrove forest	767	0.19	1581	0.32
Closed woodland	93,004	22.45	87,290	17.79
Open woodland	237,052	57.22	359,973	73.37
Plantation forest	6695	1.62	5545	1.13
Thicket	6368	1.54	9719	1.90

) mostly associated with open woodland and related forest types such as closed woodland (Figure 7a,b). Open woodland is the most significant forest type by area, with 57%, followed by closed woodland with 22%. Therefore, woodlands occupy around 79% of the forest types, spreading from central to the western part of the country and with a mosaic of lowland forest along the coast and southernmost (Figure 6).

3.3. Estimated Forest Extent in Protected Areas

The forest cover of in situ conservation strategies such as protected areas is necessary for biodiversity and ecosystem protection in Tanzania. Therefore, this study also provides forest extent in forest reserves and wildlife managed areas (

Table 9. A summary of forest extent (ha) in protected areas. Area percentage computed from a total estimated forest area (40,797,600 ha) of the classification result (Table 5).

Protected Area Category
Cover	Forest Reserve	Area (%)	Wildlife Area	Area (%)
Forest	6,911,300	17	11,339,583	28

and

Table 10. Forest types in protected areas.

Forest Type (Area (ha))
Category	Mo	Lo	Ma	CW	OW	PF	Th
Forest Reserve	379,626	643,797	58,336	1,254,700	4,426,201	127,952	20,686
Wildlife Area	268,842	1,171,620	-	3,551,964	6,201,071	-	146,087

Mo: montane; Lo: lowland; Ma: mangrove; CW: closed woodland; OW: open woodland; PF: plantation forest; Th: thicket.

) as an ecological parameter required to produce desired conservation outcomes (Figure 8a,b).

3.4. Forest Cover Change Results

The accuracy assessment will focus on the 2019 change product as this is the only full year. The 2018 change contains the product of changes that occurred between 2013 and the end of 2018. While there were not sufficient observations to confirm the changes for 2020, the 2019 result demonstrated a countrywide wall-to-wall map of forest cover change over Tanzania (Figure 9) adequately detected forest changes from the baseline map and could be used for reporting annual forest change statistics. The forest loss area estimates for 2019 from this study were compared with the global forest change analysis of Hansen et al. [14] version 1.7, (

Table 11. Forest loss area for the year 2019 (this study) compared with global forest cover loss from Hansen et al. [14] version 1.7.

	Forest Cover Loss Area (ha) and Area (%)
Class	This Study	%	Hansen et al. [14]	%
Forest Change	157,204	0.39	142,773	0.36

) and found to be comparable.

3.4.1. Accuracy Assessment

Table 12. Model performance evaluation metrics for potential forest changes in 2019 at a 95% confidence interval, compared to global forest change analysis version 1.7 of 2019 Hansen et al. [14].

Class	Measure	This Study	Hansen et al. [14]
Change	Producer accuracy(%)	96.13 ± 0.74	65.98 ± 2.04
No-change		92.74 ± 0.36	84.21 ± 0.33
Change	User accuracy(%)	71.42 ± 1.52	34.25± 1.59
No-change		99.22 ± 0.15	95.20 ± 0.37
Change	Precision	0.96	0.65
No-change		0.92	0.84
Change	Recall	0.71	0.34
No-change		0.99	0.95
Change	F1-score	0.82	0.45
No-change		0.96	0.89
	Overall accuracy(%)	93.28 ± 0.38	82.19 ± 0.59

summarises the accuracy assessment results for the changes identified in 2019. The accuracy was found to be good with an F1-score of 0.82 compared to the global forest change assessment of Hansen et al. [14] with an F1-score of 0.45, highlighting that the Hansen et al. [14] product is not capturing the full extent of change with Tanzania. For the no-change class, the accuracy was similar, with an F1 score of 0.96 for this study and 0.89 for Hansen et al. [14]. It would be expected that the no-change results will be similar, with change only representing a few per cent of the landscape. Therefore, most accuracy assessment points will be in no-change regions. Even significant errors within the change result would only result in small changes to the no-change class extent. This difference in the accuracy of the change product also demonstrates an improvement in forest change data quality. This result also highlights the importance of locally optimised analysis methods compared to being reliant on global datasets for national reporting.

3.4.2. Estimated Forest and Forest Type Change by Region

The forest change results were also summarised at a localised level to indicate regions with high deforestation rates (

Table 13. Summarise forest change extent by region in Tanzania.

Area (ha), Year
Region	2018	2019	2020
Tabora	211,412	31,587	5332
Katavi	120,629	23,171	2030
Rukwa	74,957	13,378	1360
Mtwara	30,431	12,633	2227
Mbeya	84,983	10,582	4185
Lindi	35,004	10,223	780
Singida	32,096	10,175	4268
Kigoma	76,750	8926	893
Songwe	72,716	7020	2945
Iringa	22,470	6033	664
Ruvuma	30,246	5991	638
Morogoro	16,674	5745	92
Geita	33,318	3630	338
Pwani	8120	1769	57
Shinyanga	10,449	1646	125
Njombe	5529	1384	111
Dodoma	5233	759	1054
Kagera	5823	411	57
Mwanza	3455	374	38
Kilimanjaro	27	28	-
Tanga	689	18	10
Mara	1453	18	10
Simiyu	91	10	6
Manyara	456	9	3
Dar Es Salaam	10	2	-
Arusha	19	1	-

). The majority of forest changes identified are within the open woodlands, primarily due to shifting cultivation. However, closed woodlands and lowland forests have also witnessed significant change (

Table 14. Forest type change by regions.

	Forest Type Change (Area (ha))-2019
Region	Mo	Lo	Ma	CW	OW	PF	Th
Tabora	-	-	-	212	30,815	-	416
Katavi	-	2	-	687	21,256	-	-
Rukwa	8	2	-	1765	11,225	-	-
Mbeya	8	10	-	267	9860	-	-
Singida	-	-	-	40	7677	-	2441
Kigoma	11	20	-	1003	7432	-	-
Songwe	11	13	-	113	6689	-	-
Iringa	5	11		131	5382	-	-
Lindi	-	3599	3	1221	4921	-	-
Mtwara	-	5891	-	869	4893	-	-
Ruvuma	7	154	-	843	4766	-	-
Geita	-	1	-	106	3519	-	-
Morogoro	13	1116	-	2582	1669	-	-
Shinyanga	-	-	-	13	1633	-	-
Njombe	9	151	-	53	709	-	-
Dodoma	-	-	-	83	630	-	16
Pwani	-	688	-	562	421	-	-
Kagera	-	1	-	18	390	-	-
Mwanza	-	81	-	48	239	-	-
Kilimanjaro	-	-	-	2	25	-	-
Simiyu	-	-	-	-	10	-	-
Tanga	3	7	-	5	8	-	-
Manyara	-	-	-	-	5	-	3
Mara	-	-	-	10	4	-	-
Dar Es Salaam	-	-	-	-	1	-	-
Arusha	-	-	-	-	-	-	-

Mo: montane; Lo: lowland; Ma: mangrove; CW: closed woodland; OW: open woodland; PF: plantation forest; Th: thicket.

, although it should also be noted that the extent of change within each forest type corresponds with the area of that type (i.e., open woodland, closed woodland and lowland forest are the three most extensive forest types by area in Tanzania) (Table 8).

3.4.3. Estimated Forest Change in Protected Areas

Although forest loss is more pronounced outside the protected areas, the result also highlights that forest loss is also occurring in protected areas (forest reserves and wildlife areas) (

Table 15. Estimated forest cover change in protected areas.

Area (ha), Year
Category	2018	2019	2020
Forest reserve	126,925	24,869	3600
Wildlife protected area	72,153	13,587	2842

and Figure 10). These protected areas are designed to protect and support the country’s biodiversity, and therefore changes within these areas are particularly significant. An example is shown in Figure 10, where deforestation in the protected area is occurring due to encroachment from the boundaries, with changes occurring up to 500 m from the protected area boundary. This is particularly concerning, as unchecked, these protected areas could witness further encroachment increasing the vulnerability of these important habitats and ecosystems.

3.4.4. Updating Earlier Forest Baseline

The confirmed forest changes for 2018 and 2019 were used to update the forest extent baseline to generate national forest baselines for 2018 and 2019.

Table 16. Forest loss extent comparison with Hansen et al. [14] version 1.7 for the period 2013–2019.

Forest Cover Loss Area (km${}^2$)
Class	This Study	Hansen et al. [14]
Forest change	10,462.37	11,539.11

present the comparison of forest loss from this study as compared to the global dataset of [14] version 1.7 over six years (2013–2019). It should be noted that these studies are not directly comparable for the period 2013–2019 as while the [14] product is producing an annual change product for each year, this study is using a baseline forest mask which is the product of imagery from 2013 to 2018 and not a 2013 map of forest extent. Therefore, while the accuracy assessment of 2019 suggests that Hansen et al. [14] is underestimating the true extent of change in Tanzania, this study has a lower value as it is against a different baseline. Only the 2019 annual change is directly comparable to the Hansen et al. [14] products.

Table 17. Update of forest baseline extent from the detected forest loss for the period 2013–2019.

Forest Area (km${}^2$) for Earlier Baseline, Updated Baseline and Area (%)
Class	Earlier Baseline Extent	%	Estimate Forest Loss	%	Updated Baseline Extent	%
Forest	407,976	45.76	10,462.37	2.56	397,514	43.20

summarises the updated baseline forest extent from the forest loss detected from the earlier baseline.

4. Discussion

4.1. Summary of Results

Our study demonstrated the potential application of open-free software, freely available satellite data, and advanced remote sensing techniques to provide a cost-effective method to obtain wall-to-wall information on the forest extent (Figure 5a,b and Figure 6) and associated changes (Figure 9) in Tanzania. It bridges the information gap and knowledge concerning the use of remote sensing for generating forest information about status, extent, types and location over a larger geographical area (national level) [46]. Therefore, the forest baseline (reference level) represents the extent (intact) forest area required to monitor future forest cover in Tanzania. The information is vital for developing practical, long-term plans to conserve and manage biodiversity, based on forest extent, type and composition [46].

The forest/non-forest classification model achieved an overall accuracy of 89.22% with an F1-score of 0.87 (Table 4). However, this still resulted in a 10% error in the classification. From a visual assessment of the map (Figure 5a), the likely sources of error were ascertained to be forest intermixed with edaphic areas, especially in woodland areas and the forest-grassland mosaic that remains evergreen throughout the year. Yet, the results are sufficient for the reporting of forest extent in Tanzania. This study introduced a novel method using habitat suitability modelling [21] to constrain forest type classification (Figure 5b) such that only appropriate forest types were considered for classification inappropriate geographical regions. For forest type classification, an overall accuracy of 85% with an F1-score ranging from 0.77 to 0.99 was achieved (Table 7). A particular challenge for the classification of forest types (15% error) was differentiating closed and open woodland areas (Figure 7a–c), as the boundary between these classes is based on the tree cover rather than the species composition of the woodlands. Future studies could also consider approaches that aim to retrieve associated biophysical parameters such as canopy cover. However, the result is considered the best mapping of Tanzania currently available and could be applied to other neighbouring countries (e.g., Kenya and Mozambique) which have similar ecosystems. These maps will help to establish a structure and long-term forest monitoring system in Tanzania. Forest cover information is needed to support the national forest policy to sustainably manage, conserve, restore and utilise forests and associated resources for Tanzania’s socio-economic growth and climate resilience.

A further consideration is that this study used imagery over 5 years (2013–2018) to mitigate the issues of data availability given the high level of cloud cover, particularly in the coastal areas. However, during this period, change will have occurred in the forest extent. Therefore, the forest extent and type maps represent the forest cover for the majority of the scenes within the period. Thus, the maps were updated to provide maps for 2018 and 2019 through the change detection result (Section 3.4.4).

4.2. Forests and Forest Types Extent

The present analysis estimated a forest cover area of 407,976 km${}^2$, representing 45.76% of the country landmass (Table 5). The forest type classification (Figure 5b) result indicated a prominent class of woodlands (closed, open woodland and thickets), estimated to cover approximately 336,405 km${}^2$, which make up 81.20% of the forested land (Table 8), an important ecosystem of great significance to human economies [1,2,3], mainly covering the central and western part of the country. The montane forests represent biodiversity hotspots along a chain of isolated mountain ranges (Figure 5b), supporting a diversity of endemic species [47], an area of approximately 9717 km${}^2$ representing 2.35% of the forest cover (Table 8). Similarly, montane forests harbour the world famous tropical montane rivers, including the Eastern Arc Mountains [48], feeding major rivers, floodplains and ocean. The proposed construction of the large Julius Nyerere hydro-power station across the Rufiji river will depend on forest conservation, especially upland montane forests, to reduce siltation, but also the Selous Game Reserve and Nyerere National Park, home to a wealth of flora and fauna as the long-term resource sustainability base for the nation at large. The lowland forest habitat overlapping with montane forest and woodlands with the most significant biological value and source of water supply for wildlife and people was estimated to cover 60,718 km${}^2$, representing 14.16% of the forest cover next to closed woodland (Table 8). Therefore, the result will support developing diverse conservation states for different forest types that minimise overexploitation, especially on fragile sites.

In the previous National Forest Inventory (NAFORMA) released in 2015, the area was estimated at 481,000 km${}^2$ representing 54% of the total land area [20]. Arguably, the main observational gap between these analyses should not be considered due to forest loss in Tanzania, but due to differences in timing, methodology and accuracy in the two studies, i.e., wall-to-wall mapping using Earth Observation data viz., sample-based forest inventory plots with a relative sampling error of 8.89% [20]. There are challenges reaching some areas with traditional forest inventory and pushing for sampling intensity reduction and focused sampling efforts, with few samples/plots selected from these areas [49], particularly in the mountains. The forest is evenly distributed over the country. The top three regions with extensive woodland areas include Lindi, Ruvuma and Morogoro, followed by the western part of the country, covering extensive dry miombo woodlands stretching across Tabora, Katavi, Kigoma and parts of the Rukwa region (Figure 5b). Overall, these areas occupy essential protected areas (Figure 8) in the country and are the cornerstone of forest and biodiversity conservation in Tanzania [50]. Therefore, increasing and maintaining a well-connected system of protected areas is a viable conservation strategy as a natural solution to global challenges, including climate change and deforestation [1].

Similarly, the increase in industrial forest plantations from government, communities and individual farmers supports the increase in forest extent in the southern highlands (Figure 6) [51]. For example, approximately 564,678 ha of plantation forest are found in the three regions of Iringa, Njombe and Mbeya. Accordingly, the result provides a consistent forest extent at a national level, whereby conservation policy actions can be planned and evaluate future forest changes and carbon storage assessment, e.g., Suarez et al. [52].

4.3. Scene-Based Forest Change Detection

Seasonality changes and persistent cloud cover in Tanzania create low data availability and excessive gaps (missing data) in the Landsat archive. Yet, the scene-based change detection method proposed in this study was found to overcome these challenges and achieved an overall accuracy of 82% (Table 12). However, future work could also include other pixel-based time series change detection methods such as Breaks For Additive Seasonal and Trend (BFAST; [53]), BFAST Monitor [54], Continuous Change Detection and Classification (CCDC; [55]), Jumps Upon Spectrum and Trend (JUST; [56]) and Exponentially Weighted Moving Average Change Detection (EWMACD; [57]). These methods have been demonstrated to be applicable to a large range of land cover change problems (e.g., [58,59,60,61]), primarily with Landsat but also Sentinel-2. For Tanzania, the availability of a longer time series of Sentinel-2 imagery will likely make such approaches viable at a national scale. The high levels of cloud cover, particularly in the East of the country, can make optical remote sensing difficult.

4.4. Forest Change Area Estimates

The forest cover change estimates provide essential information to guide policy formulation and implementation in protecting forests with better decision making in government programmes and other forest protection fiscal incentive projects [62]—notably, for addressing issues such as the UN Reduced Emissions from Deforestation and Forest Degradation (REDD+) and Green House Gas (GHG) as international reporting commitments and the 2030 Agenda for Sustainable Development Goals (SDGs) over fighting deforestation [16]. The result identified 157,204 ha of forest loss for 2019 in Tanzania, presenting 0.39% loss of intact forest, close to the global forest change analysis of Hansen et al. [14] (version 1.7), which mapped 142,773 ha in the same period. However, this study’s forest change analysis presented a methodology that was optimised for Tanzania and therefore had a higher degree of accuracy (Figure 9) and given that change is infrequent, even large changes in the accuracy of the change detection algorithm can result in relatively small changes in the area estimated. However, these relatively small geographic areas can be significant if found to occur in areas of importance (e.g., protected areas). Therefore, it can be considered that global forest change datasets remain suitable for providing an indicative trend of forest loss at a national scale [63,64] but locally optimised products are preferred for national and regional management decision making. Therefore, the forest change area estimates from this study provide an essential reference point in the region to which the Hansen et al. [14] product can be compared. Few countries have established a wall-to-wall forest change map or method for long-term national-scale forest monitoring.

The prime three regions with high deforestation include Tabora, Katavi and Rukwa (Table 13). The forest cover change analysis was achieved with an accuracy of 82% compared to 45% from the global forest change analysis [14] (Table 12). The changes detected by the monitoring system were also used to update the forest baseline map. Hence, the baseline map (forest mask) was updated from 407,976 to 397,514 km${}^2$ by 2019, indicating a decrease of 2.56% of the forest cover in Tanzania over six years (2013–2019) (Table 17).

The developed deforestation monitoring methodology aimed to provide the ability to respond immediately to reduce or stop the newly detected illegal deforestation situation from further expanding. In the future, the proposed change system might be considered for the formation of an early warning system. For early warning, the number of scenes used to confirm a change can be reduced (e.g., two observations) but would require visually checking to confirm the change and is recommended for forest guardians working in areas of high importance, such as protected areas under pressure from deforestation rather than a national system. For the annual reporting of forest loss for policy makers, using five observations of change to confirm a true change has been demonstrated to produce reliable results which are fit for this purpose.

4.5. Forest Management Outlook

Deforestation is detected beyond general-use land in protected areas (Table 15) and (Figure 10) typical of explicit forest conservation status, indicating a threat to the last remnant of important tropical forests in Tanzania. Protected areas and other conservation strategies support forest protection. However, they will become increasingly isolated and fragmented (Figure 11) as surrounding forest land is removed, turning to non-forest. Similarly, these protected areas will lose all of their protected status, leading to conversion to other land uses. Likewise, forest loss will increase the isolation (patches) of the protected areas impacting wildlife corridors. For example, Figure 10a at the Msaginia forest reserve supports wildlife movement between Katavi National Park and other protected areas, but the ongoing forest loss will limit this corridor.

Similarly, the detected deforestation in the western part of the country (Figure 9i), if it remains unchecked, will disrupt water flow, increase soil erosion, threaten the Malagarasi river that supports the Malagarasi-Muyovozi Ramsar site, and increase the siltation of Lake Tanganyika. It will raise the severity of the flooding that is already occurring and disrupting livelihoods [65]. Therefore, it can be argued that forest cover loss is endangering Tanzania’s economy. Increasingly, natural catastrophes such as droughts and El Niño climatic crises have influenced much agricultural productivity, power generation and transportation. Consequently, the forest monitoring system aimed to enhance forest law enforcement in protecting forests with better decision making.

5. Conclusions

Timely tropical forest monitoring is required to provide information about forest extent and changes over time, reducing data gaps necessary for forest conservation, management and responding to climate change for sustainable development. This study has provided the first consistent and robust forest cover extent with an area of 407,976 km${}^2$ (45.76% of the country’s area) and an overall accuracy of 89%. The estimated forest loss for 2019 was 157,204 ha, with an overall accuracy of 82%, contributing essential information for both science and forest management in Tanzania. These results have improved the quality of information available for Tanzania, which were previously considered inadequately in terms of quality and coverage to provide a baseline for national reporting. The forest monitoring system developed through this study is intended to link policy making on forest conservation and protection to meet national forest data requirements and integrate them into national institutions. It will enhance conservation programmes, which are rearguard efforts to save the last remnants of pristine forests remaining in Tanzania.

Innovative methods were developed through this study, including constraining forest type classification using the results of a habitat suitability analysis that achieved an overall accuracy of 85%. The methods used could be directly applied to other optical remotely sensed data, particularly Sentinel-2, which for the monitoring of change would provide an increase in the temporal frequency of the observations and therefore speed up the identification and confirmation of changes. In the future, SAR data such as Sentinel-1, ALOS-2 PALSAR-2 and NISAR data could also be integrated into the monitoring system, reducing the impacts of clouds, haze, and dust, common in Tanzania. Future studies might also consider removing the single-scene and multi-scene thresholds and merging all individual forest/non-forest classifications to create a single probability surface for forest extent. A local sensitivity analysis could be carried out to derive the optimal thresholds within biogeographic regions rather than nationally.