Assessment of land use change in the Thuma forest reserve region of Malawi, Africa

A study was conducted in Thuma area in central Malawi to quantify contemporary land cover and to explore the degree of land use change in the Thuma forest reserve area of Malawi by analysing and comparing satellite-derived land cover maps from 1997, 2007 and 2017. The study was carried out using Remote Sensing and Geographic Information System (GIS), focusing on analysis of Landsat 5 ETM and Landsat 8 ORI/TIRS satellite images. The classification was conducted for the following distinct classes; closed forest, open forest, shrubland, savanna grassland, agriculture fields, and water. The analysis revealed that closed forest diminished from 19% in 1997 to 10% in 2007 to 6% in 2017. Open forest reduced from 30% to 21% from 1997 to 2007 but increased to 22% in 2017. Agriculture area almost doubled from 37% in 1997 to 64% in 2017. The actual area from 1997 to 2017, shows that closed forest has reduced from 7,000 ha to 3,000 ha while open forest from 12,900 ha to 7800 ha. Savanna grassland has doubled from 5,900 ha to 13,000 ha. However, future studies should use modern satellites such as Sentinel and Landsat 9 for improved quantification of changes. The findings show that even the protected forest reserve (previously dominated by closed forest) is not fully protected from deforestation by local communities. Government and other stakeholders should devise measures to meet the needs of the surrounding communities and the ecological/biophysical needs of the reserves. Based on this study, issues of re-demarcation of the forest reserve and the accessed area should also be explored. This study serves as a reference for the management of Thuma Forest Reserve as a refuge for natural tree species, rivers that harbour endemic fish species (Opsaridium microlepis and Opsaridium microcephalis) and the sustainable management of endangered elephants in the reserve.


Introduction
Land cover and land use mapping have become crucial components in the assessment and management of natural resources. Assessment of land cover dynamics is an essential requirement for the sustainable management of the environment and represents a fundamental baseline through which to support activities linked to environmental studies (Tiwari and Khanduri 2011). An accurate and updated land cover database is multipurpose in its utility and can support the assessment of landscape dynamics. This can be used to facilitate the planning of environmental restoration activities (Latham 2013).
Most protected forest reserves have failed to meet the expectation of increasing forest cover due to anthropogenic interferences (Kusimi 2015). Detection of land cover changes as a result of land use modification especially those changes caused by humans is paramount in policy decision making for sustainable management of natural resources (Estoque and Murayama 2015). Current studies in most sub-Saharan countries show that there is a general reduction in forest cover in most protected forest reserves due to expansion of agricultural land or mere tree removal for other amenities (Damnyag et  The area was physically delineated to include Thuma Forest Reserve and its surrounding area based on the hydrology network of the area, especially focusing on the small tributaries around the forest (figure 1). This area supports livelihoods in Salima, Dowa, Lilongwe and Dedza districts. It is a source of several small rivers, source of energy (firewood and charcoal), habitat for a range of flora and fauna, including elephants. Communities around Thuma area practice subsistence and commercial agriculture systems. Those involved in subsistence mainly grow maize, groundnuts, beans and vegetables while those practicing commercial agriculture grow tobacco which requires firewood for curing (Crops Department Report 2017).

Work flow
The procedure, as shown in figure 2, involved imagery acquisition, pre-processing to remove atmospheric errors, to get a clear view of the satellite images, enhancement and interpretation of the images. Field validation was done at the study site to obtain photos and descriptions of the sampled points. Unsupervised and supervised classification were performed on the images based on ground truth data. Error assessment was done to determine the level of accuracy of the classification exercise and finally, imagery change detection was performed.
2.3. Data sourcing 2.3.1. Satellite images used Landsat imagery data (figure 3) was used for land cover classification because it is free, has a long history (since 1972 to date) and higher frequency of archives with an average temporal resolution of 16 days (Manandhar et al 2009). For the year 1997 and 2007, Landsat images from Landsat 5 Thematic Mapper were used, which was operational from 1984 to 2013. Landsat 5 was used because Landsat 7 went into permanent malfunctioning on 31 May 2003 due to a fault in the satellite's Scan Line Corrector (SLC) (United States Geological Survey online 2021). For 2017, Landsat 8 Operational Land Imager and Thermal Infrared Sensor were used. Landsat 8 has more bands than Landsat 5 because it carries two sensors. It has additional band 1, deep blue visible channel for viewing water and coastal areas and band 9 for detection of Cirrus clouds (Mwaniki et al 2015). The images were downloaded from the United States Geological Survey website (GloVis). We specified the time of the imagery concentration on summer months because that is the time to obtain clear images for Malawi (Gumma et al 2019). Cloud cover was specified from 0 to 10%.
For this image analysis, seven bands from Landsat 5 and Landsat 8 were used. Band 1 (Blue green) differentiates soil and rock surfaces from vegetation. Band 2 (Green) covers the green reflectance peak from leaf surfaces near the visible spectrum (Lillesand et al 2015). Band 3, is the red band chlorophyll absorption band which helps in distinguishing vegetation and monitoring vegetation health. In this study, it was useful in distinguishing between closed forest, open forest, shrubland and grassland. The Near Infra-Red (NIR) Band 4 is good for differentiating vegetation, soil and water. To view false colour in Landsat 5, a combination of bands 4,3,2 was used. This means that to view vegetation, band 4 (NIR) was used basing on its wavelength which has  high reflectance at 0.7 to 0.9 micrometres. In Landsat 8, band 5 (NIR) which ranges from 0.851 to 0.879 micrometres (figure 5) was used to view vegetation. Different types of vegetation possess different types of spectral signatures. Bands 5,4,3 were used to view false colour in Landsat 8. Vegetation is viewed in red in false colour because green and blue are absorbed and only red is reflected by vegetation. Vegetation has a unique spectral signature which enables it to be distinguished readily from other types of land cover in the near-infrared. So, in vegetation mapping, this was mainly used to distinguish between closed forest, open forest and shrubland, agriculture and savanna. Water absorbs all colours, so the rivers were appearing as dark blue.

Field validation
Field validation process is done to improve the quality and accuracy of the data base (Latham 2013). In this study, field validation was done to obtain ground truth data to be used in image training for classification. Ground truth points were sampled from a polygon of the study area. The validation samples were collected for the year 2017. The samples were randomly selected in a Stratified Random Sampling based on accessibility (along the roads/paths) of the study area and location homogeneity (Gumma et al 2019). They were generated in ArcMap and were converted to GPX format and uploaded in a Geographic Positioning System (GPS) device which was then used to locate the sampled points on the ground. A total of 60 points (See appendix A) were selected. A deliberate effort was made to sample areas where it was difficult to distinguish the land cover classes in the preliminary interpretation, especially between savanna grassland and agriculture fields or open forest and shrubland. Four photographs, in the directions of North, East, South and West were taken at each point and dominant land cover was recorded (Manandhar et al 2009). Some examples of the selected points are shown in table 1.
2.4. Image pre-processing 2.4.1. Atmospheric correction Atmospheric correction was done to influence the pixel values of the image in order to adjust the values to compensate for the atmospheric effect (Campbell and Wyne 2011). The energy that is captured by Landsat sensors is influenced by the Earth's atmosphere. These effects include scattering and absorption due to interactions of the electromagnetic radiation with atmospheric particle like aerosols, gases and water vapour. Hence the need to correct the atmospheric Landsat data to earth surface values (Young et al 2017). Atmospheric correction was done for only Landsat 5 data in Quantum GIS, using the SCAP plug in. The correction was done in all the seven bands (Sharmila 2013). There was no need to perform atmospheric correction for Landsat 8 data because its data already have in-built atmospheric corrections. At-satellite radiance was calculated from Digital Number (DN) values basing on the maximal and minimal spectral radiance value (or the gain and bias) for each band (Chander et al 2009). The following equation was used for atmospheric correction (Cui et al 2014).
Atmospheric correction, Where Lβsat = Spectral radiance at the sensor, DN cal =Digital numbers from a specific band, Gain β = Radiometric gain and Bias β =Radiometric bias

Cloud removal and masking
The 1997 Landsat 5 satellite image had clouds in southwest area. So, ArcMap spatial analyst tool was used to remove the clouds. An equal area was clipped out of all three images. Since forests are characterized by dark colour due to chlorophyll and shadows they cast in the visible band realm (Goward et al (1994); Huemmrich and Goward 1997), they are easily distinguished in image processing (Dodge and Bryant 1976). This forms the basis for setting boundaries for delineation of pixel for true colour representation in Landsat images. Pixel area delineation was done using histogram in the local image window. Due to their darker nature of forest images, they appeared at the base of the histogram. Cloud delineation was done by normalizing the temperature at each isometric zone. The details of this procedure are provided by Huang et al (2010). Then linear cloud boundaries were made to achieve non-infrared bands in the Landsat image to construct the spectral-temperature space. The shadows of the clouds were detected spectrally and geometrically. Each cloud pixel was detected with reference to location, solar illumination and cloud height. An appropriate lapse rate was used to estimate cloud height (Smithson et al 2008). Finally, a determination of a classification was made.
2.5. Image processing 2.5.1. Image enhancement The images were imported and opened in ArcMap where all the bands from 1 to 7 were stacked using the composite band tool. Image enhancement was performed on the images to improve visibility. In ArcMap the Dynamic Range Adjustment tool together with the stretched and red green and blue (RGB) composite tools were employed (Jia and Richards 2006). False colour image combination was used because it fits well with the target of viewing forest and other vegetated land covers. Because of the spectral signature of vegetation (Appendix B), the classes were distinguished visually using Near Infrared, Red and Green bands. Summarily, all the desired rigorous protocols/procedures were followed to identify the respective classes using the high-resolution images from Google Earth for the large volumes of data. The technique termed, Spectral Matching Technique has been fully described and used by Thenkabail et al (2007)

Unsupervised classification
Unsupervised classification was performed on the images so that land cover classes with similar spectral signatures could be easily identified. This process was run on a seven-band stack of all the images to gauge the performance of the software on the classification of the pixels depending on the specified number of classes without human intervention. This process helped getting a general understanding of the probable land cover classes in the area to be input as training sites in supervised classification (Liu et al 2003). The IsoCluster unsupervised method was used in ArcMap for computing the minimum Euclidean distance when assigning each candidate cell to a cluster using a modified iterative optimization clustering procedure, also called Migrating Means Technique. The algorithm separated all the cells into the user-specified number of distinct unimodal groups in the multidimensional space of the input bands (Lillesand et al 2015). The result from the unsupervised classification was nine distinct classes of which only two (water and closed forest) could relate to the situation on the ground.

Supervised classification
Supervised classification was conducted in order to extract quantitative information from remotely sensed Landsat images data. The goal of supervised classification was to assign each pixel in a study area to a class or category. In this classification, knowledge about the study area was obtained from field validation which was used to identify representative areas, or samples of each class (Feringa et al 2001). The training sites were determined by ground truth data, an average of twenty training areas were done per class. The training areas were larger for big classes like agriculture and small for small classes like shrubland and water. The classification was based on the following distinct land cover classes; closed forest, open forest, shrubland, savanna grassland, agriculture fields and water. In the present study, training sites for these classes were developed and run to classify the images of interest. A classification signature file was created for 2017 image which was then used as the basis for classifying 2007 and 1997 images. This was then used to run the maximum likelihood classification in ArcMap. The resulting product was a classified raster image with the five land cover classes for each year. The symbols and colours of the classes were edited and labelled appropriately. After observing the results of the accuracy assessment, training samples or classification parameters were adjusted to get a better result or more accurate classification (Herold et al 2008).
2.6. Error assessment Accuracy assessment was done to compare a classified image with the ground truth data to assess the reliability of the data (Congalton 1991). Classification error occurs when a raster pixel belonging to one land cover class is assigned to a different land cover class. Errors of omission occur when a feature is not included in the category being evaluated; errors of commission occur when a feature is incorrectly added in the category being evaluated. An error of commission in one category will be counted as an error in omission in another category (Feringa et al 2001). An error matrix compares the relationship between known reference data and the results of the classification, based on the classes. In this study, accuracy assessment was done for the 2017 classified raster image. A total of 180 points were randomly selected from the polygon of the study area using the sampling tool in ArcMap. This point layer was used to extract values from the classified 2017 raster.  , pivot table statistical tool was used to pivot the frequency table into an error matrix table. Calculations were done in the confusion matrix and it was established that the classification had an accuracy of 68%.

Change detection analysis
After developing a classified land cover map for each year in ArcMap, the Spatial Analyst Zonal Statistics Tool for cross tabulating areas was used to detect changes in land cover per class. The tool detects the exact area change These results were also presented as percentage changes per class per year. During field validation, only 60 points were sampled due to logistical hiccups and inaccessibility of some sites, but the ideal would be more than 50 points per land cover class (Congalton 1991). However, Google Earth image was used to validate inaccessible points to minimize errors. Finally, a land cover change study interval of 10 years was used. The confusion matrix was done for 2017 only because it was possible to compare the ground truth and classified land cover unlike the cases for 1997 and 2007. The accuracy result from the confusion matrix was 68% representing the moderate agreement between the reference data and the classified image (Jensen 2004).

Results
Premised on the purpose of this study, land cover classification for the study area was run for 1997, 2007 and 2017. The results are presented in form of maps, charts and tables. The overall land cover changes are shown in figures 6 and 7. Tables 2-5 show a detailed percentage change in ha of land cover from one class to another between the years. After the image analysis process and classification, the distinct land cover classes that could be clearly differentiated from the 30 m resolution Landsat images were; Closed Forest, Open Forest, shrubland, Savanna Grassland, Agriculture and Water. As it can be observed from the maps and graphs (figures 4-5), it is evident that generally closed forest area is dwindling from 19 percent in 1997 to 11 percent in 2007 to 6% in 2017.

Discussion
The study has revealed many changes in Thuma area. The findings may be partly explained by the theory of Agricultural Adjustment to Land Quality (AALQ) by Martha and Needle (1998). The theory proposes that people's search for land of good quality leads to the clearing of forests for agricultural fields (crop and livestock production). The theory further asserts that the shift from a piece of land to another potentially leads to the forest cover increase in the abandoned agricultural land. This results in the change from agricultural land to forest cover. There was a complex dynamism in terms of the changes in land cover.    The closed and open forests in Thuma area are under ecological threat. Even the most protected Thuma Forest Reserve (which used to be dominated by closed forest) face high tree losses as the cover class keeps diminishing across the years. This implies that there are illegal logging activities in the protected forest reserve. Tree loss in Thuma area should be viewed as a source of biodiversity loss in the area. This does not only apply to tree species but also fauna species in the area. There is a general increase in agricultural land from closed and open forest by the surrounding villagers (figures 6 and 7(A)). This should be a wakeup call to policy makers and the management of Thuma Forest Reserve that encroachment is increasing in Thuma area. This will have a direct effect on the management of wildlife especially elephants in Thuma Forest Reserve as the habitat changes. Habitat change directly modifies species range (Orlovsky et al 2006). It would be expected that cases of poaching would arise in the area due to encroachment. This study forms a source of planning tool for the four surrounding districts in their 5-10 year District Development Plans (DDPs). The level of tree cover loss and the need for agricultural land should be central in the development of the DDPs considering that Thuma area is the source of wood energy and source or transit area for rivers and streams for Dedza, Dowa, Lilongwe and Salima Districts. Knowledge on positive and negative response of the land scape is paramount for future predictions of improved land governance and ecological management (Lambin et al 2003).
The loss in shrubland could be attributed to the fact that tobacco-growing farmers around the catchment cut down the small sized trees for use in tobacco curing and drying (Manyanhaire and Kurangwa 2014). This can also be attributed to use of shrubs for burning of bricks for house construction. As the shrubland decreases, savanna grassland is increasing. This could indicate a shift in cover change due to intense removal of the shrubs. A transition of land cover change from forest to bushes and shrubland has been reported in West Africa (Kusimi 2015). Although 30.9% of agriculture land went to savanna grassland from 1997-2007, almost the same area (30.5%) was converted to agriculture from 2007-2017. This could be attributed to the population increase in the four districts surrounding Thuma area (National Statistics Office 2018). From 1997-2017, the population in Dedza almost doubled from 486,682 to 830,512. In Salima, the population doubled from 248,214 to 478,346. Lilongwe and Dowa registered increases from 440,471 to 989,318 and 411,387 to 772,569 respectively. In most African countries, increased population has led to increased farm land expansion which leads to encroachment into some protected forest area (Smith et al 2016). It is also observed that there is more opening of agriculture fields in the areas with more people i.e. major village areas. The actual changes in land cover are as a result of direct causes (opening new crop land, settlement expansion, activities associated with mining, opening of new roads) and indirect factors (government policies, population growth) (Geist and Lambin 2002). Our findings are   forest area decreased as 61% was converted to barren land from 1991 to 2015 while 2.8% of the agriculture land was turned to built-environment (Munthali et al 2019).
More than 97% of the Malawi population uses unsustainable and illegal sources of biomass for domestic cooking (figures 7(B) and (C)) and heating energy (Environmental Affairs Department 2017). This is a challenge that needs both short and long-term strategies. The dwindling forest areas in Thuma should provide policymakers and natural resource managers the need to formulate both social and biophysical/ecological remedies to control the changes. It is apparent that the communities around Thuma area depend on the forest reserve for daily energy needs. It appears the 'land sparing' approach, where land is demarcated for specific use is slowly becoming unsuccessful and possibly there is a need to pursue 'land sharing', which entails the integration of more than two interventions (land uses) on a parcel of land (Loconto et al 2018). This could be achieved by encouraging farmers to engage in large-scale agroforestry where both targeted and naturally regenerating tree species should be managed at farm level. These trees should be the source of fuelwood energy by the communities. The approach of mixing crops and trees called Taungya in Ghana has been recorded to have reduced communities' reliance and level of encroachment in forest reserves (Acheampong et al 2016, Akamani andHall 2019). This system stabilised the agricultural land cover such that from 2002 to 2017 the land lost from agriculture was equal to the land gained (Appiah et al 2021). Other options could revolve around subsidising of agricultural inputs to increase per unit production, strict bylaws, enforcement of national laws, and comanagement. Such approaches will provide time for the recovery of the deforested forest reserves. Skole et al (2021), estimated that forest deforestation in Malawi ranges from 22,410 to 38,937 ha y −1 .
Water class in this case represents surface water flowing in all the rivers and Linthipe and Linduzi rivers in the micro-catchments. Although, it was difficult to map the rivers accurately because their width was very small in the 30 m × 30 m pixel of Landsat image, the trend showed that the surface area of the rivers changed from 309.87 ha to 139.95 ha. This may have an effect on access to water by wildlife in Thuma area which may lead to migration of these game out of the area in search of water. Tree species that are endemic to water sources may also be reduced due to water loss. However, the change in the surface area of the rivers in the present study might be over or underestimated since some river portions are covered by trees, making it impossible for the satellite sensors to capture their exact spectral signature. Nonetheless, the change may be attributed to reduced cover that may have exacerbated soil erosion (that may have silted the rivers) or increased evaporation (Vargas and Omuto 2016).

Conclusion
Land cover and land use change assessment is crucial for sound decision-making at national and local level. The prominent land cover classes for Thuma forest reserve area (closed forest, open forest, shrubland, savanna grassland, agriculture and water) have shown changes from 1997, 2007 and 2017. The area under closed forest reduced from 7,000 ha to 3,000 ha while open forest shifted from 12,900 ha to 7,800 ha from 1997 to 2017 respectively. The area under savanna grassland has almost doubled from 5,900 ha to 13,000 ha. The area under agriculture has increased from 10,000 ha to 13,900 ha. By decadal percentage changes, closed forest changed from 19% in 1997 to 11% in 2007 to 6% in 2017, as agriculture area wobbled from 24% in 1997 to 21% in 2007 and 33% in 2017. The savanna area doubled from 13% in 1997 to 31% in 2017. The change from forest, shrubland, grassland to agriculture clearly shows ecosystem transition which has an impact on the biodiversity of both plants and animals. Future studies should integrate the changes with the associated drivers. The studies should also incorporate different data sets for ease of comparison. We recommend the use of Sentinel images which have a higher resolution of up to 10 m than Landsat 5 and 8. Sentinel images could help fully understand changes, especially water bodies which were not thorough in the present study. This study serves as a harbinger to government and relevant institutions for the need to put in place strategies that will buttress the successful management of Thuma forest reserve area. With the unique tree species diversity in the reserve and the fact that the reserve is currently used as a refuge for elephants, there is a need to ensure appropriate and timely remedies for the recession of the conversion of the reserve area into other land uses like farming. The reserve plays important role in maintaining the ecological status of important rivers, Linthipe and Linduzi. Linthipe River is a breeding niche for the endangered Opsaridium microlepis (locally called Mpasa) and Opsaridium microcephalis (locally called Sanjika) fish species. Based on this study, issues of re-demarcation of the area may also be explored. Programmes on forest landscape restoration through afforestation or regeneration should be considered to halt deforestation. This band is very sensitive to moisture and is therefore used to monitor vegetation and soil moisture. It is also good at differentiating between clouds and snow. Band 6 10.40-12.50 μm, thermal infrared This is a thermal band, it is used to measure surface temperature. This is primarily used for geological applications, but it is sometime used to measure plant heat stress. This is also used to differentiate clouds from bright soils since clouds tend to be very cold. One 7 2.08-2.35 μm midinfrared