Spatiotemporal Analysis of Land Cover and the Effects on Ecosystem Service Values in Rupandehi, Nepal from 2005 to 2020

: Land cover (LC) is a crucial parameter for studying environmental phenomena. Cutting ‐ edge technology such as remote sensing (RS) and cloud computing have made LC change mapping efficient. In this study, the LC of Rupandehi District of Nepal were mapped using a Landsat imagery and Random Forest (RF) classifier from 2005 to 2020 using Google Earth Engine (GEE) platform. GEE eases the way in extracting, analyzing, and performing different operations for the earth’s ob ‐ served data. Land cover classification, Centre of gravity (CoG), and their trajectories for all LC clas ‐ ses: agriculture, built ‐ up, water, forest, and barren area were extracted with five ‐ year intervals, along with their Ecosystem service values (ESV) to understand the load on the ecosystem. We also discussed the aspects and problems of the spatiotemporal analysis of developing regions. It was observed that the built ‐ up areas had been increasing over the years and more centered in between the two major cities. Other agriculture, water, and forest classes had been subjected to fluctuations with barren land in the decreasing trend. This alteration in the area of the LC classes also resulted in varying ESVs for individual land cover and total values for the years. The accuracy for the RF classifier was under substantial agreement for such fragmented LCs. Using LC, CoG, and ESV, the paper discusses the need for spatiotemporal analysis studies in Nepal to overcome the current lim ‐ itations and later expansion to other regions. Studies such as these help in implementing proper plans and strategies by district administration offices and local governmental bodies to stop the exploitation of resources.


Introduction
Land cover (LC) is an important parameter to study and track the local, regional, and global changes of the earth's surface. It has been a crucial variable for different studies such as climate change, food security, environmental studies, conservational strategies, hydrology, and landscape planning [1][2][3][4]. LC, being subjected to natural phenomena and anthropogenic causes of land use, is very dynamic [2,3]. The changes in the constitution of LC represent one of the most significant sources of system changes at local, provincial, and worldwide scales [5]. To accomplish the United Nations' Sustainable Development Goals (SDGs), timely and higher resolution LC data are critical [6].
The long-term earth observation data record is a tool for analyzing the land use and land cover changes (LCC) [7]. Remote Sensing has been a useful aid in studying the patterns of LC with a decreased cost and less effort [2]. The easily accessible and freely available satellite imageries have assisted in studying the LC patterns at a regional and global scale, which would have been difficult, if not impossible. Various global LC data have been produced over the years, such as GLC2000 [8], MCD12, and GLCNMO, GlobCover [9], IGBP DISCover [10], MODIS Collection 5 global land cover [11], Global Land Surface Satellite Global Land Cover (GLASS-GLC) [12], GLC30 [13], GLC10 [6], ESA CCI Land Cover time-series [8], ESRI 10 m [14], etc. While some researchers are concentrated on a local or regional level to map out specific phenomena, for example, meteorology [15], future prediction [16,17], urban heat islands [18], climate change [19], water [20], forest [21] and biodiversity [22], and agricultural monitoring [23], current studies are also focused on improvements such as finer resolution and classification procedures [1,[24][25][26].
Image Classification is the basis of LC, in which we assign LC classes to a multiband raster image. Image Classification in remote sensing is broadly categorized into the following three types: pixel-wise classification, sub-pixel-wise classification, and objectbased image classification. Pixel-wise classification treats each pixel as pure and is labeled as a single land use-LC type based on the spectral signature and their derivatives such as vegetation indices [27], which are further divided into an unsupervised and supervised classification. Unsupervised classification, which is based on the principle of creating clusters and assigning classes, has been used in various studies [24,25,[28][29][30]. Conversely, supervised classification is based on taking samples as training data and applying them to classify [31]. While the sub-pixel classification works to address the mixed pixel problem for coarse resolution satellite images, by treating the pixel as mixed and considering the areal proportion of each class [26,32], object-based image analysis is centered on generating image objects through image segmentation and performing the image classification on objects with different geometries rather than pixels [33][34][35].
It has been evident from the review that different types of satellite images and classification algorithms have been used and the choice is undoubtedly influenced by research objective and questions to be addressed by the work. Stephens and Diesing [36] tested the following six supervised classifiers: Classification Trees, Support Vector Machines, k-Nearest Neighbor, Neural Networks, Random Forest (RF), and Naive Bayes, and found good accuracy for tree-based methods and Naive Bayes. RF is a set of tree predictors; therefore, each tree depends on the value of an independently sampled random vector, and all the trees in the forest have the same distribution [37]. RF also performed well for the LC classification of a multi-source remote sensing and geographic data set in a study conducted by Gislason et al. [38]. Eisavi et al. [39] considered different scenarios for the LC mapping using RF, with the inclusion of land surface temperature data along with spectral data, and found the approach worked best for selected features over using just spectral signatures. Therefore, RF was used as a classification method for our work as well, in which we used spectral signatures along with derived indices and topography data: slope and elevation as input parameters.
Google Earth Engine (GEE), which is a cloud computing platform, has eased in extracting, analyzing, and performing different operations for the earth's observed data. It has been employed in different regions of the world to monitor the vegetation, water resources, and tracking the changes in land cover. As GEE is technologically robust and is an open resource containing historical images, it has become increasingly popular in LULC mapping. In the Asian region as well, it has been adopted in the countries such as India [40], Bangladesh [41], Singapore [42], China [43][44][45], Cambodia [46], Nepal [47], etc. GEE has also been applied in understanding the changes in land cover at the regional level such as in Central Asia [48] and the Tigris-Euphrates basin of the Middle East [49]. Automation in land cover mapping has also been possible due to the cloud computing feature of GEE. Xie et al. [50] proposed an automatic land cover classification method using time-series Landsat data using an International Geosphere-Biosphere Program (IGBP) classification scheme (Automatic Land-Cover Mapping).
In addition to understanding the present status of LC, it is equally important to track the motion of the center of change of the specific LC over the years. In the last years, with the growing urbanization, there have been significant changes in the types of regional land use. Therefore, the analysis of the central displacement of LC can help understand the development of the region and provide references for resource management and territorial planning [51]. For this purpose, a shift of the center of gravity (CoG) of the particular LC class has been studied [51][52][53]. The CoG models also help to explain regional shifts along the development axes [53].
Ecosystem services provide the basis for human life that are directly and indirectly contributing to human well-being and sustenance. The services range from filtering air and water by plants to decomposing wastes by microbes. These ecosystem services have been broadly grouped into the following four categories: provisioning, regulating, supporting, and cultural. With rapid development, Chinese researchers have been also using Ecosystem Service Value (ESV), an approach to assign monetary values to an ecosystem and its key ecosystem goods and services. Since Costanza et al. [54] published their article on the global value of ecosystem services, it has increasingly become more common and has been adapted in different areas modifying the coefficients of ESV [55][56][57][58]. A total of 17 types of ecosystem services were integrated into the original article by Costanza et al. [54], including climate regulation, water regulation, nutrient cycling, genetic resources, and culture. Moreover, the land is also classified into ecological land and non-ecological land, divided based on the positive or negative impacts on the ecosystem. Ecological land is the land type providing basic ecosystem services, includes non-built-up areas such as agriculture, water bodies, forest, grassland, and thus, has a high ESV. In contrast, nonecological land is an imperious, built-up area, thus weakening the ecosystem services [59]. With the population and urbanization pressure, the richest and diversified ecosystems are under a great threat. For instance, it was reported in 2018 that wetlands are declining three times faster than forests [60]. Thus, a study of ESV could help concerned stakeholders in understanding the rapid trend of urbanization, its severe environmental effects, and the load on the ecosystem.
Rupandehi, bordering its entire southern boundary with India, is a major corridor for the economic trade between Nepal and India. Being an economic hub, the district is facing the problem of managing urban growth. Home to happening towns, the district has undergone some remarkable changes and it has become extremely crucial to monitor the changes it has experienced or will experience in the future. While there have been some studies on the LC, they are generally at a national level [3,61]. The district is composed of a large percentage of agricultural land, followed by forest cover. Due to the result of unmanaged urbanization, the area has been severely affected by land fragmentation. This has resulted in improper drainage and flooding. Our study aims to provide LC and the pattern of a centroidal shift of the LC, which can help in curbing the effects, and present plans for tackling such situations. As the changes in the past and present conditions of the district are not known, it becomes extremely crucial to monitor the changes, as it can face severe drought and other climate crises in the future. Therefore, a detailed study on the spatiotemporal analysis of the LC is a must for a rapidly growing urban region such as Rupandehi.
The study focuses on classifying the LC of the Rupandehi District into five broad classes (agriculture, built-up, water bodies, forest, and barren) based on spectral information from Landsat Sensor using the RF algorithm. The major objectives are (1) to explore the spatiotemporal change of the Rupandehi District from 2005-2020, (2) to analyze the CoG trajectory of the land cover classes, (3) to calculate and analyze the ESV of individual LC classes and the whole district, and (4) to discuss the aspects and problems of the spatiotemporal analysis of developing regions. To the best of our knowledge, this is the first study that utilizes GEE, CoG, and ESV to understand the spatiotemporal analysis in the study area. Additionally, the workflow can be adopted in other places to quantity land changes and econometrics.

Study Area
Nepal, being one of the ten least urbanized countries in the world, is also one of the top ten fastest urbanizing countries [62]. Due to the concentration of educational and health facilities, and growing employment opportunities in Terai and major cities of valleys, there has been an upward trend of people migrating [63]. The Maoist Insurgency and the Madhesh Movement also had a significant impact on urbanization development [64]. This has led to uncontrolled population growth in the plains and valleys, resulting in unmanaged urbanization. Terai, being home to a large percentage of agricultural and forested land, is possessed by a great threat of haphazard development. Figure 1 shows the status of an urban area and population percentage in three different geographical regions (Himalayan, Hilly, and Terai). Without proper plans, the current trend would result in a sharp decline of arable land and deterioration of forests.  Figure 2 shows the study area, which is the Rupandehi District located in the Lumbini Province of Nepal. It is the third most populated district of Nepal. It covers an area of 1360 km 2 . Almost 60.2% of the area of the Rupandehi District is covered by agricultural land. The district lies in the southwestern part of Nepal. It shares a border with Nawalparasi District in the east, Kapilvastu District in the west, Palpa District in the north, and India in the south. The elevation of the district is between 100 and 1229 m from sea level. Of the area of Rupandehi District, 16.1% is in Churia Range and the rest is in the Terai region. It has major rivers such as Tinau and Rohini that are the source of water for its people. Most of the area of the Rupandehi District has a lower tropical climate (89.3%). It has two of the major economic hubs of Nepal: Bhairahawa and Butwal. Not only economically, but the Rupandehi District is also historically and culturally important. For the past many years, Rupandehi District has faced rapid LCC due to the migration of people from hilly to terai regions in search of flat cultivable land. It is very important to study the LCC pattern of the Rupandehi District. The spatiotemporal analysis of LCCs for the chosen study area, i.e., Rupandehi District is still not quantified. Using Landsat Imagery in the GEE platform, LCs were derived for different years. Figure 3 shows the workflow of this study. Annual cloud-free composites were obtained for the Landsat imageries for the years 2005, 2010, 2015, and 2020. Initially, unsupervised classification (clustering) was conducted to randomly create reference points that help in capturing the variability in pixels. In addition to the surface reflectance of the bands, different indices were also derived and used as input parameters for the RF. Topographic information such as elevation and slope derived from shuttle radar topography mission (SRTM) digital elevation model (DEM) was also given as input, as they also helped in accurately classifying imageries by obviating the effect of shadows from hills and buildings. To ensure the points have been assigned the right classes, it was overlayed with very high-resolution imagery from Google Earth Pro (GEP) and crosschecked with expert knowledge. The total points obtained were randomly divided into two sets of testing and validation. The classified maps were, then, assessed based on the validation set by generating a confusion matrix. Based on obtained LC areas, shifts in the CoG and ESV and spatiotemporal changes were analyzed.

Satellite Data
All the datasets used in this study are available in the Google Earth Engine (GEE) repository. The primary dataset used was the atmospherically corrected, tier 1 surface reflectance Landsat imagery. Landsat 5 was used for the years 2005 and 2010 while Landsat 8 was used for the years 2015 and 2020. Figure 4 shows the multispectral bands of Landsat 5 and 8 data. Operating at a spatial resolution of 30 m, both have a temporal resolution of 16 days. Revolving around in sun-synchronous, near-polar orbit, both data were acquired on the Worldwide Reference System-2 (WRS-2) path/row system. To acquire smooth continuous coverage of the scene without clouds, we used cloud masking and median filtering functions for the image in GEE, resulting in the annual median composite. In addition to spectral reflectance, different spectral indices (SIs) were derived and added as the input parameters for training. SIs have been used in vegetation studies, especially to assess the health of vegetation, which could provide valuable information in classification approaches. The following were the indices used for our study: a. Normalized Difference Vegetation Index (NDVI): An indicator of vegetation greenness, NDVI is the ratio of the spectral reflectance difference between the near-infrared (NIR) and red bands to the sum of the reflectances of the Near Infrared (NIR) and red bands [68]: c. Normalized Difference Built-up Index (NDBI): NDBI is the spectral index to detect built-up areas. It is calculated as normalized difference between Short Wave Infrared (SWIR) and green bands: Topographic information such as slope and elevation were also used in classification. They were derived from the SRTM, which is an international research effort in 2000 to obtain a DEM on a near-global scale. The data in GEE are SRTM V3 product (SRTM Plus) that is provided by NASA JPL at a resolution of 1 arc-second (approximately 30 m) and has undergone a void-filling process using open-source data (ASTER GDEM2, GMTED2010, and NED), as opposed to other versions that contain voids or have been void-filled with commercial sources.
GEE has the Landsat data dating back from its launch to the present day. GEE eases the way in extracting, analyzing, and performing different operations for the earth observed data. Additionally, GEP hosts the historical high-resolution data making it possible to test beforehand for each of the generated reference data for different years.

Reference Data
For generating reference data, unsupervised classification was used as an initial step. Generating a large amount of training data can be tedious. Therefore, unsupervised classification can help understand the distribution of samples and the clustering of similar features. It was applied by Wagle et al. [70] to generate their training data. Thus, we adopted ISO clustering to obtain 100 classes and applied stratified random sampling to generate 800 points. Each point was cross-checked by layering over the high-resolution imagery from GEP of the same year. To have improved representation and accuracy, additional knowledge from the author's prior knowledge of the study area (hometowns) was applied for the modification, addition, and deletion of points. Based on the ability of clusters and our ability to identify the LC, we generated the following 5 classes: agriculture, built-up, water bodies, forest, and barren (Table 1). These were the most dominant types of LC in the district. Due to the moderate resolution of Landsat, it was not advisable to further subdivide the classes, as it would have been difficult and would result in misclassification. For the images of 2010 (Landsat 5), a total of 806 reference points for all LC classes were cumulated, while for the images of 2015 (Landsat 8), the number of reference points was expanded to 819. Of the total points, 70% were randomly selected for training and the remaining 30% for validation.

Classification and Accuracy Assessment
Random Forest (RF) is a non-parametric ensemble learning algorithm combining many decision trees, each working independently of the other. It is an improvement on traditional decision trees, and it has been shown to yield better results in numerous case studies. Figure 5 shows the representation of an RF, with decision trees that consist of decision nodes, leaf nodes, and a root node. Each specific decision tree produces the individual output, which is the leaf node of the tree. The final output is chosen using a majority-voting procedure. In this situation, the final output of the system is the output chosen by the majority voting of decision trees. Decision trees are useful in sorting observations into categories based on features. They aid in differentiating between members of various groups and determining the degree of resemblance among members of the same group. The RF generates a class prediction for each tree, and the class with the highest votes becomes the model's prediction. RF chooses samples at random with a replacement, resulting in various trees. RF separates a random selection of features, resulting in more variation, as opposed to traditional decision trees, which analyze all potential features and choose the one that creates more variation. RF produces a robust result.

Decision Node
Leaf Node

Root Node
For accuracy assessment, producer's accuracy (PA), user's accuracy (UA), overall accuracy (OA), and Kappa statistics (Kappa) were calculated using the 30% validation data via confusion matrix. These statistics were used to evaluate as well as compare the classification results.

Centre of Gravity
The CoG is the point where the mass of the body seems to be concentrated. In the case of landcover type, it determines the spatial concentration of the landcover. Additionally, the CoG transfer trajectory determines the pattern of landcover change over time and helps in the study of LCC analysis. The concept of an LC CoG is derived from the idea of the population centroid and the economic centroid in geography.
CoG can be calculated as follows: where Xt and Yt are the X and Y coordinates of the CoG of the particular LC in year , Xi and Yi are the X and Y coordinates of the geometric center of the th patch, and Ati represents the area of the th patch.

Economic Value of Ecosystem Services
Due to the topographical and development stage similarities between China and Nepal, researchers have used a similar coefficient in Nepal for calculating the ESV values in major river basins of Nepal such as the Koshi river basin [71], Gandaki river basin [72], and Karnali river basin [73]. The studies had used the coefficients developed by Xie et al. [74], which had the value zero for built-up areas. Our study also used the same value for agriculture, water bodies, forest, and barren, but used the negative value of built-up areas developed by Duan et al. [75] to address the negative impacts due to growing urban areas, which was also used by Li et al. [59] in their study. The values used are shown in Table 2. The general formula for the value of ESV is as follows [76]: where ESV is the total ecosystem value (in USD), A is the area of land use type i (in hectares), and Pi represents the value coefficient for the ecosystem services provided by land type i (in USD/hectare).

Land Cover Classification and Accuracy Assessment
The LC maps of the Rupandehi District between the years 2005 and 2020 are shown in Figure 6. The district has seen some great variations and undergone some drastic changes. It can be seen that agricultural lands and forests are among the most abundant classes. An accuracy assessment was carried out on the resulting classified imagery using the PA, UA, OA, and Kappa from the confusion matrix to test the precision and accuracy of the classified imagery. As the Kappa values greater than 60% are considered as of substantial agreement, it was 62, 66, 68, and 69% for the years 2005, 2010, 2015, and 2020, respectively, showing the good agreement between the classified results and ground truth values in such fragmented LCs. The OA assesses the percentage of data correctly classified. It was found to be the highest (0.80) for 2020 and lowest (0.77) for 2005. For the years 2010 and 2015, it had the same value (0.79). Table 3 shows the accuracy assessment statistics obtained for different years.

Analysis of Centre of Gravity
After the preparation of LC maps of different years from 2005 to 2020, we calculated the centroid for each class and mapped the trajectory. The trajectory of the center is the centroid, and their lengths are shown in Figure 7 and Table 5.
Looking at the CoG transfer trajectory of agriculture, we can see that there has not been much change in the spatial pattern of the agricultural land compared to other classes over the past two decades. However, the centroid of agriculture has shifted toward the west, indicating that there was a great migration of people to the western part that started agricultural practices through various means. The transfer trajectory of barren is larger than other classes. There is a zigzag pattern in the trajectory change.

Ecosystem Value Services Analysis of Land Cover
Based on the value of ecosystem services, modified by Xie et al. [63] and Duan et al. [64] (Table 2), our study used and estimated the ESV between 2005, 2010, 2015, and 2020 in the Rupandehi District, as shown in Table 6.

Discussion and Conclusions
This study focused on deriving three major outcomes for the Rupandehi District. Based on the Landsat satellite image, we obtained LC information for the years 2005, 2010, 2015, and 2020 through RF classification in GEE, calculated the CoG for each LC type and analyzed the transfer along the years, and determined the change in ESV over the years. A cloud-based computing platform (GEE) facilitated acquiring and analyzing data for the study, including median composite and derivative indices from Landsat images and topographic information from SRTM DEM. GEE has certainly overcome the tiring effort to download each earth scene and offline analysis. This has reduced the time and cost, increasing efficiency. The study produced satisfactory results in classifying the LC with kappa statistics showing substantial agreement (>0.60) and accuracy in a good range (>0.70) for all the years. We also evaluated the ESV values for individual land cover classes and total values for each of the years, indicating the fluctuations. Lastly, at the end of the paper, we also discuss the aspects and problems of the spatiotemporal analysis of developing regions including the lack of sustainability in the developmental projects.
If we analyze the individual classes, agricultural land, which is the major percentage of LC, has been subjected to various fluctuations. Built-up areas, which were mostly centered in the two major cities, in the earlier years, have now become more spread out, inviting the problems such as land fragmentation and unmanaged urbanization. It increased from 18.14 sq. km. in 2005 to 71.15 km in 2020 (292% increment). A heavy increase in builtup areas can be seen progressing along the highway connecting major cities. Thus, the CoG of a built-up area is located around the centroid of the two cities. The decrease in water bodies from 2005 to 2015 can be explained by riverside encroachment and drying up of natural ponds due to extreme heat. However, artificial ponds for commercial fish farming have increased over the past few years, which can explain the increase in area in 2020 in a western direction. Forest area is subjected to less change in compared to other classes. It may also be the undisturbed class as the shift in the CoG is also very low. Strict rules against deforestation and the concept of community forest may be the reasons for the increasing forest area. In 2015, significant changes can be seen with decreasing forest and water bodies and increasing agriculture. The conversion of private green covers into farmland may be a contributing factor for that. Barren land increased from 2005 to 2010 and then, subsequently, decreased in the following years. This may be explained by the riverbank encroachment due to the human settlement. The barren lands in our study are generally the shores of rivers and ponds. The patches are often thin and elongated along with the water sources. Due to the limited spatial resolution of Landsat, the area of the barren land was smaller than the individual pixel most of the time, resulting in either water bodies or nearby major classes. Thus, it generated 0 values for the UA and the PA. The year 2005 had the highest ESV value due to the large area of ecological land (agriculture, water bodies, and forest). As the ESV values are directly proportional to the area of the land cover, increasing the built-up area increased negative ESVs. Due to the fluctuations in the area of other land cover classes, ESV values were also varied over the years.
We also compared our classification results of 2010 with a national LC map prepared by Integrated Mountain Development (ICIMOD) [61]. ICIMOD had divided the land into 12 classes and used object-based classification, and obtained a little higher overall accuracy compared to ours, of 85.13%, and a kappa of 0.82. Figure 8 shows the comparison of our result with that of ICIMOD, alongside high-resolution imagery from GEP for 2010. The comparison is performed by zooming in on Bhairahawa, the second-largest city in the district. While the ICIMOD landcover helps to identify the runway in the airport, our results just show certain signs of a built-up area. Both studies have identified agriculture to a quite good extent. However, our classification showed better results in identifying the water bodies and overall built-up areas, which is shown mostly as barren in the work by ICIMOD. Due to the inclusion of broad classes, ICIMOD, thus, had the advantage of identifying the peri-urban areas but had the limitation of misclassifying the major classes. For the year 2020, we compared our classified LC map with a global 10-meter resolution map of the Earth's land surface produced by ESRI, alongside high-resolution imagery from GEP for the same year. ESRI had produced the classification results with 10 classes processing over 400,000 Earth observations. Figure 9 shows the comparison between them focused at Butwal, the biggest city of the district. The classes were more generalized by ESRI, giving the smoothness, but failed to quantify the minor classes such as water bodies. It can be seen that water bodies were more accurately classified in our work compared to ESRI's LC. Similarly, the problem of the fragmentation of agricultural lands by haphazard built-up areas was observed more clearly in our work. Figure 9. Comparison of LC map around Butwal and its vicinity area for the year 2020 from top to bottom: 30-meter map obtained from our study, high-resolution imagery on Google Earth Pro (GEP), and global 10-meter map prepared by ESRI using Sentinel 2.
From the results of our study, we can discuss the following aspects of the spatiotemporal changes in developing regions: a. Sustainability is something every developmental effort should focus on. However, Nepal is often struggling to incorporate this into its strategies. District-level policies also lag in that aspect. Providing the historical and present status of LC can help in understanding how the land resources have been used over the years. The critical information obtained can be of high importance for planning more sustainably. It is highly recommended to incorporate ecological-based approaches to address disturbances in the ecosystem in the short-and long-term regulations. It must be ensured that the concepts of open spaces, green fuels, water conservation, forest restoration, food production, and watershed management are given high importance to curb the unplanned and haphazard development. b. Along with knowing the status of LC, it is extremely crucial to monitor where the particular LC class is focused. The analysis of the central displacement of LC can help understand the development of the region and provide references for resource management and territorial planning. The rapid shift of a particular LC in a direction can be an alarm to the overexploitation of another resource. Therefore, the analysis of LC along with the gravity shift of the classes is endorsed.
c. Calculating ESV is a definitive and appropriate approach for evaluating the ecosystem on a monetary basis providing the scientific ground for commanding the policies. In addition, ESV can be an effective way of communicating the results. It can provide insights into the load that we put on ecosystems and thus can help in implementing proper plans and strategies by district administration offices and local governmental bodies to stop the exploitation of resources. Integrating economic, ecological, and social dimensions in spatial planning ensures sustainable urban development plans. d. Accuracy is highly dependent on the quality of LC data and the unit value of ESV for each LC class type. Therefore, the improvement in LC data is highly recommended, along with more details and empirical studies for obtaining a unit value of ecosystem services. e. Studies such as ours are scientific efforts for developing regions, incorporating spatial econometrics and statistics with earth observation studies and GIS. Often, the topics are left out at the local level and are not considered in urban planning and development.
The study is the first of its kind in the Rupandehi District. The pattern and above discussion would be helpful to concerned authorities and stakeholders. We hope the study could be fruitful to the authorities interested in knowing the spatial changes of the land cover and also to those working on maintaining sustainable ecosystem services. Besides our effort to quantify the spatiotemporal analysis of the region, this study faced some challenges and limitations, which are as follows: i.
Stakeholders may be interested in the extraction of small classes such as road networks, types of forest, grassland, or pasture, which this study did not try to achieve due to the mid resolution of Landsat and the fragmented characteristics of the study area. Additionally, it can only be possible with a high level of spatial data. However, it may be difficult to acquire high-resolution images for long-term studies. ii.
Another problem is the class imbalance problem. While we had enough training, data were labeled as agriculture due to the large area covered by farmland, but the minority classes such as barren land and water bodies faced the problem of insufficient training data in comparison to the majority classes. However, we compensated for the problem by manually adding the points for the minority classes. iii.
Obtaining a proper ESV in the context of Nepal. We used the values derived from the Tibetan plateau, which shares a border with Nepal and has a similar economic and development stage. The values have also been used by previous researchers in the context of Nepal. However, there is certainly the need for detailed and thoroughly established ESV coefficients applicable to the whole of Nepal.
Future work can focus on resolving the limitations considering the new data, techniques, and coefficients. Either classifying more LC classes on Landsat or applying data fusion techniques with Sentinel imagery could be conducted to obtain finer and more precise LC classes. The class imbalance problem can be addressed using resampling techniques (either an under-sample majority class or an oversample minority class). Trying a variety of algorithms can be beneficial with imbalanced datasets. Besides RF, ensemble methods such as boosting, and bagging could be implemented to evaluate accuracy and robustness. For improvement in regional-based ESV factors, local and federal level collaborative research is recommended along with the inclusion of questionnaires for local ecologists and community-based surveys in the method.
Author Contributions: Conceptualization, Aman KC and Tri Dev Acharya; data curation, Aman KC and Nimisha Wagle; formal analysis, Aman KC and Tri Dev Acharya; investigation, Nimisha Wagle; methodology, Aman KC, Nimisha Wagle and Tri Dev Acharya; resources, Tri Dev Acharya; software, Nimisha Wagle; supervision, Tri Dev Acharya; validation, Nimisha Wagle; visualization, Tri Dev Acharya; writing-original draft, Aman KC; writing-review and editing, Nimisha Wagle and Tri Dev Acharya. All authors have read and agreed to the published version of the manuscript.