An Analysis of Urban Land Use / Land Cover Changes in Blantyre City, Southern Malawi (1994–2018)

: Rapid and unplanned urban growth has adverse environmental and social consequences. This is prominent in sub-Saharan Africa where the urbanisation rate is high and characterised by the proliferation of informal settlements. It is, therefore, crucial that urban land use / land cover (LULC) changes be investigated in order to enhance e ﬀ ective planning and sustainable growth. In this paper, the spatial and temporal LULC changes in Blantyre city were studied using the integration of remotely sensed Landsat imageries of 1994, 2007 and 2018, and a geographic information system (GIS). The supervised classiﬁcation method using the support vector machine algorithm was applied to generate the LULC maps. The study also analysed the transition matrices derived from the classiﬁed map to identify prominent processes of changes for planning prioritisation. The results showed that the built-up class, which included urban structures such as residential, industrial, commercial and public installations, increased in the 24-year study period. On the contrary, bare land, which included vacant lands, open spaces with little or no vegetation, hilly clear-cut areas and other fallow land, declined over the study period. This was also the case with the vegetation class (i.e., forests, parks, permanent tree-covered areas and shrubs). The post-classiﬁcation results revealed that the LULC changes during the second period (2007–2018) were faster compared to the ﬁrst period (1994–2007). Furthermore, the results revealed that the increase in built-up areas systematically targeted the bare land and avoided the vegetated areas, and that the vegetated areas were systematically cleared to bare land during the study period (1994–2018). The ﬁndings of this study have revealed the pressure of human activities on the land and natural environment in Blantyre and provided the basis for sustainable urban planning and development in Blantyre city.


Introduction
Studying land use/land cover (LULC) is vital for enhancing our understanding of global environmental change and sustainability [1]. In recent times, LULC changes have remarkably intensified due to increased anthropogenic processes such as urbanisation [2]. In 2018, 55 per cent of the world's population lived in urban areas, a proportion that is anticipated to reach 68 per cent by 2050. Almost 90 per cent of this expected growth will occur in Asia and Africa, especially in medium and small-sized cities [3]. In Africa alone, the urban population was 42.9 per cent in 2018 and is projected to reach 56 per cent by 2050 [4]. Moreover, in sub-Saharan Africa, where over 200 million people in urban areas reside in informal settlement, a higher urbanisation rate at 4.5 per cent annually has been reported [5]. The growth of urban areas has a significant influence on the global and regional environments, including LULC changes, and has implications for environmental, social and economic sustainability [6]. Rapid and unplanned urbanisation has had dire consequences, such as a reduction in vegetation cover and loss of biodiversity, as habitats for species become fragmented through the conversion of land for infrastructure development [7].
Like other anthropogenic-environmental interactions, urban LULC changes are due to a myriad of factors, as no single factor can account for these changes. The interactions are different in every region, but most scholars agree that most LULC changes are influenced by specific economic, demographic, socio-political and environmental conditions [8,9]. These factors are usually interrelated. For example, the economic and social advantages found in urban areas compared to rural areas attract many people to cities, leading to rapid population growth that contributes to the over-exploitation of natural resources for settlement and livelihoods.
Malawi is among the world's least developed and one of the most densely populated countries in Africa. Like other developing countries in sub-Saharan Africa, Malawi has been experiencing progressive urban growth since it attained its independence on 6th July 1964. Malawi's system of governance changed from a one-party rule dictatorship  to a multiparty democratic system from 1994 to the present. During the democratic rule, the urban population increased 2-fold from 1,095,419 in 1998 to 2,115,867 in 2018 in its four main cities of Lilongwe, Blantyre, Mzuzu and Zomba, where over 70 per cent lives in informal settlements [10,11]. Presently, Malawi is ranked among the top 10 countries in the world projected to have the largest population increase in both its rural and urban areas [4]. The urban population as a percentage of the national population was 16 per cent in 2018 and it is projected to rise to 30 per cent by 2030 and 50 per cent by 2050 [10][11][12].
Despite that several studies have been done on LULC changes in Malawi [13][14][15][16][17][18], comprehensive studies on urban LULC change in Malawi's cities remain scarce, as such, the understanding of urbanisation is primarily based on population figures only. This inadequacy of LULC change information constrains effective economic and environmental planning, resulting in uninformed policy decisions [16]. This is particularly prominent in low-income countries like Malawi, which commits most of its resources to address urgent needs such as poverty reduction at the expense of maintaining a vibrant LULC system [16].
Over the past three decades, advances in remote sensing technologies have expedited LULC change studies. By obtaining satellite imagery over a period of time, remote sensing methods can be utilised to analyse historical LULC changes [19]. Regardless of some shortfalls, such as spatial and spectral confusion of the urban areas [20], remote sensing remains a reliable source of data to support LULC studies [21]. Therefore, studies of urban LULC changes using remote sensing data, such as Landsat, are essential for land management and urban land use planning, especially in developing countries where they can provide fundamental and cost-effective information that is not available from other sources [22,23].
Although cities in Malawi actively engage in planning-including the development of master plans that provide an overview of spatial and infrastructural intentions, strategic plans that outline the broad ways in which thematic issues are to be addressed and investment plans that list priority infrastructure-the plans are barely recognised by the public in urban jurisdictions and often lack effective implementation mechanisms. This is usually blamed on a lack of resources [24]. Additionally, the existence of multiple institutions in the land administration in cities causes the disorder, which arises from unclear authority and mandates over land including development control [24]. Furthermore, cities in Malawi neither exercise control over key sectors (such as utility providers) nor have the authority to make other agencies align their plans to citywide coordination, and thus are unable to facilitate cross-sectoral planning [12]. As a result, cities continue to experience various challenges, such as rapid urban growth, which can lead to irreversible LULC changes, the proliferation of unplanned settlements and environmental degradation.
In Blantyre city, the second-largest city in Malawi by population and size, the urban LULC situation is unclear and has not yet been quantified. There is also not enough empirical data on urban LULC changes over time. Therefore, the aim of this study is to analyse urban LULC changes from 1994 to 2018 in Blantyre city using remotely sensed Landsat data in order to support sustainable urban planning. The year of 1994 was chosen as a start date because it is when Malawi changed its governance system from a one-party system to a multiparty system . The year of 2018 was chosen as  the end date because it is when the field verification was taken, while the year of 2007 was chosen as  an intermediate date to illustrate the rates of change. The study also identified random and systematic transitions derived from the classified maps to focus on prominent signals of change in the study area. The study period 1994-2018 was chosen in order to understand how the democratic governance systems influenced the LULC changes in urban environments of Blantyre city. The significance of the study is threefold: (1) the study will reveal the underlying human processes in the urban environment and their interactions; (2) the information generated can assist with managing the pressures of human activities and urban developments on the land, and (3) the results of this study can be used as baseline information to determine future urban land use and for setting policy priorities to promote inclusive and equitable urban development. Ultimately, this will help to realise well-balanced urban growth for citizens and the environment in Blantyre city.

Study Area
Blantyre is the oldest and second-largest city in Malawi. It was established by the Scottish missionaries in the 1870s and declared a planning area in 1897. It is the main commercial city hosting most of the private sector headquarters in the country. The city lies within the latitude 15 • 47 10" S and longitude 35 • 00 20.10" E, as shown in Figure 1. It covers an area of 240 km 2 [11]. The city is located in rugged terrain with multiple hills and river systems which have a direct effect on the development of the city [25]. Its topography ranges from an elevation of about 780 to 1612 m above sea level (m.a.s.l), as shown in Figure 1. The main mountains, such as Ndirande, Soche, Bangwe and Michiru, are also designated forest reserves [25].
In 1966, Blantyre was the prime city in Malawi with 109,461 people, seconded by Zomba, the then capital city, with only 19,666 inhabitants. With the transfer of the capital city from Zomba to Lilongwe, as part of a deliberate government policy to redistribute the urban population and implement spatially balanced development across the country [26], Lilongwe started to grow rapidly. The population grew from a mere 19,425 in 1966 to nearly 100,000 in a decade. By the year 2008, Lilongwe had overtaken Blantyre by about 7700 people. Lilongwe city continued to grow ahead of Blantyre city in 2018, registering a population of 989,318 against 800,852 in Blantyre [11]. Blantyre city currently hosts almost 5 per cent of the national population and has the highest population density in the country with 3328 people/km 2 [11]. Over 91 per cent of the city residents are below the age of 45 years, and over 65 per cent reside in the informal settlements [11].

Remote Sensing Data
The remotely sensed data used for this study were Landsat Thematic Mapper (TM) and Landsat Operational Land Imager (OLI) with 30 m spatial resolution and a 16 day repeat cycle [27]. The cloud-free Landsat images covering the study area (path 167, row 71) were downloaded from the United States Geological Survey (USGS) Earth Explorer https://earthexplorer.usgs.gov/, as shown in Table 1. The acquisition quality of these images was high (09, meaning no quality issues/errors detected) [28]. The data were acquired during the dry season, which starts from May to October, to best distinguish the spectral signature of different land cover types, and near-anniversary dates were chosen for consistency between the two time periods. Table 1. Landsat data used for the analysis of land use/land cover (LULC) changes in the study area.

Image Processing
The Landsat images were an L1T product (systematically, geometrically and topographically corrected). The study area was clipped as a vector file onto the Landsat image in ArcGIS Pro software (ESRI, Redland, CA, USA). In order to reduce confusion between spectral features and improve overall image classification, the radiometric and atmospheric corrections were applied using the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) Model in ENVI 5.2 software (Harris Geospatial, Broomfield, CA, USA) prior to image classification [29].

Image Classification
The support vector machine (SVM) classifier in the supervised classification scheme in ArcGIS Pro software (ESRI, Redland, CA, USA) was employed to derive four major LULC classes, namely bare land, built-up area, vegetation and water, as shown in Table 2. The classes were adopted from [30] with minor modifications to suit the study area. The training samples for each LULC class were determined by comparing the false colour composites using spectral bands 5, 4, 3, 2 and 1 for Landsat 5 TM and 6, 5, 4, 3 and 2 for Landsat 8 OLI as shown in Figure 2. Reference was also made to the Google Earth archived images and the ground control points collected during the study area field visit in August 2018 when collecting the training samples.

Accuracy Assessment
Accuracy assessment is a prerequisite step in the LULC classification process. This aims to quantitatively determine how effectively pixels are grouped into the correct feature classes in the area under investigation [31]. Accuracy assessment, therefore, compares the classified results to geographically referenced data that are presumed to be true. This comparison is typically achieved in a statistically rigorous fashion using error matrices. For this study, the reference data for 1994, 2007 and 2018 maps were collected from Google Earth image archives for the respective years. A set of at least 100 reference points were collected using stratified random sampling based on the sizes of the LULC classes for 1994, 2007 and 2018 classified images. The producer's, user's and overall accuracies, as well as a non-parametric Kappa index, were calculated [32,33]. The Kappa coefficient was computed using Equation (1).
where N is the total number of observations in the error matrix, r is the number of rows and columns in the error matrix, X ii is the number of observations in row i and column i (i.e., the diagonal elements), X +i is the marginal total of row I, X i+ is the marginal total of column i. A Kappa value greater than 0.80 indicates excellent agreement, while a value between 0.4 and 0.80 indicates moderate agreement and a value less than or equal to 0.4 indicates poor agreement between classification categories [32].

Annual Rates of Change
For an improved description of LULC changes, we calculated the annual rate of change. This process measures the amount of LULC change per year and was calculated according to the approach proposed by [34] and described in Equation (2).
where R is the rate of change, A 1 and A 2 are the area in km 2 at the beginning and end of the analysis period, and t 1 and t 2 correspond to the time in years from start to finish, respectively.

Land Use/Land Cover Change Detection
LULC change detection was carried out by a post-classification comparison of bi-temporal maps. This is perhaps the most common approach for change detection [35] and has been successfully used by many studies, such as [1,19,[36][37][38][39][40][41]. The cross-tabulation matrix for the LULC changes for the first period (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007), the second period (2007-2018) and the overall period (1994-2018) were generated. From the transitional matrices, the diagonal values in each matrix indicate persistency between LULC classes from initial time t 1 , and the later time t 2 , while the off-diagonal entries indicate the transition from one LULC class to the other. The gain was also calculated through the difference between the total value of each LULC class for the later period P + j , and persistency P jj , while the loss was the difference between the total for the initial time P i+ , and persistency P ii . On the other hand, swapping tendencies were calculated as two times the minimum value of the gains and losses for each LULC class [42]. Furthermore, the concept of systematic and random change was applied to the transition matrix in order to detect significant inter-category transitions [42]. An abrupt or unique process of LULC change is regarded as a random transition while that characterised by a permanent, stable or common process of change is regarded as a systematic transition. In order to determine whether a transition is random or systematic, the expected value was compared with the observed value in the matrix. A transition is said to be random if the difference between the observed and expected values is zero, while it is a systematic transition if it is not near zero. The expected values in terms of gain, G ij , and loss, L ij , were calculated using Equations (3) and (4) for (∀i j).
where P + j − P jj is the observed total gross gain of category j, P i + is the size of category i in the start date, 1994 in our case. J i=1, i j P j+ is the sum of the sizes of all categories, excluding categories of j, in the start date (1994).
where (P i+ − P ii ) is the observed total gross loss of category i, P + j is the size of category j in the last date, 2018 in our case. J j=1, j j P + j is the sum of all sizes of all categories, excluding j, in the later date (2018).
Transitions with observed values that were larger than expected values by 1 per cent point under both random processes of gain and losses were identified as systematic transitions [43,44].

Land Use/Land Cover Classification
The spatiotemporal patterns of LULC classes considered in this study for 1994, 2007 and 2018 derived from the Landsat images using the SVM classifier are shown in Figure 3. Table 3 Table 4 shows the results of the accuracy assessment for the 1994, 2007 and 2018 classified LULC maps. The overall accuracies were 89.5, 87.5 and 86.6 per cent for 1994, 2007 and 2018 images, respectively. These accuracies met the minimum overall accuracy set out by USGS (i.e., greater than 85 per cent). Hence, the classified results can be used as a data source for post-classification comparisons and further analyses. The Kappa coefficient of 0.79, 0.76 and 0.75 for the respective maps showed a good agreement between the classified map and reference data [33].  Table 5 shows the LULC transition matrix derived from the classified maps for the first period (1994-2007), the second period (2007-2018) and the overall period (1994-2018). In Table 5a-c, the proportions of the LULC classes that were persistent are provided in the diagonal entries of each matrix. During the first period, as shown in Table 5a, only 18 per cent of the landscape changed from one class to the other compared to about 24 per cent in the later period (2007-2018) as shown in Table 5b, affirming that changes were more rapid in the later period than the former. Further analysis of the LULC transition matrix for the entire period from 1994 to 2018 revealed that about 74 per cent of the landscape remained stable, meaning that only 26 per cent of the study area exhibited the transitions from one LULC class to a different LULC class.

Post-Classification Change Detection
The matrices were also analysed to understand the transition budget by calculating the gains, losses, total change, swap and net change as shown in Table 6. These calculations suggest the degree of interactions among the LULC classes in the landscape. The analysis of Table 6 revealed that the higher proportion of bare land and vegetation classes experienced a swap type of change throughout the study period. For example, 87 and 71 per cent of the total changes for the bare land class swapped locations during the first and second intervals, respectively. Similarly, almost 91 and 74 per cent of the total change for the vegetation class also swapped locations during the respective intervals. The swapping tendencies for bare land and vegetation classes could be attributed to the simultaneous reforestation and deforestation activities within the landscape. Unlike the bare land and vegetation classes, the built-up class experienced a larger proportion of the net change during all the intervals, taking up 55 and 70 per cent of the total change in the class during the first and second study interval, respectively.
The losses and the resulting swap proportions for the built-up class would not be expected but it could have resulted from the errors of commission and omission during the classification process. The error of commission and omission are the complement of the user's and the producer's accuracies, respectively, as shown in Table 4. In addition, the Landsat spatial resolution of 30 m and the edge effects in the determination of the changes could also result in some misclassifications along with the difficulties with spectral confusion between the built-up and bare land classes in the urban environment [19][20][21][22][23].

Detection of Random and Systematic Transitions
This study further analysed the transition matrix for 1994 and 2018 in order to identify the most significant changes among the LULC classes. The expected gains under the random process of gain were calculated using Equation (3). The expected gains and the differences between the observed and expected proportions are given in Table 7a,b. In Table 7b, numbers closer to zero indicate a random transition between LULC classes, while numbers further from zero implied systematic transition [42]. In this study, only transitions that were larger than 1 per cent point were identified as systematic transitions [43,44]. The differences between the observed and expected gains under the random processes of gain from the bare land class to the built-up class (1.20) and the vegetation class to the bare land class (1.64) was positive and above 1 per cent point, meaning they occurred systematically rather than randomly.
The differences between the observed and expected gains for vegetation to the built-up area (−1.19) and built-up to bare land (−1.61) were negative and above 1 per cent point. This implies that when vegetation gains, new vegetation systematically avoids gaining from built-up areas. Similarly, when bare land gains, it systematically avoids gaining from the built-up class [42].
The equivalent relationship between the different LULC classes, which examines the actual and expected losses under the random process of loss, is shown in Table 8a, b. The difference between observed and expected losses under the random process of loss for the bare land class to the built-up class and the vegetation class to the bare land class was 1.01 and 1.09 per cent points, respectively. This infers that the bare land class and the vegetation class lost systematically to the built-up and the bare land classes, respectively. Likewise, the difference between the observed and expected losses for vegetation to the built-up area (−1.10) was above 1 per cent point but negative, implying that vegetation systematically avoided losing to the built-up class [42].
In a nutshell, there has been a systematic conversion of the bare land class to the built-up class as well as a systematic degradation of the vegetation class to the bare land class in Blantyre city over the 24-year study period. This was demonstrated by the concurrent incidences of the systematic gains and losses [44]. This means that 13.66 per cent transition of bare land to the built-up class and the degradation of 7.31 per cent of vegetation to bare land were due to a systematic process of change, as shown in Table 5c. The LULC map showing systematic changes from bare land to built-up area and vegetation to bare land, as well as persistence and other random changes, is shown in Figure 4.

Discussion
The results indicated that there were increased anthropogenic-induced urban LULC changes in Blantyre city over the past 24 years. This was substantiated by the observed accelerated increase in built-up area from 4.1 per cent during the initial period to almost 7 per cent in the later period, with an overall annual change at 5.3 per cent. This is similar to the accelerated declining rates observed for the vegetation and bare land classes (Section 3.1). The increase in the built-up class is comparable to other studies conducted in the sub-Saharan African cities, including the Dakar metropolitan area in Senegal, Nairobi city in Kenya, and Harare city in Zimbabwe, which experienced an increase in built-up areas at the annual rates of 9.6, 9.5 and 4.7 per cent, respectively, between the years 1990 and 2014 [45]. The observed increase in the built-up class and the decline in vegetation and bare land classes, respectively, have several implications for sustainable urban planning of the area.
Firstly, the observed annual growth rate of 5.3 per cent for the built-up class in the study area surpassed the urban population growth rate estimated at 2.3 per cent between 1998 and 2018 [10,11]. This means that the urban growth in Blantyre City, with a land consumption ratio of 2.3, is becoming more expansive than compact [46]. This kind of growth creates profound repercussions for environmental sustainability in the city and also prevents the city from enjoying high social interaction due to close integration of communities, and easy access to social-economic facilities [47]. Further to that, the excessive increase in the built-up area indicates an increment of more impervious surfaces in the city. The increase of impervious surfaces causes a decrease in the groundwater recharge as well as high reflection of solar radiation back to the atmosphere, hence contributing to environmental problems such as urban flooding and urban heat islands [48]. This development puts the city at risk and calls for better management to guarantee sustainable urban development.
Secondly, the vegetation loss observed in this study signifies the loss of green spaces (such as forests and parks) [49]. The vegetation loss results in declining ecosystem services, such as air and water purification, flood mitigation services, noise reduction and climate regulation, including urban cooling [50]. It also causes soil degradation [51], which leads to the formation of gullies and derelict landscapes. In addition, such losses increase residents' vulnerability to environmental stress due to the loss of non-material benefits that people obtain from ecosystems through spiritual enrichment, cognitive development, reflection, recreation and aesthetic experience, as well as their role in supporting knowledge systems, social relations and aesthetic values [52].
Lastly, the decline in bare land exerts enormous pressure on land suitable for urban development. Recently, over 10,000 families were reported to have built their houses illegally in areas such as wetlands, steep slopes and river/stream buffer zones in the city [53]. Such unplanned settlements are prone to multiple hazards such as floods and landslides, which have increased in frequency and intensity in recent times due to climate change and climate variability [54]. Similar situations were also observed in Mzuzu City, northern Malawi, where people have encroached into areas prone to multiple hazards [55].
Furthermore, the analysis of the transitional matrix derived from the 1994 and 2018 classified maps, as shown in Table 5c, has revealed the two-way systematic pathway of changes as illustrated in Figure 5.
These systematic processes suggest that the gain in built-up targets bare land and avoids vegetation. This can be explained by the fact that most land under vegetated areas, such as Ndirande, Soche and Kanjedza forest reserves, as well as the Mudi catchment area, are protected leaving bare land vulnerable to conversion to built-up. Additionally, most vegetated areas have unsuitable land for development, such as wetlands, steep slopes and river/stream buffer zones. Likewise, the gain in bare land systematically targets the vegetation class and avoids the built-up class. This systematic process has caused the loss of vegetative cover in the study area. This loss could be explained by anthropogenic activities, such as wood extraction for firewood, brick burning and lapses in management [56][57][58]. For example, the majority of Malawians have no access to a reliable energy source. The proportion of households with access to electricity in Malawi has increased from 8 per cent in 2010 to 11 per cent in 2017 [59]. Nevertheless, the most common source of cooking fuel in the country is firewood with 81 per cent, followed by charcoal (16 per cent), electricity (2 per cent) and crop residue at 1 per cent. Although 62.9 per cent of Blantyre city residences had access to electricity in 2017, only 15.5 per cent use it for cooking and heating [59]. While acknowledging that 90 per cent of the firewood supply in urban areas originates from rural areas, the remaining 10 per cent still comes from urban forests and other vegetation types [60] and may have contributed to the vegetation loss observed in this study. Secondly, brick burning also contributes to vegetation loss. Construction industries rely heavily on wood for burnt brick production, which remains the main walling construction material in Malawi. It is estimated that over 0.7 metric tons of firewood ae required to create 1000 bricks [61], leading to massive forest degradation [62]. Traditionally, bricks are usually made at the construction site because when made off-site, transportation increases their costs substantially [63]. As such, it is logical to attribute the production of burnt bricks as a possible cause of vegetation loss in Blantyre city during the study period. From a governance perspective, the authorities prior to the democratic era had protected forest reserves from poachers and squatters, however, this authoritative control was lost during the transitional period in the early 1990s. This was largely due to the austerity measures that were introduced under the structural adjustment program. These measures led to the reduced number of forest guards who were protecting forest reserves. This change, coupled with the severe economic downturn around the same time, induced people to overexploit forest resources to support their livelihoods [64]. By 1995, the Ndirande Mountain Forest Reserve was largely deforested [56], and other forest reserves in Blantyre, such as Kanjedza, Soche, Michiru and Bangwe, were not spared. Local residents also invaded the Mudi Dam catchment area owned by the Blantyre Water Board of the city, destroying its vegetation cover that resulted in the siltation of Mudi dam [65].

Limitations of the Study
This study used Landsat data with a spatial resolution of 30 m, meaning that the changes below this pixel size might have been missed during classification processes. Therefore, the use of high-resolution satellite data could have provided more detailed information about LULC changes with higher accuracy.

Conclusions
This study provides the details of urban LULC changes in Blantyre city from 1994 to 2018 characterised by an increase in built-up area and the decline of vegetation and bare land areas, respectively. The study further reveals a two-step temporal transition, firstly from the vegetation class to the bare land class, and from the bare land class to the built-up class. This clearly demonstrates the existence of inefficiencies in the management of urban growth in the city. Based on this information, urban development stakeholders can make policy and planning priorities to ensure sustainable urban development. In order to ensure sustainable urban development in Blantyre city, this study suggests that the authorities should expedite the allocation of all suitable land for development while safeguarding unauthorised development in risky areas. This study, further, calls for all planning authorities in Malawi at the city, district, regional and national levels to review physical plans more often and help to facilitate the timely supply of serviced land across the country.
Regarding the observed vegetation cover loss, some innovative policy measures, such as identifying crucial vegetation areas for protection from anthropogenic forces, could be implemented. There is also a need to enhance cooperation among urban stakeholders and local residents, including providing incentives that could encourage locals to conserve vegetation and allow natural regeneration on bare hills and other vegetation reserves within the city boundary.
Author Contributions: J.M., T.W. and R.A. designed the study, analysed the datasets, drafted the manuscript and addressed the reviewer's comments. All authors have read and agreed to the published version of the manuscript.
Funding: Japan International Cooperation Agency (JICA) provided financial help in ABE 4th Batch Scholarship for JM and in the filed survey for JM and TW. Part of the APC was covered by JICA.