Utilizing Sentinel-2 Satellite Imagery for LULC and NDVI Change Dynamics for Gelephu, Bhutan

: Gelephu, located in the Himalayan region, has undergone significant development activities due to its suitable topography and geographic location. This has led to rapid urbanization in recent years. Assessing land use land cover (LULC) dynamics and Normalized Difference Vegetation Index (NDVI) can provide important information about urbanization trends and changes in vegetation health, respectively. The use of Geographic Information Systems (GIS) and Remote Sensing (RS) techniques based on various satellite products offers a unique opportunity to analyze these changes at a local scale. Exploring Bhutan’s mandate to maintain 60% forest cover and analyzing LULC transitions and vegetation changes using Sentinel-2 satellite imagery at 10 m resolution can provide important insights into potential future impacts. To examine these, we first performed LULC mapping for Gelephu for 2016 and 2023 using a Random Forest (RF) classifier and identified LULC changes. Second, the study assessed the dynamics of vegetation change within the study area by analysing the NDVI for the same period. Furthermore, the study also characterized the resulting LULC change for Gelephu Thromde, a sub-administrative municipal entity, as a result of the notable intensity of the infrastructure development activities. The current study used a framework to collect Sentinel-2 satellite data, which was then used for pre-and post-processing to create LULC and NDVI maps. The classification model achieved high accuracy, with an area under the curve (AUC) of up to 0.89. The corresponding LULC and NDVI statistics were analysed to determine the current status of the LULC and vegetation indices, respectively. The LULC change analysis reveals urban growth of 5.65% and 15.05% for Gelephu and Gelephu Thromde, respectively. The NDVI assessment shows significant deterioration in vegetation health with a 75.11% loss of healthy vegetation in Gelephu between 2016 and 2023. The results serve as a basis for strategy adaption required to examine the environmental protection and sustainable development management, and the policy interventions to minimize and balance the ecosystem, taking into account urban landscape.


Introduction
The evolution and transformation of land use land cover (LULC) has emerged as one of the drivers of global environmental change [1].The LULC changes are mainly driven by urbanization thrusts, population growth, depletion of natural resources and climate change impacts [2].Remote sensing (RS) technology is a powerful tool and offers a more convenient and efficient approach to obtaining information about the Earth's surface [1].The fusion of Geographic Information System (GIS) and RS has gained widespread prominence in the context of LULC mapping and monitoring, effectively supporting sustainable development programs in different geographic contexts [3].The focus on the Himalayan region has been a breeding ground for studies using GIS and RS techniques to examine LULC dynamics and provide valuable insights for sustainable planning and management efforts [4] and supported policy and decision-making processes in different regions, e.g., [5].Undoubtedly, the LULC maps are among the most important documents providing information for various applications in land use policy development, agricultural monitoring, urban planning, ecosystem services, nature conservation, and dynamic assessment [6].Furthermore, the application of artificial intelligence (AI) has been notable in the areas of earth monitoring and earth system sciences [7].Machine learning [8,9] and deep learning techniques [10][11][12] have emerged, amplifying the effectiveness of GIS and RS technologies.In addition, GISbased image classification methods have made significant progress in LULC studies [13,14], including the assessment and identification of changes in dynamics [15,16].In particular, image classification for geospatial mapping has received significant attention given the availability of medium to high-resolution multispectral satellite imagery [17], while the global products available are very broadly suited to local planning and management.
Sentinel-2 satellite imagery has been widely used for LULC studies [18][19][20][21].As of 2015, the Sentinel-2A Multispectral Imager (MSI) offers 13 optical bands: four bands at 10 m spatial resolution, six bands at 20 m spatial resolution, and three bands at 60 m spatial resolution and the 10 m spatial resolution can be used to improve the spatial resolution of other 20 or 60 m bands [22,23].Furthermore, the spatial and temporal resolution of Sentinel-2 imagery offers high potential for the generation of accurate datasets [24], including new opportunities for specific applications and investigations into lithological classification [25], mapping nitrogen-phosphorus ratio in feed [26], crop type mapping and identification [27,28], identification of ancient artifacts [29] among others.In particular, the combination of Sentinel-2 satellite' imagery and an appropriate classification method has been shown to provide better results, especially for the heterogeneous environments of the Himalayas [17].Similarly, Landsat datasets are also widely used for RS studies as they provide spatially resolved decadal time series over 30 years, despite offering 30 m resolution and a 15 m panchromatic band.Other sources of data for LULC studies are Satellite Pour l'Observation de la Terre (SPOT), Synthetic Aperture Radar (SAR), and Moderate Resolution Imaging Spectroradiometer (MODIS).
The SCP in QGIS provides a set of tools for geospatial analysis.The SCP is a Python tool for downloading and processing remote sensing data [30], which was developed in the QGIS environment [31].The SCP plugin works semi-automatically and provides a platform to download several RS products, facilitates the processing of satellite imagery, and automatically implements unsupervised and supervised classification by setting some input parameters.The SCP plugin implements the traditional pixel-based approach to image classification.Other methods also include object-based image analysis, which differs in the way pixels are grouped or how homogeneous pixels are converted into objects [22,32,33], and the accuracy of the results often depends on the geographical characteristics of the study area [34].The automated approaches are simple [35,36], and cost-effective without the need for physical interventions [37].Some of the research studies used the SCP for application in various land cover classifications, e.g., [38].In particular, in the heterogeneous Himalayan region of Bhutan, an SCP-based random forest classifier was first employed to great effect to classify satellite imagery for geohazard mapping [17].
In Bhutan, the absence of substantive RS-based studies exacerbates research and development gaps in Earth observation and monitoring.The limitations of such studies have led to knowledge gaps in the decision-making process and pose greater challenges in achieving Sustainable Development Goals (SDGs).To fill this gap, this study explores how RS techniques and image classification algorithms can be used efficiently, especially for LULC mapping, while maintaining the ecological benefits.Examination of existing literature reveals few studies on LULC, for example, one at the national level that presented the LULC map of Bhutan for 2010 [39] and 2016 [40,41] and the other for the LULC map of Thimphu district of 2018 [41].To capitalize on RS technology and GIS-based classification techniques, the proposed study examined the Gelephu region for which land cover information is lacking.Likewise, NDVI mapping and change detection are also performed to evaluate and monitor the health and density of vegetation in Gelephu.The NDVI index is widely used for various applications such as climate change response [42], vegetation degradation dynamics [43], and geo-environmental monitoring and land degradation assessment [44].
In summary, the purpose of the study is to examine the LULC and NDVI change at the local level.This study leverages Sentinel-2 satellite imagery and efficiently utilizes the GIS-based machine learning (RF) algorithm on SCP and provides the first local-scale LULC and NDVI maps for Gelephu.The expansion of urban infrastructure in Gelephu Thromde is also closely examined.The study highlights the change dynamics of LULC and NDVI that can be used for planning and managing sustainable development activities.

Study Area
Gelephu (Figure 1) is one of the local government administrations (gewogs) in Sarpang district and covers an area of 54.50 km 2 .According to the census 2017, the actual resident population is 4461 and consists of 1027 households [45].Geographically, Gelephu lies between, 25 to 27 • north and 90 to 91 • east with an elevation between 168 m in the southern floodplain and 1799 m towards the north.Gelephu consists of one of the administrative units of the municipality, Gelephu Thromde, and covers an area of 10.66 km 2 .The LULC dynamics in the study area have changed considerably due to recent urban planning and subsequent infrastructure development activities such as road network, building construction, automobile repair zones, and others, with most of the infrastructure development activities being concentrated in Gelephu Thromde.

Data Collection and Processing
To assess spatiotemporal LULC and NDVI and perform change detection analysis for Gelephu, the author conducted a study using satellite imagery data.The availability of Sentinel-2 satellite data between 2016 and 2023 was examined and downloaded from the USGS website https://ers.cr.usgs.gov/(accessed on 2 February 2022 and 3 June 2023, respectively).The Sentinel-2 images used in the current study for 2016 and 2023 correspond to 24 October 2016, and 21 January 2023, respectively.The selection of satellite images took into account a cloud coverage range of 0-5%, specifically targeting the drier months to minimize cloud interference.This approach was chosen to ensure data availability and to improve the quality of the image classification used for the analysis.The Sentinel-2 provides 13 spectral bands which consist of 10 m, 20 m, and 60 m spatial resolution products and the three red-edge bands that capture the strong reflectance of the vegetation in the near-infrared region of the electromagnetic spectrum.The Sentinel-2 multispectral imagery is used as it provides better accuracy for LULC as the 10 m RGB band sets and the B8 nearinfrared band provide much higher spatial resolution compared to other satellite products.Sentinel-2A images are classified as medium-resolution multispectral images.The most popular open-access remote sensing data for image classification are the satellite products Landsat and Sentinel-2.In the current study, Sentinel-2 imagery is used because it offers a relatively high spatial resolution of 10 to 20 m.A high spatial resolution was required to enable accurate image classification in the heterogeneous state of the study area.GIS is utilized in the current project to determine the coordinate reference system (CRS) of the study area using EPSG:32646-WGS 84/UTM-Zone 46N.The UTM Zone 46N includes parts of countries such as India, Bangladesh, Bhutan, Myanmar, and China.This CRS is based on the World Geodetic System 1984 (WGS 84) ellipsoid, a popular reference ellipsoid for worldwide positioning and mapping applications.Within the designated UTM zone, this coordinate reference system provides a localized and accurate basis for a range of areas for GIS operations such as mapping, spatial analysis, and data integration.
The methodological flow chart implemented for the image classification workflow is shown in Figure 2, which shows the data download, pre-processing, and post-processing of LULC and NDVI including the change detection and statistical evaluation.Two sets of Sentinel-2A satellite data from 2016 and 2023 are archived.Each data set is processed in the SCP in QGIS to first convert the Sentinel-2 band set raster from calibrated digital numbers (DN) to reflectance.Based on this degree of reflection, the colour composites required for image classification are then generated.The colour composites are essentially helpful in creating Regions of Interest (ROI) that helped formulate training data when implementing the RF algorithm during the image classification process [17].In particular, the colour composites such as band sets 4-3-2 (true colour), 5-6-2 (healthy vegetation), 5-4-3 (infrared colour, vegetation), and 7-6-4 (false colour, urban) were analysed and used as shown in Figure 3. GIS is utilized in the current project to determine the coordinate reference system (CRS) of the study area using EPSG:32646-WGS 84/UTM-Zone 46N.The UTM Zone 46N includes parts of countries such as India, Bangladesh, Bhutan, Myanmar, and China.This CRS is based on the World Geodetic System 1984 (WGS 84) ellipsoid, a popular reference ellipsoid for worldwide positioning and mapping applications.Within the designated UTM zone, this coordinate reference system provides a localized and accurate basis for a range of areas for GIS operations such as mapping, spatial analysis, and data integration.
The methodological flow chart implemented for the image classification workflow is shown in Figure 2, which shows the data download, pre-processing, and post-processing of LULC and NDVI including the change detection and statistical evaluation.Two sets of Sentinel-2A satellite data from 2016 and 2023 are archived.Each data set is processed in the SCP in QGIS to first convert the Sentinel-2 band set raster from calibrated digital numbers (DN) to reflectance.Based on this degree of reflection, the colour composites required for image classification are then generated.The colour composites are essentially helpful in creating Regions of Interest (ROI) that helped formulate training data when implementing the RF algorithm during the image classification process [17].In particular, the colour composites such as band sets 4-3-2 (true colour), 5-6-2 (healthy vegetation), 5-4-3 (infrared colour, vegetation), and 7-6-4 (false colour, urban) were analysed and used as shown in Figure 3.

The SCP and Random Forest Classifier
The SCP is a powerful QGIS plugin and offers a range of tools for geospatial analysis.The SCP is a Python graphical user interface (GUI) tool developed by [30] in the QGIS environment for downloading and processing remote sensing data.Typically, SCP currently supports Landsat, Sentinel, MODIS, and GOES satellite data and can process any size of corresponding RS data and essentially does not require a high-end computer.The SCP plugin works on most major operating systems as QGIS is a cross-platform opensource software compatible with Windows, macOS, and Linux.The SCP 7.10.11was used in the current study.The SCP plugin provides a semi-automatic workflow to download multiple RS products, facilitates the processing of satellite imagery, and implements unsupervised and supervised classification.The SCP plugin implements the traditional

The SCP and Random Forest Classifier
The SCP is a powerful QGIS plugin and offers a range of tools for geospatial analysis.The SCP is a Python graphical user interface (GUI) tool developed by [30] in the QGIS environment for downloading and processing remote sensing data.Typically, SCP currently supports Landsat, Sentinel, MODIS, and GOES satellite data and can process any size of corresponding RS data and essentially does not require a high-end computer.The SCP plugin works on most major operating systems as QGIS is a cross-platform open-source software compatible with Windows, macOS, and Linux.The SCP 7.10.11was used in the current study.The SCP plugin provides a semi-automatic workflow to download multiple RS products, facilitates the processing of satellite imagery, and implements unsupervised and supervised classification.The SCP plugin implements the traditional pixel-based approach to image classification.The SCP is also integrated with the Sentinel Application Platform (SNAP) application tools, enabling SCP to implement the Random Forest Classifier, a supervised machine learning algorithm.
Efficient classifier selection is required for successful image classification to classify the spectral signature features.RF classifier is widely applied as it offers efficient and high performance resulting in high accuracy of the classification model [46].In machine learning, the typical workflow involves dividing the dataset into two parts: the training set and the testing set.The training set is used to train the machine learning model, while the testing set is used to evaluate its performance and generalization ability on unseen data that executes an ensemble learning method based on a decision tree, combined with massive ensemble regression and classification trees [8].However, in the context of the QGIS-based SCP, the focus is on classification tasks for the ROI classes using RF.The vegetation, water, built-up, and barren land ROIs were created with polygons using composites grouped under each class with similar spectral signatures and tested for preview.These ROIs can be reviewed until the classification reaches sufficient accuracy.In addition, higher thresholds such as the number of training samples and decision trees to run the RF classifier are set to achieve higher performance of the classification model.
SCP uses the training data iteratively to improve the classification instead of using a separate testing set.The workflow is highlighted below.It is important to note that while SCP does not strictly follow the traditional training and test set approach, it does leverage custom training data to build a classification model and iteratively refine it to allow for accurate image classification and land cover mapping.
Colour composites (Figure 3) were utilized during the creation of training data to enhance the interpretation and discrimination of land classes.By combining different bands, colour composites provide a visually appealing representation of the land cover features, allowing for better identification and selection of ROIs.Essentially, the different colour composite images highlight unique spectral signatures that correspond to the different land cover classes, allowing for accurate delineation and representation of each class within the training dataset.In particular, the use of colour compositions facilitated the identification of representative spectral signatures of each class for training the supervised RF classifier.The RF uses the labelled training data resulting from the generated ROIs to accurately train a model to classify the associated pixels in the image under each predefined land class.

LULC Mapping
The current study focused on image classification for four primary LULC classes, namely water bodies, vegetation, built-up areas and barren land.The barren land typically includes sand and gravel river banks, landslide areas, and excavated construction areas.The RF classifier was used for two separate image classification schemes for the period between 2016 and 2023 and enabled comparative statistical analysis of LULC change dynamics occurring over seven years, including mapping of LULC change detection.LULC focuses on categorizing and classifying the types of land use activities and surface covers that indicate the spatial distribution of various land features.Statistical data was extracted from the image classification reports and a comparative LULC analysis was performed for Gelephu in 2016 and 2023.
Further to validate the LULC maps, the current study used a rapid validation method using the true colour composite from the 2023 Sentinel-2 image and Open Street Map (OSM) data.The building footprint of the built-up class was particularly taken into account during the validation.To validate the LULC map, two sets of building polygons were used.In the absence of OSM data in the southern zone of the study area, the building polygons from the true colour composite of the 2023 Sentinel-2 image have been meticulously drawn to accurately represent the locations and sizes of the buildings.For the

LULC Mapping
The current study focused on image classification for four primary LULC classes, namely water bodies, vegetation, built-up areas and barren land.The barren land typically includes sand and gravel river banks, landslide areas, and excavated construction areas.The RF classifier was used for two separate image classification schemes for the period between 2016 and 2023 and enabled comparative statistical analysis of LULC change dynamics occurring over seven years, including mapping of LULC change detection.LULC focuses on categorizing and classifying the types of land use activities and surface covers that indicate the spatial distribution of various land features.Statistical data was extracted from the image classification reports and a comparative LULC analysis was performed for Gelephu in 2016 and 2023.
Further to validate the LULC maps, the current study used a rapid validation method using the true colour composite from the 2023 Sentinel-2 image and Open Street Map (OSM) data.The building footprint of the built-up class was particularly taken into account during the validation.To validate the LULC map, two sets of building polygons were used.In the absence of OSM data in the southern zone of the study area, the building polygons from the true colour composite of the 2023 Sentinel-2 image have been meticulously drawn to accurately represent the locations and sizes of the buildings.For the northeast zone, OSM building footprint extraction was used for validation.The main objective of this validation process was to compare the identified building shapes with the corresponding pixels categorized as built-up areas on the Gelephu 2023 LULC map.The procedure involved overlaying the building shapes onto the LULC map allowing for a comparison between these two sets of data.

NDVI Mapping
To further delve into the vegetation monitoring tools, vegetation change assessment was carried out separately using NDVI.NDVI serves a specific purpose of monitoring vegetation density and health status.
Gelephu, a scenic town in the southern part of Bhutan, is known for its harmonious natural landscapes and rich vegetation.The vegetation change dynamics is a critical component for effective environmental monitoring and conservation for sustainable land management.In this study, the author assessed the NDVI to identify vegetation change dynamics in Gelephu using the same Sentinel-2 satellite images for 2016 and 2023, taking into account the consistency of spatiotemporal and spectral resolution.The NDVI, derived from the remote sensing satellite imagery serves as a valuable indicator of vegetation health and density.Using the QGIS software 2.28.6, the NDVI is calculated as the ratio of the difference between the near-infrared (NIR) and red spectral bands normalized by their sum.This index efficiently captures vegetation density and vitality, providing insights into changes in plant health and growth.By applying the NDVI formula to band 8 (NIR) and band 4 (red), we obtained NDVI values for each pixel defined in Equation ( 1), which were then mapped to visualize the spatial distribution of vegetation in Gelephu.
To assess the changes between the 2016 and 2023 NDVI maps, a reclassification process was performed to categorize the NDVI values into different vegetation classes.The reclassification aimed to provide a more meaningful representation of the vegetation density and the dynamics of change.

Accuracy Assessment Methods
SCP is designed to train the classification model (RF) for machine learning by iteratively improving the input training data set consisting of ROIs (polygons) for each class.The classification can be revised by refining the training area without test data.However, to verify the level of accuracy achieved, the current study used additional testing data consisting of 51 random points created with Google Earth.These points were used to perform the accuracy assessment through an Area Under the Curve (AUC) approach for the RF classification model for both LULC maps.
In the current study, the sample points used are typically well-known ground truth labels that illustrate the presence or absence of a particular land cover class.The ArcSDM (Spatial Data Modeller) tool in ArcGIS was employed to construct the receiver operating characteristic curve (ROC) curves.A ROC curve is created by plotting the true positive rate (TPR) against the false positive rate (FPR) for various classification thresholds (spectral signatures of respective classes).The TPR and FPR are computed using the equations Equation (2) and Equation (3), respectively.The TPR represents the proportion of positive instances correctly classified, while the FPR indicates the proportion of negative instances incorrectly classified as positive.By systematically adjusting the classification threshold, we can observe how the model's TPR and FPR change, resulting in a curve that summarizes its discriminatory power across different operating points.The ROC curve provides a visual representation of the model's trade-off between sensitivity (recall) and specificity.A perfect model would achieve a TPR of 1 and an FPR of 0, resulting in an ROC curve that hugs the top-left corner of the graph.The closer the curve comes to this ideal point, the better the model's performance.The diagonal line from the bottom left to the top right represents the random classifier, indicating an equal chance of correctly classifying positive and negative instances.Any model that falls below this line is considered worse than random.
The True Positive Rate (TPR), also known as Sensitivity or Recall, is defined as: Similarly, the False Positive Rate (FPR) is defined as: The area under the ROC curve, commonly referred to as the Area Under the Curve (AUC), quantifies the overall performance of a classification model.A perfect model has an AUC of 1. Higher AUC values indicate better discriminative ability.AUC serves as a reliable metric for comparing different models, as it remains unaffected by class imbalance.

LULC Maps
Figure 4 shows the spatial distribution and visual representation of the LULC changes that have occurred over seven years.In Figure 4a, the 2016 LULC map predominantly shows vegetation cover and also identifies pockets of urban development, indicating the early phases of urbanization.Conversely, the LULC map for 2023 (Figure 4b) indicates an obvious change in the landscape of Gelephu.The growth of urban areas has been substantial, with an increase in buildings, roads, and other infrastructure, indicating a rapid pace of urbanization and migration that has taken place in Gelephu in recent years.A prominent shift from green space to urbanized areas is notable as demand for residential and commercial space soared.In addition, the implementation of local area plans (LAP) and policies were reflected in the LULC map for 2023.Gelephu LULC maps from 2016 and 2023 showed major land-use changes, concentrated in the southern zones, which comprise the Gelephu Thromde sub-administrative unit.The implementation of the LAP falls into this area.In order to specifically assess the urban extent of the Gelephu Thromde, the study additionally focused on this administrative boundary.Gelephu Thromde LULC maps for 2016 and 2023 are shown in Figure 5.

LULC Change Detection
Additionally, change detection on the LULC was performed (see Figure 6) using the postprocessing tools in SCP in QGIS.The main goal of change detection mapping is to identify the change or transition of each land-use class to each other in the study area.Gelephu's LULC change detection map is shown in Figure 6a.The colors of the legend in Figure 6 indicate quantification of the current state of the land classes in 2023 (Figure 6a) that have already transitioned within the assessment period of the 2016 study, with some of the classes with no changes.Further, the LULC change detection map is rendered in Figure 6b to visually depict the delineation of land-use class change spots and represent the actual land-use transition through the year 2023.In Figure 6b, a white background is assigned to the land-use class to represent the no-change land-use class.

LULC Change Detection
Additionally, change detection on the LULC was performed (see Figure 6  The LULC change detection results from Gelephu show several significant transitions (Figure 7).The largest land cover class, vegetation, was fairly stable, with an area of 34.95 km 2 .This provides information about the notable change in the vegetation areas over just a seven-year study period.The changes between water and vegetation are relatively small.The water-to-vegetation transition covers an area of 0.04 km 2 , indicating slight dispersal of vegetation into water bodies.Similarly, transitions involving water and built-up areas are relatively small in scale.Conversely, the transition from vegetation to water covers an area of 0.05 km 2 , indicating a slight decrease in vegetation cover and an increase in water bodies.A total area of 2.62 km 2 shows a transition from barren land to vegetation.The conversion of vegetation into built-up areas covers an area of 4.82 km 2 , indicating urban expansion and reflecting the conversion of natural vegetation into builtup infrastructure, such as residential, commercial, or industrial facilities.The conversion of vegetation to barren land covers an area of 4.48 km 2 , suggesting possible land degradation or land-use changes.This transition between barren land and built-up areas illustrates urban expansion into previously unproductive or degraded areas.The conversion of barren land to built-up areas covers 1.16 km 2 , indicating land utilization for urban development purposes.Conversely, the conversion of developed land to barren land covers an area of 1.34 km 2 .The barren land cover class remained relatively stable, with an area of 3.38 km 2 , indicating the persistence of unproductive or degraded areas within the study area.The LULC change detection results from Gelephu show several significant transitions (Figure 7).The largest land cover class, vegetation, was fairly stable, with an area of 34.95 km 2 .This provides information about the notable change in the vegetation areas over just a seven-year study period.The changes between water and vegetation are relatively small.The water-to-vegetation transition covers an area of 0.04 km 2 , indicating slight dispersal of vegetation into water bodies.Similarly, transitions involving water and built-up areas are relatively small in scale.Conversely, the transition from vegetation to water covers an area of 0.05 km 2 , indicating a slight decrease in vegetation cover and an increase in water bodies.A total area of 2.62 km 2 shows a transition from barren land to vegetation.The conversion of vegetation into built-up areas covers an area of 4.82 km 2 , indicating urban expansion and reflecting the conversion of natural vegetation into built-up infrastructure, such as residential, commercial, or industrial facilities.The conversion of vegetation to barren land covers an area of 4.48 km 2 , suggesting possible land degradation or land-use changes.This transition between barren land and built-up areas illustrates urban expansion into previously unproductive or degraded areas.The conversion of barren land to built-up areas covers 1.16 km 2 , indicating land utilization for urban development purposes.Conversely, the conversion of developed land to barren land covers an area of 1.34 km 2 .The barren land cover class remained relatively stable, with an area of 3.38 km 2 , indicating the persistence of unproductive or degraded areas within the study area.

LULC Statistics
The LULC statistics for Gelephu are presented in Table 1.In 2016 the area covered by vegetation was 44.31 km 2 making up 81.30% of the total area.However, by 2023 the vegetation area had shrunk to 37.62 km 2 , which represents a decrease of, about 11.99%.This means that the percentage of land covered by vegetation in relation to the area dropped to around 69.31%.As for the water area in 2016, it accounted for 0.52 km 2 or about 0.95% of the total area due to being occupied by a river.However, in 2023 this water area decreased significantly to about 0.02 km 2 indicating a decline of roughly 0.91%.This change can be attributed to alterations in the watercourse itself.Consequently, the proportion of water coverage in the area decreased to approximately just, 0.04%.The built-up area in 2016 was 4.03 km 2 , representing 7.39% of the total area.By 2023, the built-up area increased to 7.08 km 2 , showing a significant increase of 5.65%.The portion of the built-up area in the total area increased to 13.04%.Likewise, the barren land coverage was 5.64 km 2 in 2016, accounting for 10.35% of the total area.The barren land area increased to 9.55 km 2 in 2023, representing a significant increase of 7.25%.The proportion of barren land in the total area increased to 17.60%.Overall, between 2016 and 2023, vegetation and water cover decreased while built-up areas and barren land increased significantly.
Additionally, Table 1 lists the LULC statistics for Gelephu Thromde.The LULC statistics between 2016 and 2023 show significant transformations in Thromde's land use dynamics.Vegetation, which accounted for 67.04% of the total area in 2016, saw a significant decrease of 29.93%.The vegetation area decreased from 7.15 to 3.96 km 2 , with a corresponding decrease in percentage cover from 67.04% to 37.11%.While water is a smaller land cover category, it saw a slight increase in both absolute area and percentage coverage.The water area increased from 0.03 to 0.04 km 2 with a corresponding increase in percentage cover from 0.25% to 0.37%.On the other hand, both the built-up area and the undeveloped land showed significant growth during the period under study.The built-up area increased by 15.05%, from 1.73 to 3.34 km 2 , with the percent coverage rising from 16.24% to 31.28%, indicating rapid urbanization.Similarly, the area of barren land increased by 14.76%, from 1.76 to 3.33 km 2 , with the percent coverage increasing from 16.47% to 31.23%.

LULC Statistics
The LULC statistics for Gelephu are presented in Table 1.In 2016 the area covered by vegetation was 44.31 km 2 making up 81.30% of the total area.However, by 2023 the vegetation area had shrunk to 37.62 km 2 , which represents a decrease of, about 11.99%.This means that the percentage of land covered by vegetation in relation to the area dropped to around 69.31%.As for the water area in 2016, it accounted for 0.52 km 2 or about 0.95% of the total area due to being occupied by a river.However, in 2023 this water area decreased significantly to about 0.02 km 2 indicating a decline of roughly 0.91%.This change can be attributed to alterations in the watercourse itself.Consequently, the proportion of water coverage in the area decreased to approximately just, 0.04%.The built-up area in 2016 was 4.03 km 2 , representing 7.39% of the total area.By 2023, the built-up area increased to 7.08 km 2 , showing a significant increase of 5.65%.The portion of the built-up area in the total area increased to 13.04%.Likewise, the barren land coverage was 5.64 km 2 in 2016, accounting for 10.35% of the total area.The barren land area increased to 9.55 km 2 in 2023, representing a significant increase of 7.25%.The proportion of barren land in the total area increased to 17.60%.Overall, between 2016 and 2023, vegetation and water cover decreased while built-up areas and barren land increased significantly.Additionally, Table 1 lists the LULC statistics for Gelephu Thromde.The LULC statistics between 2016 and 2023 show significant transformations in Thromde's land use dynamics.Vegetation, which accounted for 67.04% of the total area in 2016, saw a significant decrease of 29.93%.The vegetation area decreased from 7.15 to 3.96 km 2 , with a corresponding decrease in percentage cover from 67.04% to 37.11%.While water is a smaller land cover category, it saw a slight increase in both absolute area and percentage coverage.The water area increased from 0.03 to 0.04 km 2 with a corresponding increase in percentage cover from 0.25% to 0.37%.On the other hand, both the built-up area and the undeveloped land showed significant growth during the period under study.The built-up area increased by 15.05%, from 1.73 to 3.34 km 2 , with the percent coverage rising from 16.24% to 31.28%, indicating rapid urbanization.Similarly, the area of barren land increased by 14.76%, from 1.76 to 3.33 km 2 , with the percent coverage increasing from 16.47% to 31.23%.

Accuracy Assessment of LULC Maps
The LULC raster files are imported into ArcGIS and 51 testing sample points (Figure 8) are overlayed.A TP occurs when a test sample representing a particular land cover class is correctly classified as that class in the LULC map.Conversely, FP results arise when a test sample representing a particular land cover class is incorrectly classified as another class in the LULC map.TP indicates the accuracy of the classification model (RF) in identifying a particular land cover class.The higher the number of TPs, the more accurate the classification is for that class.

Accuracy Assessment of LULC Maps
The LULC raster files are imported into ArcGIS and 51 testing sample points (Figure 8) are overlayed.A TP occurs when a test sample representing a particular land cover class is correctly classified as that class in the LULC map.Conversely, FP results arise when a test sample representing a particular land cover class is incorrectly classified as another class in the LULC map.TP indicates the accuracy of the classification model (RF) in identifying a particular land cover class.The higher the number of TPs, the more accurate the classification is for that class.The ROC for image classification models is shown in Figure 9.The RF classification model resulted in 0.883 AUC for 2016 LULC mapping (Figure 9a) and 0.890 AUC for 2023 (Figure 9b) for Gelephu.The similarity in AUC values between the two time periods (2016 and 2023) indicates the robustness and consistency of the RF model's performance.The model's ability to maintain a high level of accuracy across different periods suggests that it is reliable and can effectively handle temporal changes in the landscape.The ROC for image classification models is shown in Figure 9.The RF classification model resulted in 0.883 AUC for 2016 LULC mapping (Figure 9a) and 0.890 AUC for 2023 (Figure 9b

NDVI Change Detection and Statistics
Analysis of the data for NDVI change detection between 2016 and 2023 (Table 3) shows significant shifts in land cover types within the study area (Figure 11).There was a significant decrease in water bodies: the area decreased from 0.96 to 0.06 km 2 and the percentage of water cover fell from 1.76% to 0.11%.This 1.65% decrease in water bodies was mainly observed due to a change in the watercourse outside the study area zone.The category of "No vegetation" showed the most significant increase in both absolute area and percent coverage.The area of no vegetation extended from 2.87 to 8.08 km 2 , with the percent coverage increasing from 5.27% to 14.83%.This significant growth of 9.56% suggests a loss of vegetation cover and raises concerns about potential ecological impacts, such as soil erosion and reduced biodiversity.The class of "Low vegetation" also accords with substantial growth, with the area increasing from 2.07 to 24.54 km 2 and the percent coverage rising from 3.80% to 45.04%.This notable expansion of 41.24% indicates vegeta-

NDVI Change Detection and Statistics
Analysis of the data for NDVI change detection between 2016 and 2023 (Table 3) shows significant shifts in land cover types within the study area (Figure 11).There was a significant decrease in water bodies: the area decreased from 0.96 to 0.06 km 2 and the percentage of water cover fell from 1.76% to 0.11%.This 1.65% decrease in water bodies was mainly observed due to a change in the watercourse outside the study area zone.The category of "No vegetation" showed the most significant increase in both absolute area and percent coverage.The area of no vegetation extended from 2.87 to 8.08 km 2 , with the percent coverage increasing from 5.27% to 14.83%.This significant growth of 9.56% suggests a loss of vegetation cover and raises concerns about potential ecological impacts, such as soil erosion and reduced biodiversity.The class of "Low vegetation" also accords with substantial growth, with the area increasing from 2.07 to 24.54 km 2 and the percent coverage rising from 3.80% to 45.04%.This notable expansion of 41.24% indicates vegetation changes, probably due to natural or human-induced factors such as afforestation or agricultural activities.The class of "Healthy vegetation" exhibited moderate growth, with the area increasing from 7.13 to 21.27 km 2 , and the percent coverage expanding from 13.09% to 39.04%.This positive trend of 25.95% suggests an increase in vegetation density and health within the study area.Surprisingly, the class of "Very healthy vegetation" experienced a significant decline.The area of very healthy vegetation decreased dramatically from 41.46 to 0.53 km 2 , and the percent coverage declined from 76.09% to 0.98%.This significant 75.11% decline in very healthy vegetation raises questions about possible environmental degradation, such as deforestation or land-use changes.Gelephu has seen significant LULC changes in recent years.According to this study, within a short period of seven years, urban settlement area change increased by 5.65% (13.03% in total) and vegetation cover decreased by 11.99% due to land use practices, a situation that is in line with the main trend of LULC, which is urbanization in developing countries [47].The results underline the complex interplay between urban development and environmental changes in the region.In this regard, Aryal et al. [48] indicated that the results could vary depending on the geographical characteristics of the study area.Although several factors influence urbanization, one of those detected in our research is migration, which played a crucial role in determining the Gelephu's LULC dynamics.With the growing population and increasing economic activity in Gelephu, there has been a remarkable transformation of agricultural and rural areas into urbanized ones.A related context has been identified, indicating that population growth and economic activities play a significant role in influencing urban expansion and the arrangement of built-up areas [49,50].One of the biggest LULC changes in Gelephu is the conversion of wet farmland to dry land for building purposes.Due to increasing demand for residential and commercial infrastructure, agricultural land has been repurposed for the development of buildings, roads, and other urban facilities.These multiple factors have led to the changes, impacted traditional farming practices, and contributed significantly to the decline in farmland.As expected, the key changes in the LULC dynamics come from the implementation of Local Area Plans (LAP) and policies.The construction of road networks has further augmented the impact on the dynamics of LULC in Gelephu.As the city is one of the main gateways to Bhutan and borders Assam in India, transport infrastructure has been upgraded to facilitate connectivity and trade.This has resulted in significant changes in land cover structure, including the clearing of vegetation and remodelling of surrounding areas to accommodate transport networks, in addition to the construction of buildings and other service infrastructures.An in-depth analysis of land cover changes within Gelephu Thromde over the same period confirms the overarching trends observed in the wider Gelephu area.The significant decrease in vegetation cover by 29.93%, from 67.04% to 37.11%, raises concerns about environmental sustainability.At the same time, the significant growth of the built-up area by 15.05% and barren land by 14.76% highlights the impact of rapid urbanization on the local landscape.Overall, these results suggest a transformative shift in land use patterns that may pose challenges for future environmental protection and land management strategies.
Other factors that cause the land cover changes in Gelephu are due to natural factors such as flooding and landslides.Due to the topography and monsoon climate, the region is prone to natural hazards leading to river course changes, erosion, and sediment deposition.These events have contributed significantly to the loss of vegetation and changes in water bodies.It is also important to note that the LULC changes observed in the maps are not solely due to human activity.The region's vulnerability to natural hazards such as floods and landslides was reflected in changes LULC map in certain areas resulting in vegetation loss and changes in the watercourse.To this end, the LULC maps for Gelephu provide important evidence and valuable insight into the dynamic nature of LULC and illustrate the significant transitions that have occurred within a relatively short time frame.In addition, the use of finer RS products allows differentiation of more subtle geometric differences in land cover types than coarse or medium spatial resolution images [51], thereby leveraging the use of Sentinel-2 satellite imagery for current studies.

NDVI Change Dynamics
Analysis of NDVI change detection data reveals complicated shifts in vegetation cover, enabling a more nuanced understanding of observed land cover changes.The notable 9.56% increase in the category of "No vegetation" underscores areas facing degradation, while the expansive 41.24% growth in "Low vegetation" signals potential ecological disturbances.
The staggering 75.11% decrease in "Very healthy vegetation" emphasizes the severity of vegetation ecology.These NDVI results also complement the land cover observations and provide additional insights into the varying degrees of vegetation health and environmental dynamics within the study area.The NDVI indices also represent one of the approaches to monitoring urban built-up areas, which are generally used to optimize land use patterns [52].
The observed dynamics of NDVI change from 2016 to 2023 have implications for the impacts of climate change in the study area.These NDVI change dynamics are consistent with broader discussions about the influence of climate change on vegetation patterns.The findings serve as a call for continuous monitoring and adaptation strategies to manage the evolving environmental landscape and mitigate the potential consequences of climate change impacts in the region.This particularly highlights the need to further investigate the impacts on land surface temperature, which is directly related to NDVI change dynamics [53,54], through the utilization of other suitable satellite products.In summary, the NDVI change statistics essentially provide useful information and quantification on vegetation health and density.The statistics are valuable for monitoring changes in plant growth, stress, and environmental conditions, and provide insights into the health of ecosystems, potentially helping to conserve biodiversity [55] in the context of the impacts of climate change.

Validation
Analysis revealed that a significant majority of the building shapes perfectly aligned with the pixels classified as built-up areas on the LULC map (Figure 12).This agreement between both datasets provided evidence regarding the accuracy and dependability of the LULC map generated by utilizing the RF classifier.Validating the LULC map through comparison, with building shapes demonstrates the effectiveness of using this classification method.It provides confidence in the ability of the random forest classifier to accurately identify and classify land cover types, particularly built-up areas.This validation process contributes to the overall reliability of the LULC map and enhances its suitability for various applications, such as urban planning, resource management, and environmental monitoring in the Gelephu region.the NDVI change statistics essentially provide useful information and quantification on vegetation health and density.The statistics are valuable for monitoring changes in plant growth, stress, and environmental conditions, and provide insights into the health of ecosystems, potentially helping to conserve biodiversity [55] in the context of the impacts of climate change.

Validation
Analysis revealed that a significant majority of the building shapes perfectly aligned with the pixels classified as built-up areas on the LULC map (Figure 12).This agreement between both datasets provided evidence regarding the accuracy and dependability of the LULC map generated by utilizing the RF classifier.Validating the LULC map through comparison, with building shapes demonstrates the effectiveness of using this classification method.It provides confidence in the ability of the random forest classifier to accurately identify and classify land cover types, particularly built-up areas.This validation process contributes to the overall reliability of the LULC map and enhances its suitability for various applications, such as urban planning, resource management, and environmental monitoring in the Gelephu region.

Conclusions
The study performed spatiotemporal LULC and NDVI mapping for Gelephu.In addition, the LULC maps for Gelephu Thromde have also been closely evaluated.Both the LULC and NDVI maps provide valuable information for decision-making processes related to land management, urban planning, and environmental conservation.
The results of the study highlight significant changes in Gelephu's land cover types between 2016 and 2023.The vegetation cover witnessed a decrease of 11.99%, resulting in

Figure 1 .
Figure 1.Study area.(a) Geographical location of Bhutan; (b) Digital elevation model of the study area.

Figure 1 .
Figure 1.Study area.(a) Geographical location of Bhutan; (b) Digital elevation model of the study area.

Figure 2 .
Figure 2. Implementation of the methodology.Bands B2 to B8A and B11-B12 were used for image classification (LULC) in SCP (a total of 10 optical bands with 10 m and 20 m resolution are included in the image processing.B1, B11, and B12 with a resolution of 60 m are excluded).For NDVI, only the two optical bands B4 and B8 are used for analysis.

Figure 2 .
Figure 2. Implementation of the methodology.Bands B2 to B8A and B11-B12 were used for image classification (LULC) in SCP (a total of 10 optical bands with 10 m and 20 m resolution are included in the image processing.B1, B11, and B12 with a resolution of 60 m are excluded).For NDVI, only the two optical bands B4 and B8 are used for analysis.

•
User-defined training areas: The user manually selects representative samples from the image, called ROI or training polygons.These areas should represent the different land cover classes present in the image.• Feature extraction: SCP extracts a set of features (attributes) from the image data within the training areas.These features could include spectral values from different bands, texture measures, indices, or other relevant information.• RF training: SCP utilizes the ROIs and the associated extracted features to train an RF classifier.RF is an ensemble machine learning algorithm that constructs a multitude of decision trees and combines their predictions for classification.• Classification: Once the RF model is trained, SCP applies it to the entire image to classify the remaining pixels.The model uses the extracted features from the image to assign each pixel to a specific land cover class.• Iterative improvement: SCP assesses the initial classification results visually (preview).If the results are unsatisfactory, the user can refine the training areas, extract additional features, or adjust the classification parameters to improve the classification accuracy.
) using the postprocessing tools in SCP in QGIS.The main goal of change detection mapping is to identify the change or transition of each land-use class to each other in the study area.Gelephu's LULC change detection map is shown in Figure 6a.The colors of the legend in Figure 6 indicate quantification of the current state of the land classes in 2023 (Figure 6a) that have already transitioned within the assessment period of the 2016 study, with some of the classes with no changes.Further, the LULC change detection map is rendered in Figure 6b to visually depict the delineation of land-use class change spots and represent the actual land-use transition through the year 2023.In Figure 6b, a white background is assigned to the land-use class to represent the no-change land-use class.

Figure 6 .
Figure 6.LULC change detection map of Gelephu.(a) Gelephu LULC map of 2023 quantifying the transition area of each land class to one another, including those classes with no changes.The quantifications are represented by the legend colors (Green: Vegetation; Blue: Water; Red: Built-up; Yellow: Barren land), (b) Change detection map of Gelephu delineated with 'No change classes' on a white background.This map represents transitional classes only.

Figure 6 .
Figure 6.LULC change detection map of Gelephu.(a) Gelephu LULC map of 2023 quantifying the transition area of each land class to one another, including those classes with no changes.The quantifications are represented by the legend colors (Green: Vegetation; Blue: Water; Red: Built-up; Yellow: Barren land), (b) Change detection map of Gelephu delineated with 'No change classes' on a white background.This map represents transitional classes only.

Figure 7 .
Figure 7. LULC transition pattern for Gelephu.The colors represent the land class color scheme used for LULC mapping (Green: Vegetation; Blue: Water; Red: Built-up; Yellow: Barren land).
land Barren land to Built-up Built-up to Barren land Area (x km 2 ) Vegetation (no change = 34.95Water (no change) = 0.06 Built-up (no change) = 1.04 Barren land (no change) = 3.38

Figure 7 .
Figure 7. LULC transition pattern for Gelephu.The colors represent the land class color scheme used for LULC mapping (Green: Vegetation; Blue: Water; Red: Built-up; Yellow: Barren land).

Figure 8 .
Figure 8. Testing data set.The confusion matrices (Table 2) of ROC analysis for the 2016 and 2023 image classification model by RF classifier are shown in Table 2.The actual classes represent the true positive (TP) and the predicted classes represent the false positive (FP).The point where TP and FP intersect in the confusion matrices provides pivotal information on the accuracy of the classification model.The classification model for 2016 and 2023 indicates a consistent TP result across both years, indicating stability in positive class predictions.The vegetation and barren land classes show a higher positive class prediction among the other two classes, least of all for the water classes.

Figure 8 .
Figure 8. Testing data set.The confusion matrices (Table 2) of ROC analysis for the 2016 and 2023 image classification model by RF classifier are shown in Table 2.The actual classes represent the true positive (TP) and the predicted classes represent the false positive (FP).The point where TP and FP intersect in the confusion matrices provides pivotal information on the accuracy of the classification model.The classification model for 2016 and 2023 indicates a consistent TP result across both years, indicating stability in positive class predictions.The vegetation and barren land classes show a higher positive class prediction among the other two classes, least of all for the water classes.

Figure 10
Figure 10 shows the NDVI maps of Gelephu.The NDVI is primarily a measure of the greenness and health of vegetation and is the most commonly used vegetation index [45].The NDVI value range is between −1 to +1.The positive index value indicates a vegetated area, and the negative value suggests a non-vegetated area.In Figure 10, the 2016 NDVI map shows an NDVI range between −0.35 to 0.86 and −0.07 to 0.71 in 2023.The shift in the

Figure 10
Figure 10 shows the NDVI maps of Gelephu.The NDVI is primarily a measure of the greenness and health of vegetation and is the most commonly used vegetation index [45].The NDVI value range is between −1 to +1.The positive index value indicates a vegetated area, and the negative value suggests a non-vegetated area.In Figure 10, the 2016 NDVI map shows an NDVI range between −0.35 to 0.86 and −0.07 to 0.71 in 2023.The shift in the NDVI paradigm from −0.35 to −0.07 suggests an increase in the non-vegetation area and thus a decline in healthy vegetation patterns with a change in index value from 0.86 to 0.71.The reclassified NDVI maps (Figure 11) consist of five categorical vegetation classes including water bodies identified with NDVI values less than zero that have no vegetation.Other classes include no vegetation, low vegetation, high vegetation, and very high vegetation, indicating denser and healthy vegetation, respectively.Areas with NDVI values between 0 and 0.2 were classified as bare.NDVI values between 0.2 and 0.4 were labelled as low vegetation, indicating areas with limited vegetation cover.High vegetation was defined for NDVI values between 0.4 and 0.6, representing areas with moderate vegetation density.Finally, NDVI values greater than 0.6 were categorized as very high vegetation, indicating areas of dense and lush vegetation cover.This reclassification process allowed clear differentiation of vegetation classes and enabled a comprehensive detection of the changes that occurred between the 2016 and 2023 NDVI maps.By comparing the distribution and extent of each vegetation class, patterns of vegetation change or shifts in vegetation density are identified and analysed.Appl.Sci.2024, 14, x FOR PEER REVIEW 17 of 23

Table 1 .
Summary of LULC statistics for Gelephu and Gelephu Thromde respectively.

Table 1 .
Summary of LULC statistics for Gelephu and Gelephu Thromde respectively.

Table 2 .
Confusion matrices for ROC analysis for the 2016 and 2023 image classification models, respectively.

Table 2 .
Confusion matrices for ROC analysis for the 2016 and 2023 image classification models, respectively.

Table 3 .
NDVI change detection statistics for Gelephu.

Table 3 .
NDVI change detection statistics for Gelephu.