Analysis of land cover dynamics in Mozambique (2001–2016)

Land cover change (LCC) is a complex and dynamic process influenced by social, economic, and biophysical factors that can cause significant impacts on ecological processes and biodiversity conservation. The assessment of LCC is particularly relevant in a country like Mozambique where livelihood strongly depends on natural resources. In this study, LCC was assessed using a point-based sampling approach through Open Foris Collect Earth (CE), a free and open-source software for land assessment developed by the Food and Agriculture Organization of the United Nations. This study aimed to conduct an LCC assessment using CE for the entire Mozambique, and according to three different land classifications: administrative boundaries (provinces), ecoregions, and protected vs unprotected areas. A set of 23,938 randomly selected plots, with an area of 0.5 hectares, placed on a 4 × 4 km regular grid over the entire country, was assessed using CE. The analysis showed that Mozambique has gone through significant loss of forest (− 1.3 Mha) mainly to the conversion to cropland. Deforestation is not occurring evenly throughout the country with some provinces, such as Nampula and Zambezia, characterized by higher rates than others, such as Gaza and Niassa. This result can be explained considering a combination of ecological and socio-economic factors, as well as the conservative role played by the protected areas. Our study confirmed that LCC is a complex phenomenon, and the augmented visual interpretation methodology can effectively complement and integrate the LCC analyses conducted using the traditional wall-to-wall mapping to support national land assessment and forest inventories and provide training data for environmental modeling.


Introduction
Land cover change (LCC) is a complex and dynamic process influenced by social, economic, and biophysical factors that can cause significant impacts on land productivity, biodiversity conservation, and greenhouse gases emissions and removal (Geist and Lambin 2002;Kleemann et al. 2017;Woollen et al. 2016). In Mozambique, LCC has been considered a tool to monitor the achievements of several international commitments including the Aichi Targets of the United Nations Convention on Biological Diversity (UNCBD), the REDD + (Reducing Emissions from Deforestation and forest Degradation) mechanism of the United Nations Framework Convention on Climate Change (UNFCCC), the Land degradation neutrality (LDN) initiative of the United Nations Convention to Combat Desertification (UNCCD), and the Sustainable Development Goals (SDG), particularly the SDG target 15.3 dedicated to the restoration of degraded land and soil. Moreover, LCC has been used to support the assessment of ecosystem services and their spatial and temporal patterns (Atumane and Cabral 2021;Mahamane et al. 2017;Niquisse et al. 2017). This is particularly relevant in a country like Mozambique where livelihood strongly depends on natural resources, including wood harvesting for construction, firewood, charcoal, collection of edible and medicinal plants, livestock grazing, and subsistence farming (Nhantumbo et al. 2001;Woollen et al. 2016;Nicosia et al. 2022). Such human pressure on natural resources will be exacerbated by population growth and the increased demand for agricultural land, especially large commercial farms, and an increase in wood harvesting (Ryan et al. 2014;Filipe and Norfolk 2017;Instituto Nacional de Estatística 2020). To monitor such dynamics, the LCC assessment can be conducted using a wall-towall or a point-based sampling approach. The former can benefit from readily available global products such as the Climate Change Initiative of the European Space Agency (ESA CCI) land cover map (European Space Agency 2017) or from the development of maps supported by tools such as Google Earth Engine cloud computing capacities (Bey et al. 2020;Bey and Meyfroidt 2021;Lopes et al. 2019;Mahamane et al. 2017;Mananze et al. 2020;Tokura et al. 2020). The latter can be based on the use of Open Foris Collect Earth (CE), a free and open-source software for land assessment developed by the Food and Agriculture Organization of the United Nations (Bey et al. 2016), which has been already used to assess at global level dryland forests (Bastin et al. 2017), at national level the forest degradation processes in Papua New Guinea (Gamoga et al. 2021), and at a local level the spatial pattern of tree cover and density in Southern African national parks (Messina et al. 2018). This study aimed to conduct a LCC assessment using CE for the entire Mozambique, analyzing and comparing the results according to three different land subdivision systems: administrative boundaries, ecoregions, and protected vs. unprotected areas. Data and results provided by this analysis can be used to inform and support conservation and land use planning initiations at a national scale.

Study area
Mozambique is located in southeastern Africa at latitude 10°27′S to 26°52′S and longitude of 40°51′E to 30°12′E (Fig. 1) and has a total surface area of about 800,000 km 2 (Instituto Nacional de Estatística 2020). Approximately 70% of its surface is covered by forests and woodlands, and 29% of the Country's land surface and 2% of the marine areas are currently designated as protected areas (UNEP-WCMC 2022). The Zambezi River crosses Mozambique from east to west flowing toward the Indian Ocean and cuts the country roughly into two main regions: the southern region, dominated by lowlands with an altitude ranging from 0 to 200 m, and the northern region, which consists of a large plateau ranging between 200 and 600 m and reaching up to 1000 m. In the highland region, Mount Biga, the highest peak in the country, reaches 2436 m (Ministry for the Coordination of Environmental Affairs -MICOA 2014). The country's soils can also be subdivided into two main classes: sedimentary soils, characterizing the southern and the east/ central regions, and the ancient granite rock basement of Africa which underlies most of the northern and west-central regions. Most of the country falls in the tropical climate region, with a sub-tropical climate in the south, characterized by two main seasons: a cooler dry season, from April to October, and a warmer humid season from October to April (Barbosa et al. 2001). The average annual temperature is about 24 °C, with the warmest temperatures measured near the coast and a colder climate further inland. Rainfall distributions fluctuate widely through the country, following a north-south gradient. Annual average rainfall varies between 800 and 1200 mm along the coast, and from 600 to 1000 mm in inland areas (World Bank 2017).

Data collection
The LCC assessment was performed using Collect Earth (CE) version 1.5 (Open Foris n.d.), which facilitates the access to multiple freely available archives of satellite imagery, including those with very high spatial resolution imagery (Google Earth and Bing Maps) and those with very high temporal resolution imagery (Google Earth Engine Explorer and Code Editor). Collectively, these archives offer free access to information on current and past land dynamics for any location in the world. Therefore, CE draws upon these archives and the synergy of imagery at multiple resolutions to enable an innovative method for land monitoring, referred to as augmented visual interpretation (Bey et al. 2016). The assessment units were constituted by square plots, each with an area of 0.5 hectares (approx. 70 × 70 m). The plots were placed on a 4 × 4 km regular grid, covering the country, for a total number of 47,876 plots. The plot size of 0.5 ha was chosen to be consistent with the FAO definition of forest, which has a tree cover ≥ 10% spanning an area of more than 0.5 ha that is not predominantly used for agriculture or urban activity, as well as areas in which tree cover is temporarily < 10% but is expected to recover (FAO-FRA 2020). This grid has been defined in collaboration with the National Directorate of Forestry and is the same used for the National Forest Inventory (NFI) of Mozambique. A subset of 23,938 randomly selected plots was used for the assessment. Satellite imagery available for the interpretation had a spatial resolution of ≤ 1 m (VHR) for ∼96% of plots, ≤ 10 m (high resolution, e.g., SPOT, RapidEye, and Sentinel 2) for ∼4% of plots, and > 10-100 m (medium resolution, e.g., Landsat) for less than 1% of plots, and was complemented by the analysis of the Normalized Difference Vegetation Index (NDVI) computed using Landsat 8 data (Bey et al. 2016).
Through the timeline tool of Google Earth, all the available high-resolution images of the plot area can be assessed in the same software. The acquisition dates of the high-resolution satellite images available for each plot (considering both Google Earth and Bing Maps depositories) spanned from 2001 to 2016, although most of them (> 87%) were acquired in 2013 or later. Due to the absence of high-resolution satellite images before the year 2000, this year has been chosen as the lower temporal limit of the assessment.
The plots were randomly assigned to 50 different grids (504 plots per grid) and subsequently randomly assigned to different operators. Each grid contained plots placed in all provinces, to compensate for possible bias due to differences in the assessment approach by the single operators. The assessment was carried out in three separate time periods, each followed by a quality control assessment of the collected data. To help the visual evaluation, each plot was subdivided into 49 (7 × 7) subplots, each representing ~ 2% of the plot area. This subdivision is used to assess the percentage cover of different elements (trees, shrubs, grasses, buildings, huts, roads, paths, agricultural fields, internal water bodies, sand, rock, and bare ground) when the plot is placed on a land area (Fig. 2). The land cover elements in the sampled plots were assessed through 14 classes of minimum percentage cover for each different element (0%, 2%, 4%, 6%, 8%, 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, and 90%).
Recognition of different elements is not always simple and straightforward, and to distinguish one element from another, knowledge of local land use practices is essential. While infrastructure is generally well recognizable, the distinction between natural elements such as trees and shrubs or grassland and cultivated fields can be tricky.
As a general rule, trees are identified by a recognizable shadow, while shrubs do not project a recognizable shadow and are generally under 3 m tall. Plants with height > 3 m are considered trees and can be identified by comparison with surrounding elements. Crop fields are generally recognized based on their regular shape and traces of tillage or the presence of symmetrically distributed elements.
To calculate the occupied area by each land cover class, the following formula was used: where A i is the area (ha) for each i land use class, N i is the total number of plots assessed for each land use class, and E f is the expansion factor, resulting from the ratio of the total surveyed area A (Mozambique: 801,537 km 2 ) and the total plots assessed N. In a sampling survey, the expansion factor is used to expand the response for a single unit to an estimate for the whole set of units. In this case, the single unit is represented by the single plot (0.5 ha) and the whole set of units by the total area of the country (A). The expansion factor therefore indicates the amount of area that each plot represents in the survey, in this case being 3348.4 ha (A/N).
Once the plots have been assessed and the number N i of plots assigned to each land cover class i has been obtained, an estimate of the total area of the country covered by each land cover class i is obtained by multiplying the expansion factor by the number of plots N i .
Of the 23,938 plots initially selected, 148 have been discarded due to the unavailability of VHR images, which made it impossible to assess land cover. This resulted in a total of 23,790 plots, distributed in the ten provinces of the country as reported in Table 1.
To define the previous land uses and monitor the changes occurred in the past, the timeline tool of Google Earth allows to navigate through the archive of past images for the same plot area, while the grid interface provided by Collect Earth for each plot it is maintained also for the past images, allowing the repetition of the assessment procedures for any available images. The operator has also access to a series of accessory data obtained from Google Earth Engine, as NDVI charts and MODIS fire observation. The total corpus of data allows the operator to monitor even land use changes not directly visible through the comparison of present and past high-resolution images.

Classification systems
The LCC was conducted according to the following land cover classification systems: (1) IPCC (IPCC 2006) (2) FAO Global Forest Resources Assessment (FAO-FRA 2020).
However, in this study, only the results obtained according to the IPCC classification are reported.
Six main land cover classes were used in a hierarchical order in the assessment, following the IPCC categories (IPCC 2006): Settlement > Cropland > Forest > Grassland/ Savannah > Waterland > Other land. The classification of each plot has been obtained analyzing the distribution and functional relationship between land elements (trees, shrubs, grasses, buildings, huts, roads, paths, agricultural fields, internal water bodies, sand, rock, and bare ground) using the 49 subplots grid as a reference system. Each subplot has been classified in one of the main classes depending on which element it covered (ex. Subplot over tree-> Forest; Subplot over meadow-> Grassland; Subplot over road-> Settlement) (Martínez and Mollicone 2012). Each plot was then categorized as belonging to the first class to reach a 20% subplots' threshold (IPCC 2006) following the hierarchical order mentioned above: starting from Settlement and following the hierarchical order, the plot was classified in the highest class to reach ten subplots (20% of the total 49).
As mentioned above, available high-resolution satellite images from Google Earth and Bing Map to assess each plot spanned from 2001 to 2016. Therefore, following the described methodology for baseline year (oldest available image) and target year (most recent available image), LCC dynamics were assessed within this time framework.
To identify and discuss specific patterns in LCC, the data collected have been analyzed and compared against three land subdivision systems: Administrative division-Mozambique is divided into 11 administrative areas called provinces, including the special region of Maputo City (Fig. 1), which retain a certain amount of autonomy in the planning and management of the land (Informação Geral-Portal Do Governo de Moçambique).
Ecoregions-To determine if there has been a selective pressure causing different LCC patterns on distinct ecosystems and biological communities in the country, an evaluation was conducted over a simplified version of the ecoregions' map defined in Dinerstein et al. 2017. The original 13 classes have been reduced to 5 (Fig. 3)  Conservation status-As much as 29% of Mozambique's land is covered with different types of protected areas (Fig. 4), ranging from the National Parks to Game Reserves (Conservation Areas-ANAC, n.d.). We selected protected areas whose focus is connected to environmental or biodiversity conservation. These cover a total surface of 111,442 km 2 and include Special Reserves, Sanctuaries, National Reserves, National Parks, Forest Reserves, Environmental Protection Areas, and Buffer Zones. A comparison of the trends of deforestation inside and outside of those protected areas has been conducted, as a proxy measure of their effectiveness in tackling deforestation, particularly in connection with the increasing production of charcoal at a national level (Sedano et al. 2016).

Quality control
Quality control measures were adopted to reduce the random error as much as possible (e.g., the misinterpretation of land cover or involuntary wrong input in the survey form) and systematic bias (e.g., users overestimating a particular land cover class) due to user interpretation. In particular, 5% of the total plots (selected randomly) were duplicated and reassessed by two different users. After each of the assessment phases, the duplicated plots were compared with the corresponding quality control plots and reassessed if inconsistencies were found. If recognizable systematic error patterns emerged, a full review of plots was recommended to the user(s). The duplicate plots used for quality control were not considered in the calculations to estimate land cover areas. The sampling error of the proportion of a given class i of land cover was calculated from the proportion of plots assessed with CE falling in the land cover i weighted by the total number of plots assessed. The confidence interval related to the area of class i is then calculated as follows: 1 3 where CI stands for Confidence Interval, Atot is the total area of the assessment (Mozambique: 801,537 km 2 ), Pi for the proportion of the i class in the total plots assessed, (1-Pi) as the proportion of plots not classified as i in the total plots assessed, n is the total number of plots assessed, and Ai is the area of the class i.
Data quality was improved through a semi-automated cleaning procedure (Bastin et al. 2017) which highlighted potential inconsistencies (e.g., presence of tree cover with no trees recorded, high tree cover with low tree count) that needed to manually be reassessed. The identification of inconsistencies was based on the structure of the survey created for the assessment to make possible errors more evident through the creation of incompatibilities easy to isolate through an automated routine.

Land cover change
In Mozambique, a change in land cover occurred in all the analyzed classes from 2001 to 2016 (Table 3 and Table 4). Analyzing their initial extent, forest and grassland/savannah areas decreased by 2.7% and 5.7%, respectively, and 6.2% of the forest cover existing in 2001 has been lost. On the other hand, some forest gain due to the encroachment of grassland/savannah has been observed. Wetlands also decreased by 1.9%. The increase of Cropland (18.2%) was mainly at the expense of the forest cover, while the expansion of Settlement (21.1%) was at the detriment of cropland and forest areas. When looking at land cover change data aggregated at the province level (Fig. 6), the same trend characterized at the national scale is shown by almost all provinces with a reduction of forest and grassland/savannah and an expansion of settlement and cropland. The reduction of the forest is in large part found in the provinces of Nampula and Zambezia. Only in Gaza and Niassa forest increased due to encroachment of the grassland/savannah. The land cover analysis at the ecoregion level showed that the greatest change occurred in the Miombo woodland and coastal areas with a significant reduction of the forest and grassland/savannah and an increase of cropland (Fig. 7). Settlement increased mainly in the coastal area. Mopane woodlands showed an increase of forest cover due to tree and shrub encroachment in plots previously classified as grassland/savannah protected areas, instead, had a more conservative effect on land cover changes (Table 5). In 2001, the 11.1 Mha of protected areas were mostly covered by forest (67.60%) and grassland/savannah (25.85%). Cropland, wetland, settlement, and other land covered 3.39%, 2.11%, 0.12%, and 0.93%, respectively, of the land cover. In 2016, land cover changes remained below the threshold of 1.00% (+ 0.78%, − 0.09%, + 0.06%, + 0.03%) for cropland, wetland, settlement, and other land. However, forest and grassland/    savannah areas have changed, with an increase in the forest surface of 1.58% and a decrease in the grassland of 2.36%.

Uncertainty analysis
The uncertainty analysis performed on the IPCC classification data highlighted very low error levels (approx. 3%) for the most represented land cover classes (forest, cropland, and grassland/savannah). The error was above 10% only for the least common land cover classes (settlement and other land); still, in these cases, the corresponding confidence interval is narrow, due to the small area covered by these classes.

Discussion
The point-based method applied in our study has proven to be an effective tool to evaluate the LCC dynamics in Mozambique over 15 years, with a maximum error of 3% for the main land cover classes. The analysis showed that Mozambique has gone through significant changes in land cover in recent years (Table 4). Forest still represents the most widespread land cover in the country (58.2%); however, about 3.7 Mha of the original forests were lost, and an overall reduction of 1.3 Mha (2.7%) in forest cover has occurred, mainly due to the conversion of forest land to cropland (Table 3). The results confirm the findings of other assessments conducted on a national scale that identified Mozambique as one of the countries worldwide most affected by deforestation (Bodart et al. 2013;FAO 2020). Moreover, deforestation is also often accompanied by land degradation caused by fire ( Fig. 5), which has been recognized as the main disturbance factor for forests in eastern and southern Africa (Artes et al. 2019;Giglio et al. 2018). Furthermore, deforestation is not occurring evenly throughout the country and higher rates were recorded in the provinces of Nampula and Zambezia (Fig. 6). This result can be explained considering a combination of factors. These two provinces are located in coastal areas where the landscape is mainly dominated by the Miombo woodland. This specific ecoregion is characterized by higher precipitation than the Mopane woodland and is more suitable for shifting cultivation and the establishment of more permanent crop types. Moreover, in Mozambique, the long civil war (1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992) and recurrent droughts in the hinterlands have resulted in increased migration of the population toward urban and coastal areas, particularly Nampula and Zambezia, where the population has almost doubled over the past 10 years (Instituto Nacional de Estatística 2020). In turn, the growing population has increased the demand for food production and the need for large-scale agricultural investments (LAIs). Recent studies have highlighted how LAIs indirectly and directly cause deforestation on a much larger spatial and temporal scale than shifting cultivation (Di Matteo and Schoneveld 2016). Moreover, the conversion of large portions of land to intensive farming areas results in reduced access to rivers and fertile soils for local small-scale farmers (Zaehringer et al. 2018), triggering a feedback mechanism that increases the dependence of the population on large-scale production due to the lower feasibility of house farming. At the national level, these dynamics highlight how synergetic land and resource grabbing can be fostered by socio-environmental crises and in turn lead to new, more severe ones (Bruna 2019;Clements and Fernandes 2013) through the loss of valuable ecosystem services provided by the large portions of land converted. All these pressure factors have a particularly heavy impact on the Miombo woodland, which has been identified as one of the five global wilderness areas that should be prioritized for conservation (Mittermeier et al. 2003). On the other hand, Gaza and Niassa were the only two provinces characterized by an increase in forest cover (Fig. 6). This result can be explained considering that the two provinces have a lower population density (Instituto Nacional de Estatística 2020) and host the largest protected areas in the country, respectively, Bahnine, Zinave, and Limpopo National Parks in Gaza and the Niassa National Reserve in Niassa. The presence of protected areas combined with lower population pressure, as well as the higher incidence of small-scale agriculture, seems to have had a positive effect in preventing the conversion of forests into croplands that is characterizing the rest of the country (Table 5.). The increase in the forest cover in these two provinces was mainly due to the woody encroachment of grassland/savannah that is globally characterizing this type of ecosystem (Devine et al. 2017;Marston et al. 2017;Stevens et al. 2016). It has been hypothesized that woody encroachment is driven by the increase of CO2 in the atmosphere, and on a local scale, the encroachment is favored by both the large presence of protected areas in Gaza and Niassa, which lower the level of disturbances, and the low elephant population density due to intense poaching.

Conclusion
The results obtained in this study are in line with the existing literature on LCC dynamics in Mozambique and consistent with the socio-economic processes that have characterized the country in recent decades. The augmented visual interpretation of VHR satellite images using the Collect Earth software proved to be a reliable methodology to assess land cover changes, providing a comprehensive dataset that helps understand the drivers and causes of change and elaborate future action plans. The approach allows for a more precise evaluation of each plot, but it is also more time-consuming and labor-intensive than that of land cover mapping based on the automatic classification of satellite imagery. However, CE augmented visual interpretation can effectively complement and integrate the LCC analyses conducted using wallto-wall mapping to support national land assessment and forest inventories and provide training data for environmental modeling.
Acknowledgements This work has been conducted within the framework of SECOSUD II Project [Conservation and equitable use of biological diversity in the SADC region: from geographic information system (GIS) to Spatial Systemic Decision Support System (SSDSS)-AID: 9695] founded by the Italian Agency for Development Cooperation. The contents of this publication are the exclusive responsibility of the authors and do not necessarily represent the point of view of the Agency.
Funding Open access funding provided by Università degli Studi di Roma La Sapienza within the CRUI-CARE Agreement.

Data availability
The data that support the findings of this study are available from the corresponding author, Silvio Cianciullo, upon request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.