Preprint
Article

Utilizing Sentinel-2 for Recent Urbanization Trends on LULC Using a Semi-Automatic RF Classifier and Vegetation Change Dynamics for Gelephu, Bhutan

Altmetrics

Downloads

170

Views

72

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

28 November 2023

Posted:

28 November 2023

You are already at the latest version

Alerts
Abstract
Gelephu, located in the Himalayan region, has experienced one of the strongest development activities due to the suitable topography and geographic location, which has led to rapid urbanization. The main objective of the study is to perform land use land cover (LULC) mapping for Gelephu for 2016 and 2023 by a Random Forest (RF) classifier, using the Semi-automatic Classification Plugin (SCP) in QGIS and identify LULC changes. Furthermore, the study assessed the vegetation change dynamics within the study area by analysing the Normalized Difference Vegetation Index (NDVI) for 2016 and 2023. Additionally, the study 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 Sentiniel-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. This research provides the first version of LULC and NDVI maps of Gelephu. The analysis results indicate possible future implications on sustainable development management in Gelephu Bhutan.
Keywords: 
Subject: Environmental and Earth Sciences  -   Remote Sensing

1. 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 process 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, GIS-based image classification methods have made significant progress in LULC studies [13], [14], including the assessment and identification of changes 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 a 30-year period, 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 provide 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 allows the user to download several RS products, facilitates the processing of satellite imagery and automatically implements unsupervised and supervised classification by setting some parameters in the user interface. 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 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 [42]. To capitalize from 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 is 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 [43], vegetation degradation dynamics [44], and geo-environmental monitoring and land degradation assessment [45].
In summary, the purpose of the study is to perform efficient image classification using a GIS-based machine learning (RF) algorithm on SCP 7.10.11 using Sentinel-2 satellite imagery. This aims to generate reliable and precise LULC maps for Gelephu, while also extending the GIS and RS application to map NDVI and perform vegetation change assessments. The study mainly aimed to identify LULC and NDVI changes in Gelephu to reveal change dynamics that can be used for sustainable development plans and management. The study employed the open-source software QGIS 3.28.6.

2. Materials and Methods

2.1. Study Area

Gelephu (Figure 1) is one of the local government administrations (gewogs) in Sapang district and covers an area of 54.50 km2. According to the census 2017, the actual resident population is 4,461 and consists of 1,027 households[46]. Geographically, Gelephu is lies between, 25 to 27º north and 90 to 91º east with an elevation between 168 m in the southern floodplain and 1,799 m towards the north. Gelephu consists of one of the administrative units of the municipality, Gelephu Thromde, and covers an area of 10.66 km2. The LULC dynamics in the study area has 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.

2.2. Data Collection

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/. The Sentinel-2 images used in the current study for 2016 and 2023 corresponds to October 24, 2016, and January 21, 2013, 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 provide 13 spectral bands which consist of 10 m, 20 m, and 60 m spatial resolutions 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 near-infrared 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.

2.2.1. Data Projection and Processing

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. In particular, the colour composites such as band set 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.

2.2.2. 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 latest version of SCP 7.10.11 was used in the current study. The SCP plugin works semi-automatically and allows the user to download multiple RS products, facilitates the processing of satellite imagery and automatically implements unsupervised and supervised classification by setting some parameters in the user interface. 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 [47]. 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. SCP is a plugin for QGIS, designed specifically for image classification and land cover mapping which implements a semi-automatic approach where the user interacts with the software to guide the classification process. Instead of using a separate testing set, SCP uses the training data iteratively to improve the classification results as follows.
  • 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 a 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 allows the user to assess 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.
This iterative process continues until the user is satisfied with the classification results. The focus is on refining the training data and improving the classification model by incorporating user knowledge and feedback, usually preview of image classification scheme. Since the user is actively involved in the classification process, SCP does not require a separate testing set. 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.
In the current study, the process of developing a training dataset involved creating multiple ROI polygons for each specific land class. In particular, the utilization of multiple ROIs for each land class, combined with the RF classifier, aimed to improve the precision and reliability of the image classification outcomes. 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 the users 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 current study focused on image classification for four primary LULC classes namely water bodies, vegetation, built-up area and barren land. 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 classes. The RF classifier was used for two separate image classification schemes for the period between 2016 and 2023 and allowed for a comparative statistical analysis of the LULC change dynamics that occurred over a seven-year period.

2.2.3. NDVI Mapping

To further delve into the vegetation monitoring tools, vegetation change assessment was carried out separately using NDVI. NDVI serve a specific purpose of monitoring vegetation density and health status, while LULC focuses on categorising and classifying the types of land use activities and surface cover and indicating the spatial distribution of various land features.
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 Normalized Difference Vegetation Index (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 GIS software, 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 bend 8 (NIR) and band 4 (red), we obtained NDVI values for each pixel, which were then mapped to visualize the spatial distribution of vegetation in Gelephu.

3. Results and Discussion

3.1. LULC Mapping and Change Detection

Gelephu has seen significant LULC changes in recent years. It is evident that a variety of factors are major contributors to the rapid urbanization including migration, which played a crucial role in determining Gelephu's LULC dynamics. With the growing population and increasing economic activity, there has been a remarkable transformation of agricultural and rural areas into urbanized ones. 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 comes 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 clearing of vegetation and remodelling of surrounding areas to accommodate transport networks, in addition to construction of buildings and other service infrastructures. 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.
Figure 5 shows the spatial distribution and visual representation of the LULC changes that have occurred over a seven-year period. Statistical data was extracted from the image classification reports and a comparative LULC analysis was performed for Gelephu in 2016 and 2023. 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) indicate 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. 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. 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 in Figure 4 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.
Additionally, change detection on the LULC was performed (see Figure 5) using the postprocessing tool 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 5a. The LULC change detection map is further rendered in Figure 5b to visually depict the delineation of land-use class change spots and represent the actual land-use transition through the year 2023. In Figure 5b, white background is assigned to the land-use class to represent the no change land-use class.
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 6. The Gelephu Thromdes 2023 LULC map in Figure 6b shows a relatively high concentration of construction activities resulting from the construction of buildings and road networks compared to the 2016 LULC map. The change in vegetation cover up to 2023 is clearly pronounced in the land use class, suggesting a widespread shift in land-use dynamics toward barren land and built-up areas, including minor observations of water body cover.

3.2. Accuracy Assessment of LULC Maps

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 (Figure 7) 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 10.8.2 was employed to construct the ROC curves. The LULC raster files are imported into ArcGIS and the testing sample points 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.
Therefore, a receiver operating characteristic curve (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 class). The TPR and FPR is computed using the equations Eq. 1 and Eq. 2, 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 (Figure 8). 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:
TPR = TP / (TP + FN)
Similarly, the False Positive Rate (FPR) is defined as:
FPR = FP / (FP + TN)
The area under the ROC curve, commonly referred to as 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. The RF classification model resulted in 0.883 AUC for 2016 LULC mapping (Figure 8a) and 0.890 AUC for 2023 (Figure 8b) 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 time periods suggests that it is reliable and can effectively handle temporal changes in the landscape.

3.3. LULC statistics

The LULC statistics for Gelephu are presented in the Table 1. In 2016 the area covered by vegetation was 44.31 km2 making up 81.30% of the total area. However, by 2023 the vegetation area had shrunk to 37.62 km2, 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 km2 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 km2 indicating a decline of roughly 0.91%. This change can be attributed to alterations in the watercourse itself. Consequently, the proportion of water coverage to the area decreased to approximately just, around 0.04%. The built-up area in 2016 was 4.03 km2, representing 7.39% of the total area. By 2023, the built-up area increased to 7.08 km2, showing a significant increase of 5.65%. The portion of built-up area in the total area increased to 13.04%. Likewise, the barren land coverage was 5.64 km2 in 2016, accounting for 10.35% of the total area. The barren land area increased to 9.55 km2 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 area and barren land increased significantly.
Additionally, the Table 1 lists the LULC statistics for Gelephu Thromde. The LULC statistics between 2016 and 2023 shows 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 km2, 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 km2 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 km2, 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 km2, with the percent coverage increasing from 16.47% to 31.23%.
The LULC change detection results from Gelephu show several significant transitions (Figure 9). The largest land cover class, vegetation, was fairly stable, with an area of 34.95 km2. This provides information about the notable change of the vegetation areas over just seven-year study period. The changes between water and vegetation are relatively small. The water to vegetation transition covers an area of 0.04 km2, 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 km2, indicating a slight decrease in vegetation cover and an increase in water bodies. A total area of 2.62 km2 shows a transition from barren land to vegetation. The conversion of vegetation into built-up areas covers an area of 4.82 km2, indicating urban expansion and reflecting the conversion of natural vegetation into built-up infrastructure, such as residential, commercial, or industrial facilities. Conversion of vegetation to barren land covers an area of 4.48 km2, suggesting possible land degradation or land-use changes. This transitions between barren land and built-up areas illustrate urban expansion into previously unproductive or degraded areas. The conversion of barren land to built-up areas covers 1.16 km2, indicating land utilization for urban development purposes. Conversely, the conversion of developed land to barren land covers an area of 1.34 km2. The barren land cover class remained relatively stable, with an area of 3.38 km2, indicating the persistence of unproductive or degraded areas within the study area.

3.4. NDVI Calculations, Change Detection and Statistics

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, 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. 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 aim of the reclassification was to provide a more meaningful representation of the vegetation density and the dynamics of change.
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 is identified and analysed.
Analysis of the data for NDVI change detection between 2016 and 2023 (Table 2) shows significant shifts in land cover types within the study area. There was a significant decrease in water bodies: the area decreased from 0.96 to 0.06 km2 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 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 km2, 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 km2 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 km2, 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 km2, 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. In summary, the NDVI change statistics in Table 3 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 on the context of the impacts of climate change.

3.5. Validation

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 foot print 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 foot print 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 (Figure 12). Analysis revealed that a significant majority of the building shapes perfectly aligned with the pixels classified as built-up areas on the LULC map. This agreement between both datasets provided evidence regarding the accuracy and dependability of the LULC map generated by utilizing RF classifier. By validating the LULC map through comparison, with building shapes it demonstrates 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.

4. 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 a decline in percent cover from 81.30% to 69.31%. The built-up area, on the other hand, showed significant growth and increased by 5.65% from 7.39% to 13.04%. In addition, the barren land area saw a significant increase of 7.25%, with the percentage cover increasing from 10.35% to 17.60%. Likewise, the analysis of land cover changes in Gelephu Thromde between 2016 and 2023 shows significant changes in the study area. Key results indicate a significant decrease in vegetation, with a 29.93% decrease in area and a corresponding decrease in percentage coverage from 67.04% to 37.11%. The built-up area and barren land saw momentous growth, increasing by 15.05% and 14.76%, respectively. Significant growth in built-up areas and barren land indicates rapid urbanization over a seven-year period and may pose future environmental and land management challenges.
The analysis of NDVI change detection data between 2016 and 2023 in the study area reveals significant shifts in vegetation cover. The category of "No vegetation" increased by 9.56%, "Low vegetation" expanded by 41.24%, "Healthy vegetation" grew by 25.95%, and "Very healthy vegetation" decreased by 75.11%. These findings highlight the dynamic changes in vegetation cover and the potential ecological implications within the study area.

Author Contributions

Conceptualization, K.T. and M.I.; methodology, K.T.; software, K.T.; validation, K.T., A.A. and T.; formal analysis, K.T.; investigation, M.I.; resources, M.I. and A.A.; data curation, T.; writing—original draft preparation, K.T.; writing—review and editing, A.A., M.I., T. and K.T.; visualization, K.T.; supervision, M.I.; project administration, M.I. and K.T.; funding acquisition, M.I. All authors have read and agreed to the published version of the manuscript.

Funding

The APC was funded by University of South Africa.

Data Availability Statement

All data generated or analyzed during this study are included in this published article.

Acknowledgments

The author thanks the Office of the Vice-Chancellor of the Royal University of Bhutan (RUB), Bhutan, for generous support of the research through the Annual University Research Grant (AURG).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Pandey, P.C., Koutsias, N., Petropoulos, G.P., Srivastava, P.K., Ben Dor, E.: Land use/land cover in view of earth observation: data sources, input dimensions, and classifiers—a review of the state of the art. Geocarto Int. 36, 957–988 (2021). [CrossRef]
  2. Sam, S.C., Balasubramanian, G.: Spatiotemporal detection of land use/land cover changes and land surface temperature using Landsat and MODIS data across the coastal Kanyakumari district, India. Geod Geodyn. 14, 172–181 (2023). [CrossRef]
  3. Chen, S., Guo, Q., Li, L.: Sustainable Land Use Dynamic Planning Based on GIS and Symmetric Algorithm. Advances in Civil Engineering. 2022, (2022). [CrossRef]
  4. Mishra, P.K., Rai, A., Rai, S.C.: Land use and land cover change detection using geospatial techniques in the Sikkim Himalaya, India. Egyptian Journal of Remote Sensing and Space Science. 23, 133–143 (2020). [CrossRef]
  5. Sari, I.L., Weston, C.J., Newnham, G.J., Volkova, L.: Land cover modelling for tropical forest vulnerability prediction in Kalimantan, Indonesia. Remote Sens Appl. 32, 101003 (2023). [CrossRef]
  6. Thanh Noi, P., Kappas, M.: Comparison of Random Forest, k-Nearest Neighbor, and Support Vector Machine Classifiers for Land Cover Classification Using Sentinel-2 Imagery. Sensors (Basel). 18, 1–20 (2017). [CrossRef]
  7. Chaminé, H.I., Pereira, A.J.S.C., Teodoro, A.C., Teixeira, J.: Remote sensing and GIS applications in earth and environmental systems sciences. SN Appl Sci. 3, 2–4 (2021). [CrossRef]
  8. Talukdar, S., Singha, P., Mahato, S., Pal, S.: Land-Use Land-Cover Classification by Machine Learning Classifiers for Satellite Observations—A Review. (2020). [CrossRef]
  9. Hosseiny, B., Abdi, A.M., Jamali, S.: Urban land use and land cover classification with interpretable machine learning – A case study using Sentinel-2 and auxiliary data. Remote Sens Appl. 28, 100843 (2022). [CrossRef]
  10. Campos-Taberner, M., García-Haro, F.J., Martínez, B., Izquierdo-Verdiguier, E., Atzberger, C., Camps-Valls, G., Gilabert, M.A.: Understanding deep learning in land use classification based on Sentinel-2 time series. Sci Rep. 10, 1–12 (2020). [CrossRef]
  11. Sertel, E., Ekim, B., Ettehadi Osgouei, P., Kabadayi, M.E.: Land Use and Land Cover Mapping Using Deep Learning Based Segmentation Approaches and VHR Worldview-3 Images. Remote Sens (Basel). 14, 1–20 (2022). [CrossRef]
  12. Viet Du, Q.V., Nguyen, H.D., Pham, V.T., Nguyen, C.H., Nguyen, Q.H., Bui, Q.T., Doan, T.T., Tran, A.T., Petrisor, A.I.: Deep learning to assess the effects of land use/land cover and climate change on landslide susceptibility in the Tra Khuc river basin of Vietnam. Geocarto Int. 38, (2023). [CrossRef]
  13. Abdi, A.M.: Land cover and land use classification performance of machine learning algorithms in a boreal landscape using Sentinel-2 data. GIsci Remote Sens. 57, 1–20 (2020). [CrossRef]
  14. Belenok, V., Noszczyk, T., Hebryn-Baidy, L., Kryachok, S.: Investigating anthropogenically transformed landscapes with remote sensing. Remote Sens Appl. 24, 100635 (2021). [CrossRef]
  15. Pasha, S.V., Reddy, C.S., Jha, C.S., Rao, P.V.V.P., Dadhwal, V.K.: Assessment of Land Cover Change Hotspots in Gulf of Kachchh, India Using Multi-Temporal Remote Sensing Data and GIS. Journal of the Indian Society of Remote Sensing. 44, 905–913 (2016). [CrossRef]
  16. Duarte, L., Teodoro, A.C., Cunha, M.: A semi-automatic approach to derive land cover classification in soil loss models. In: Proc. SPIE 11156, Earth Resources and Environmental Remote Sensing/GIS Applications X, 111560B. pp. 1–13 (2019). [CrossRef]
  17. Tempa, K., Aryal, K.R.: Semi-automatic classification for rapid delineation of the geohazard-prone areas using Sentinel-2 satellite imagery. SN Appl Sci. 4, (2022). [CrossRef]
  18. Abida, K., Barbouchi, M., Boudabbous, K., Toukabri, W., Saad, K., Bousnina, H., Sahli Chahed, T.: Sentinel-2 Data for Land Use Mapping: Comparing Different Supervised Classifications in Semi-Arid Areas. Agriculture (Switzerland). 12, (2022). [CrossRef]
  19. Benhammou, Y., Alcaraz-Segura, D., Guirado, E., Khaldi, R., Achchab, B., Herrera, F., Tabik, S.: Sentinel2GlobalLULC: A Sentinel-2 RGB image tile dataset for global land use/cover mapping with deep learning. Sci Data. 9, 1–20 (2022). [CrossRef]
  20. Xu, F., Heremans, S., Somers, B.: Urban land cover mapping with Sentinel-2: a spectro-spatio-temporal analysis. Urban Informatics. 1, 1–18 (2022). [CrossRef]
  21. Yousefi, S., Mirzaee, S., Almohamad, H., Dughairi, A.A. Al, Gomez, C., Siamian, N., Alrasheedi, M., Abdo, H.G.: Image Classification and Land Cover Mapping Using Sentinel-2 Imagery: Optimization of SVM Parameters. Land (Basel). 11, (2022). [CrossRef]
  22. Bouaziz, M., Eisold, S., Guermazi, E.: Semiautomatic approach for land cover classification: a remote sensing study for arid climate in southeastern Tunisia. EuroMediterr J Environ Integr. 2, 1–7 (2017). [CrossRef]
  23. Wang, Q., Blackburn, G.A., Onojeghuo, A.O., Dash, J., Zhou, L., Zhang, Y., Atkinson, P.M.: Fusion of Landsat 8 OLI and Sentinel-2 MSI Data. IEEE Transactions on Geoscience and Remote Sensing. 55, 3885–3899 (2017). [CrossRef]
  24. Belgiu, M., Csillik, O.: Sentinel-2 cropland mapping using pixel-based and object-based time-weighted dynamic time warping analysis. Remote Sens Environ. 204, 509–523 (2018). [CrossRef]
  25. Ge, W., Cheng, Q., Tang, Y., Jing, L., Gao, C.: Lithological classification using Sentinel-2A data in the Shibanjing ophiolite complex in Inner Mongolia, China. Remote Sens (Basel). 10, (2018). [CrossRef]
  26. Gao, J., Liu, J., Liang, T., Hou, M., Ge, J., Feng, Q., Wu, C., Li, W.: Mapping the forage nitrogen-phosphorus ratio based on Sentinel-2 MSI data and a random forest algorithm in an alpine grassland ecosystem of the Tibetan Plateau. Remote Sens (Basel). 12, (2020). [CrossRef]
  27. Azzari, G., Jain, S., Jeffries, G., Kilic, T., Murray, S.: Understanding the requirements for surveys to support satellite-based crop type mapping: Evidence from sub-saharan Africa. Remote Sens (Basel). 13, (2021). [CrossRef]
  28. Schulthess, U., Rodrigues, F., Taymans, M., Bellemans, N., Bontemps, S., Ortiz-Monasterio, I., Gérard, B., Defourny, P.: Optimal Sample Size and Composition for Crop Classification with Sen2-Agri’s Random Forest Classifier. Remote Sens (Basel). 15, (2023). [CrossRef]
  29. Rangzan, K., Kabolizadeh, M., Mousavi, S.S., Karimi, D., Rashnoei, A.: Assessing the potential of Sentinel-2 imagery and spectroscopy for determining the origin of ancient artifacts in Khuzestan, Iran. Egyptian Journal of Remote Sensing and Space Science. 26, 455–476 (2023). [CrossRef]
  30. Congedo, L.: Semi-Automatic Classification Plugin: A Python tool for the download and processing of remote sensing images in QGIS. J Open Source Softw. 6, 3172 (2021). [CrossRef]
  31. QGIS Development Team: QGIS geographic information system, (2021).
  32. Furukawa, F., Morimoto, J., Yoshimura, N., Kaneko, M.: Comparison of Conventional Change Detection Methodologies Using High-Resolution Imagery to Find Forest Damage Caused by Typhoons. Remote Sens (Basel). 12, 3242 (2020). [CrossRef]
  33. Höhle, J.: Automated mapping of buildings through classification of DSM-based ortho-images and cartographic enhancement. International Journal of Applied Earth Observation and Geoinformation. 95, 102237 (2021). [CrossRef]
  34. Whiteside, T.G., Boggs, G.S., Maier, S.W.: Comparing object-based and pixel-based classifications for mapping savannas. International Journal of Applied Earth Observation and Geoinformation. 13, 884–893 (2011). [CrossRef]
  35. Xie, S., Liu, L., Zhang, X., Yang, J., Chen, X., Gao, Y.: Automatic land-cover mapping using landsat time-series data based on google earth engine. Remote Sens (Basel). 11, (2019). [CrossRef]
  36. Gonzalez-ollauri, A., Mickovski, S.B.: A simple GIS-based tool for the detection of landslide-prone zones. 10, 1–15 (2021). [CrossRef]
  37. Gašparović, M., Zrinjski, M., Gudelj, M.: Automatic cost-effective method for land cover classification (ALCC). Comput Environ Urban Syst. 76, 1–10 (2019). [CrossRef]
  38. Arekhi, M., Goksel, C., Sanli, B.F., Senel, G.: Comparative Evaluation of the Spectral and Spatial Consistency of Sentinel-2 and Landsat-8 OLI Data for Igneada Longos Forest. International Journal of Geo-Information. 8, 1–13 (2019). [CrossRef]
  39. Gilani, H., Shrestha, H.L., Murthy, M.S.R., Phuntso, P., Pradhan, S., Bajracharya, B., Shrestha, B.: Decadal land cover change dynamics in Bhutan. J Environ Manage. 148, 91–100 (2015). [CrossRef]
  40. FRMD: Land Use and Land Cover Assessment of Bhutan2016, Technical Report. , Thimphu, Bhutan (2017).
  41. FRMD: Land use and land cover of Bhutan 2016, Maps and Statistics. , Thimphu, Bhutan (2017).
  42. Wang, S.W., Munkhnasan, L., Lee, W.K.: Land use and land cover change detection and prediction in Bhutan’s high altitude city of Thimphu, using cellular automata and Markov chain. Environmental Challenges. 2, (2021). [CrossRef]
  43. Ghebrezgabher, M.G., Yang, T., Yang, X., Eyassu Sereke, T.: Assessment of NDVI variations in responses to climate change in the Horn of Africa. Egyptian Journal of Remote Sensing and Space Science. 23, 249–261 (2020). [CrossRef]
  44. Fokeng, R.M., Fogwe, Z.N.: Landsat NDVI-based vegetation degradation dynamics and its response to rainfall variability and anthropogenic stressors in Southern Bui Plateau, Cameroon. Geosystems and Geoenvironment. 1, 100075 (2022). [CrossRef]
  45. Kumar, B.P., Babu, K.R., Anusha, B.N., Rajasekhar, M.: Geo-environmental monitoring and assessment of land degradation and desertification in the semi-arid regions using Landsat 8 OLI / TIRS, LST, and NDVI approach. Environmental Challenges. 8, 1–11 (2022). [CrossRef]
  46. Gelegphu, http://www.sarpang.gov.bt/gewogs/gelegphu, last accessed 2023/11/15.
  47. Piao, Y., Jeong, S., Park, S., Lee, D.: Analysis of land use and land cover change using time-series data and random forest in north korea. Remote Sens (Basel). 13, (2021). [CrossRef]
Figure 1. Study area (a) Geographical location of Bhutan; (b) Digital elevation model of the study area.
Figure 1. Study area (a) Geographical location of Bhutan; (b) Digital elevation model of the study area.
Preprints 91635 g001
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. 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.
Preprints 91635 g002
Figure 3. Colour composites of the study area (a) Band set 4-3-2 (true colour); (b) 5-6-2 (healthy vegetation); (c) 5-4-3 (colour infrared, vegetation); (d) 7-6-4 (false colour, urban).
Figure 3. Colour composites of the study area (a) Band set 4-3-2 (true colour); (b) 5-6-2 (healthy vegetation); (c) 5-4-3 (colour infrared, vegetation); (d) 7-6-4 (false colour, urban).
Preprints 91635 g003
Figure 4. LULC maps of Gelephu (a) LULC map of 2016; (b) LULC map of 2023.
Figure 4. LULC maps of Gelephu (a) LULC map of 2016; (b) LULC map of 2023.
Preprints 91635 g004
Figure 5. LULC change detection map of Gelephu.
Figure 5. LULC change detection map of Gelephu.
Preprints 91635 g005
Figure 6. Subset of Gelephu LULC map (a) Gelephu Thromde LULC map of 2016; (b) Gelephu Thromde LULC map of 2023.
Figure 6. Subset of Gelephu LULC map (a) Gelephu Thromde LULC map of 2016; (b) Gelephu Thromde LULC map of 2023.
Preprints 91635 g006
Figure 7. Testing data set.
Figure 7. Testing data set.
Preprints 91635 g007
Figure 8. Model accuracy assessment. (a) ROC for 2016 LULC, (b) ROC for 2023 LULC
Figure 8. Model accuracy assessment. (a) ROC for 2016 LULC, (b) ROC for 2023 LULC
Preprints 91635 g008
Figure 9. LULC transition pattern for Gelephu
Figure 9. LULC transition pattern for Gelephu
Preprints 91635 g009
Figure 10. NDVI maps for Gelephu. (a) 2016 NDVI, (b) 2023 NDVI
Figure 10. NDVI maps for Gelephu. (a) 2016 NDVI, (b) 2023 NDVI
Preprints 91635 g010
Figure 11. Reclassified NDVI maps; (a) 2016 NDVI, (b) 2023 NDVI
Figure 11. Reclassified NDVI maps; (a) 2016 NDVI, (b) 2023 NDVI
Preprints 91635 g011
Figure 12. Validation maps. (a) South, (b) North-east
Figure 12. Validation maps. (a) South, (b) North-east
Preprints 91635 g012
Table 1. Summary of LULC statistics for Gelephu and Gelephu Thromde respectively.
Table 1. Summary of LULC statistics for Gelephu and Gelephu Thromde respectively.
Class Gelephu Gelephu Thromde
Area (km2) % Coverage %
Change
Area (km2) % Coverage % Change
2016 2023 2016 2023 2016 2023 2016 2023
Vegetation 44.31 37.62 81.30 69.31 -11.99 7.15 3.96 67.04 37.11 -29.93
Water 0.52 0.02 0.95 0.04 -0.91 0.03 0.04 0.25 0.37 0.12
Built-up area 4.03 7.08 7.39 13.04 5.65 1.73 3.34 16.24 31.28 15.05
Barren land 5.64 9.55 10.35 17.60 7.25 1.76 3.33 16.47 31.23 14.76
Total 54.50 54.27 100.00 100.00 10.66 10.66 100.00 100.00
Table 2. NDVI change detection statistics for Gelephu
Table 2. NDVI change detection statistics for Gelephu
Class Area (km2) % Coverage % Change
2016 2023 2016 2023
Water body 0.96 0.06 1.76 0.11 -1.65
No vegetation 2.87 8.08 5.27 14.83 9.56
Low vegetation 2.07 24.54 3.80 45.04 41.24
Healthy vegetation 7.13 21.27 13.09 39.04 25.95
Very healthy vegetation 41.46 0.53 76.09 0.98 -75.11
Total 54.50 54.50
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated