# Evaluation of Spectral Indices and Global Thresholding Methods for the Automatic Extraction of Built-Up Areas: An Application to a Semi-Arid Climate Using Landsat 8 Imagery

Yassine Harrak, Ahmed Rachid, Rahim Aguejdad

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Yassine Harrak, Ahmed Rachid, Rahim Aguejdad. Evaluation of Spectral Indices and Global Thresholding Methods for the Automatic Extraction of Built-Up Areas: An Application to a Semi-Arid Climate Using Landsat 8 Imagery. Urban Science, 2025, 9 (3), pp.78. 10.3390/urbansci9030078. hal-04988707

HAL Id: hal-04988707

https://hal.science/hal-04988707v1

Submitted on 12 Mar 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
Article

# Evaluation of Spectral Indices and Global Thresholding Methods for the Automatic Extraction of Built-Up Areas: An Application to a Semi-Arid Climate Using Landsat 8 Imagery

Yassine Harrak <fnref n="1" />,*, Ahmed Rachid <fnref n="1" /> and Rahim Aguejdad <fnref n="2" />,†

<figure><img src="image_3.png" /></figure>

check for
updates

Academic Editors: Chiara Chiatti and
Roberta Jacoby Cureau

Received: 24 December 2024

Revised: 17 February 2025

Accepted: 7 March 2025

Published: 11 March 2025

**Citation:** Harrak, Y.; Rachid, A.; Aguejdad, R. Evaluation of Spectral Indices and Global Thresholding Methods for the Automatic Extraction of Built-Up Areas: An Application to a Semi-Arid Climate Using Landsat 8 Imagery. *Urban Sci.* 2025, 9, 78.
https://doi.org/10.3390/urbansci9030078

**Copyright:** © 2025 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license
(https://creativecommons.org/
licenses/by/4.0/).

**Abstract:** The rapid expansion of built-up areas (BUAs) requires effective spatial and temporal monitoring, being a crucial practice for urban land use planning, resource allocation, and environmental studies, and spectral indices (SIs) can provide efficiency and reliability in automating the process of BUAs extraction. This paper explores the use of nine spectral indices and sixteen thresholding methods for the automatic mapping of BUAs using Landsat 8 imagery from a semi-arid climate in Morocco during spring and summer. These indices are the Normalized Difference Built-Up Index (NDBI), the Vis-red-NIR Built-Up Index (VrNIR-BI), the Perpendicular Impervious Surface Index (PISI), the Combinational Biophysical Composition Index (CBCI), the Normalized Built-up Area Index (NBAI), the Built-Up Index (BUI), the Enhanced Normalized Difference Impervious Surfaces Index (ENDISI) and the Built-up Land Features Extraction Index (BLFEI). Results show that BLFEI, SWIRED, and BUI maintain high separability between built-up and each of the other land cover types across both seasons, as evaluated via the Spectral Discrimination Index (SDI). The lowest SDI values for all three indices were observed for bare soil against BUAs, with BLFEI recording 1.21 in the wet season and 1.05 in the dry season, SWIRED yielding 1.22 and 1.08, and BUI showing 1.21 and 1.08, demonstrating their robustness in distinguishing BUAs from other land covers under varying phenological and soil moisture conditions. These indices reached overall accuracies of 93.97%, 93.39% and 92.81%, respectively, in wet conditions, and 91.57%, 89.17% and 89.67%, respectively, in dry conditions. The assessment of thresholding methods reveals that the Minimum method resulted in the highest accuracies for these indices in wet conditions, where bimodal medium peaked histograms were observed, whereas the use of Li, Huang, Shanbhag, Otsu, K-means, or IsoData was found to be the most effective under dry conditions, where more peaked histograms were observed.

**Keywords:** land use land cover (LULC); built-up; spectral index; seasonal variation; automatic thresholding; Landsat 8

## 1. Introduction

The increasing demand for housing and infrastructure due to population growth and urbanization is reshaping the environment at an unprecedented pace. Natural landscapes

<footnote n="1">¹ Laboratory of Physics and Chemistry of Processes and Materials, Faculty of Sciences and Techniques, Hassan First University, Settat 26002, Morocco; y.harrak@uhp.ac.ma (Y.H.); ahmed.rachid@uhp.ac.ma (A.R.)</footnote>

<footnote n="2">² Unité Mixte de Recherche (Joint Research Unit) Territoire Environnement Télédetection Information Spatiale, 34000 Montpellier, France</footnote>

<footnote>* Correspondence: rahim.aguejdad@cnrs.fr</footnote>
are being replaced by anthropogenic impervious surfaces (ISAs), a rapidly expanding phenomenon that is exerting profound impacts on ecosystems worldwide and usually implies the alteration of land use patterns by reducing agricultural and green spaces, which challenges the equitable distribution of water, energy and transportation networks [1,2]. This expansion is also associated with multiple environmental complications in the cityscape, causing more vulnerability to water and air pollution, stormwater runoff and flooding, urban heat islands, and more [3–5]. Therefore, the accurate and adaptive spatial and temporal monitoring of built-up areas (BUAs) has become presently essential for city planners to effectively manage urban expansion and support the development of smarter cities through optimized management of population growth, transportation networks, and public services accordingly, while at the same time enhancing climate resilience in the face of a rapidly changing climate [6–8]. By leveraging technological innovation and real-time monitoring systems with climate-conscious urban planning strategies, cities can achieve sustainable growth while mitigating climate risks.

BUAs can be defined as human-made surfaces that prevent water infiltration into the soil, such as concrete buildings, roads, streets, and parking lots [9]. The emergence of remote sensing and the launch of earth observation satellite missions marked a new era in land use and land cover (LULC) monitoring due to its spatially comprehensive, precise, and time-cost-effective approach [10,11]. However, still, the extraction of BUAs based on remote sensing data is still not easy, considering that these areas usually have heterogeneous spectral properties, which leads to sensors capturing the electromagnetic radiation of a mixture of land covers where similarities with bare land and water bodies in particular are considered a big challenge [12,13]. As a result, methods and techniques employed for this task have experienced substantial improvements over the years. In the early times of remote sensing, the visual interpretation of early satellite images such as those from Landsat-1 was the primary method [14]. The progress in computation and image processing unveiled new prospects in the development of multiple classification algorithms such as the Maximum Likelihood Classifier (MLC) [15], and the establishment of image differencing techniques from multi-spectral imagery such as Object-Based Image Analysis (OBIA) [16], Spectral Mixture Analysis (SMA) [17] and Spectral Indices (SIs) [18,19]. Nevertheless, mapping BUAs from satellite imagery was extensively addressed by means of semantic segmentation using pixel-based supervised and unsupervised methods in the past few years, and machine learning algorithms such as decision trees, Support Vector Machines (SVM), clustering algorithms and more recently Artificial Neural Networks (ANN) have further revolutionized LULC classification in general [20–22]. Even though these techniques are fundamental in LULC studies, they can show several limitations regarding the automation process. Supervised classification necessitates the continuous and thorough collection of training data, requiring a considerable amount of time and expert knowledge. Correspondingly, unsupervised classification also requires cluster interpretation and labeling [23]. These limitations can be avoided with the help of SIs, serving as practical and easy-to-implement tools for separating land covers based on their spectral properties [24].

A wide range of SIs for mapping BUAs from multi-spectral imagery are presented in the literature and are usually categorized following the regions of the spectrum they use within their mathematical formulation. The Vis-red/green-NIR Built-up Indices (VrNIR-BI and VgNIR-BI) [25], the Perpendicular Impervious Surface Index (PISI) [26] and the Combinational Biophysical Composition Index (CBCI) [27] for example use the visible and Near Infra Red (NIR) bands, whereas indices like the Normalized Difference Built-Up Index (NDBI) [28], the Normalized Built-up Area Index (NBAI) [29], the Built-Up Index (BUI) [30], the Enhanced Normalized Difference Impervious Surfaces Index (ENDISI) [31] and the
Built-up Land Features Extraction Index (BLFEI) [32] all use the Short Wave Infra Red (SWIR) bands additionally and are the most common for BUAs and ISAs extraction. Other forms of SIs include the use of Panchromatic and Thermal bands [33,34], Nighttime Light datasets [35] or other SIs combining these methods. To the extent of classifying pixels into built-up and non-built-up from index-based maps, the most common and simple techniques in image processing are thresholding methods [36]. The resulting continuous pixel maps of SIs are divided into a foreground and a background, enabling the construction of masks highlighting BUAs. In contrast to adaptive thresholding, where multiple thresholds are computed for different regions of an image, global thresholding is a straightforward operation that ensures consistency and computational efficiency and requires no parameter passing, which can contribute to the automation outcome in the process of extracting BUAs. Many of these methods are documented in existing research and are used in the context of SIs for various LULC applications, including histogram-based methods such as Otsu's method and clustering-based methods such as K-means, among others [31,37–39].

The evaluation of SIs for the extraction of BUAs has been the focus of numerous studies, with a number of indices and earth observation data examined in different climate conditions and seasons, utilizing various thresholding methods. Li et al. [40] and Almohamad and Alshwesh [41] have identified some indices such as PISI and CBCI to be effective both in dry and wet conditions, although the former study was conducted in a humid subtropical region and the latter in dry arid climate, conversely, some indices were found susceptible of inferring BUAs from bare land and water, such as the Index-based Built-up Index (IBI), ENDISI and NDBI. However, Chen et al. [42] found ENDISI particularly robust in various geographical regions, and Sekertekin et al. [43] had better results using IBI than other indices. This indicates that the choice of SIs can be highly circumstantial. In addition, Landsat and Sentinel have emerged as the two main satellite platforms for built-up and ISAs extraction, Xi et al. [44] showed that both data sources provided comparable results despite differences in spatial resolutions and spectral response but found Sentinel to be better suited to finer urban features. The same authors also emphasized the effectiveness of using SWIR bands, which are available in both satellites. Some of the mentioned studies also tested a number of thresholding methods, another critical aspect of the extraction process. Li et al. [40] and Almohamad and Alshwesh [41] concluded that both manual thresholding and global automated methods were effective, with manual methods offering slightly better stability and accuracy; however, automated methods like IsoData were also found to provide consistency and are especially valuable for large-scale or real-time monitoring applications.

Despite the extensive use of SIs in multiple BUA and ISA-related applications, their performance and suitability to climate conditions of semi-arid regions as to both seasonal variation and thresholding techniques are not well explored. In semi-arid regions where phenology, soil moisture content, and atmospheric conditions are closely tied to rainfall patterns, climate change is manifested by altering precipitation regimes, increasing temperatures, and intensifying droughts. Consequently, spectral responses of vegetation, soil, impervious surfaces, and, therefore, spectral indices become more dynamic and less predictable [45–47]. Given that the effectiveness of SIs is highly circumstantial, investigating the performance of SIs in such regions provides insights into their reliability to pronounced phenological shifts. Moreover, the effectiveness of indices also relies on setting appropriate thresholds as it can lead to alteration in accuracies [48]. To address these issues, we tested in this study the performance of nine SIs (NDBI, VRNIR, SWIRED, BUI, ENDISI, PISI, CBCI, BLFEI, and NBAI) using Landsat 8 imagery captured during the spring and summer of 2021 in the municipality of Settat in Morocco, then we subjected the computed spectral maps to sixteen global thresholding methods. The main purpose of this evaluation is to support the
choice of SIs based on their ability to differentiate between BUAs and each of the other land covers, and to quantify their adaptability to changes in phenology and soil moisture in accordance with the spectral response of each land cover. This study also aims to support the choice of global thresholding methods, being the main tools for automatically extracting BUAs from index-based maps.

## 2. Materials and Methods

### 2.1. Materials

#### 2.1.1. Study Area

Our study was carried out in the municipality of Settat, the capital of Settat province in Morocco. It is in the central-eastern part of the Casablanca-Settat region, which is the most populated region in the country. The municipality is positioned at 33°18′–33°57′ N and 7°33′–7°48′ W. It occupies an area of 54.67 km² and is characterized by a semi-arid climate with hot and dry summers, 'BSh' according to the Köppen–Geiger classification [49]. Precipitations in the province are distributed generally in the rainy period from October to April, with maximum values fluctuating around 50 mm for December and January, averaging an annual rainfall of 300 mm. Mean temperatures range between 5 °C and 15 °C during the rainy period and between 35 °C and 45 °C during the summer [50,51]. However, the country as a whole has experienced unprecedented successive drought years between the years 2019 and 2024 since at least the 1950s, with a registered rainfall deficit of 48% in 2023 compared to the normal mean over the period of 1981–2010 and 2023 being the warmest year ever on record in Moroccan history [52]. The city of Settat is located 370 m above sea level and is surrounded by mostly rain-fed agricultural foothills. Based on the latest Moroccan census conducted in September 2024 by the Higher Planning Commission (HCP), the municipality has a population of 171,556 against 142,250 in 2014, registering an annual growth of 1.7% [53].

#### 2.1.2. Datasets

Two Landsat 8 Collection 2 Level-2 Surface Reflectance scenes, courtesy of the U.S. Geological Survey (Path: 202, Row: 037), were downloaded from the earth explorer (https://earthexplorer.usgs.gov/ (accessed on 10 May 2024)). USGS Landsat collection 2 level 2 products are certified by the Committee on Earth Observation Satellites as Analysis Ready Products, which means that no additional correction was needed [54]. Table 1 shows details about the collected scenes. These scenes were chosen with consideration to cloud cover, phenology during spring and dryness during summer in the study area. The spectral bands used consist of 6 total bands: Blue (0.450–0.51 µm), Green (0.53–0.59 µm), Red (0.64–0.67 µm), NIR (0.85–0.88 µm), SWIR 1 (1.57–1.65 µm) and SWIR 2 (2.11–2.29 µm), all available at 30 m spatial resolution. Training and validation processes were based on Google Earth imagery, and the data used for the basic map of the study area (Figure 1) was based on the OpenStreetMap project.

<table><caption>**Table 1.** Details about Landsat 8 scenes used in this study.</caption><thead><tr><td>Date of Capture</td><td>Scene ID</td><td>Cloud Cover</td></tr></thead><tbody><tr><td>24 March 2021 (Spring)</td><td>LC08_L2SP_202037_20210324_20210402_02_T1</td><td>0.02%</td></tr><tr><td>14 July 2021 (Summer)</td><td>LC08_L2SP_202037_20210714_20210402_02_T1</td><td>0.21%</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Study area. (a) Casablanca-Settat region in Morocco; (b) Settat municipality and the province of Settat; (c) Basic map of Settat.</figcaption></figure>

## 2.2. Methods

The proposed methodology (Figure 2) consists of first performing a LULC classification by training and deploying a Support Vector Machine (SVM) classifier in order to extract accurate spectral profiles and to thoroughly examine the separability potential of each spectral index.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Flowchart of the followed methodology.</figcaption></figure>
This classification was also used to generate stratified random points for the accuracy assessment process. On the other hand, SIs were computed and fed to sixteen global thresholding algorithms, most of which are modified versions of Python 3.9.18 codes of the AutoThreshold plugin of ImageJ 1.42m freeware (https://imagej.net/plugins/auto-threshold (accessed on 20 August 2024)). Manual thresholding was also performed through a trial-and-error process for comparison purposes.

### 2.2.1. LULC Classification and Spectral Profiles

Land covers in the study area were divided into five main classes: built-up, bare land, trees, water, and agriculture. Here, agriculture refers to croplands or areas with grass or other types of vegetation except trees, and bare land refers to rock, sand, or soil with very sparse to no vegetation. Support vector machines are robust non-parametric machine learning methods used in multiple image-processing applications, including LULC classification, due to their ability to handle high-dimensional data and to ensure high classification accuracies with small amounts of training data [55,56].

Using the graphical user interface of ArcGIS Pro 3.2.2, we trained and deployed an SVM model on a composite image containing the total six OLI bands from the spring of 2021 scenes, as shown by Figure 3. Training data were collected in the form of polygons based on a geo-referenced high-resolution Google Earth image of the study area, acquired in February 2021. To assess the accuracy of both our model and SIs, 1210 points were generated in ArcGIS Pro using the stratified random method. They were then converted to a KML file and were labeled on Google Earth Pro 7.3.6 accordingly.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** SVM classification of LULC in the study area: (a) Collected training data; (b) LULC classification map; (c) Reference points for accuracy assessment.</figcaption></figure>

The confusion matrix of the classification is illustrated in Table 2, the model reached an overall accuracy of 93.88%. Spectral profiles of spring and summer scenes were constructed by averaging reflectance values for each band and for each LULC class, obtained after the application of scaling factors on digital numbers (DN). Figure 4 shows the resulting spectral profiles.
<table><caption>**Table 2.** Confusion matrix of SVM classification.</caption><thead><tr><th>Class</th><th>Agriculture</th><th>Built-Up</th><th>Bare Land</th><th>Trees</th><th>Water</th><th>UA (%)</th></tr></thead><tbody><tr><td>Agriculture</td><td>630</td><td>14</td><td>7</td><td>4</td><td>0</td><td>96.18</td></tr><tr><td>Built-up</td><td>3</td><td>286</td><td>3</td><td>0</td><td>0</td><td>97.95</td></tr><tr><td>Bare Land</td><td>20</td><td>17</td><td>196</td><td>0</td><td>0</td><td>84.12</td></tr><tr><td>Trees</td><td>4</td><td>0</td><td>2</td><td>14</td><td>0</td><td>70.00</td></tr><tr><td>Water</td><td>0</td><td>0</td><td>0</td><td>0</td><td>10</td><td>100.00</td></tr><tr><td>PA (%)</td><td>95.89</td><td>90.22</td><td>94.23</td><td>77.78</td><td>100.00</td><td></td></tr></tbody></table>

UA and PA are, respectively, the user's and the producer's accuracies.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Spectral profiles of surface reflectance in the study area: (a) Spring; (b) Summer.</figcaption></figure>

### 2.2.2. Spectral Indices

The evaluated SIs were chosen because of their different design approaches, the spectral properties of where they were developed, the purpose of deployment, and the level of complexity.

#### Normalized Difference Built-up Index (NDBI)

The Normalized Difference Built-Up Index [28] was first tested on Landsat TM sensor reflectance products to map urban areas in the city of Nanjing in eastern China. It is given by the following equation.

$$NDBI = \frac{SWIR1 - NIR}{SWIR1 + NIR} \qquad (1)$$

#### Vis-red/green-NIR Built-up Indices (VrNIR-BI and VgNIR-BI)

Vis-red-NIR and Vis-green-NIR indices use the visible red and green bands, respectively, in combination with the NIR band. They were proposed by Estoque and Murayama [25] and were tested on Landsat 7 and Landsat 8 imagery of Metro Manila in the Philippines. In this study, Vis-red-NIR (VrNIR) was tested. It can be calculated using Equation (2).

$$VrNIR = \frac{RED - NIR}{RED + NIR} \qquad (2)$$

#### Swir Red index (SWIRED)

Capolupo et al. [57] proposed a new algorithm based on two new simple SIs, namely the SwirTirRed (STRed) and the SWIRED index. STRed was used to detect water, mining
areas, and sparse and dense vegetation, while SWIRED was used to detect BUAs. It is computed using Equation (3).

$$SWIRED = \frac{SWIR1 - RED}{SWIR1 + RED} \qquad (3)$$

### Normalized Built-up Area Index (NBAI)

To address the problem of spectral similarities between bare soils and BUAs, Waqar et al. [29] proposed three indices, one of which is the Normalized Built-up Area Index (NBAI). These indices were derived from the spectral profiles of the city of Islamabad, Pakistan, and its reserved green spaces. In the case of NBAI, a combination of green, SWIR1, and SWIR2 bands was deduced as shown in Equation (4).

$$NBAI = \frac{SWIR2 - \left( \frac{SWIR1}{GREEN} \right)}{SWIR2 + \left( \frac{SWIR1}{GREEN} \right)} \qquad (4)$$

### Built-up Land Features Extraction Index (BLFEI)

BLFEI was designed and tested on OLI imagery by Bouhennache et al. [32], based on the spectral properties of the city of Algiers in Algeria. It uses the SWIR, red, and green bands, which resulted in good separability between BUAs and other land covers. It is calculated using the following formula:

$$BLFEI = \frac{\frac{GREEN+RED+SWIR2}{3} - SWIR1}{\frac{GREEN+RED+SWIR2}{3} + SWIR1} \qquad (5)$$

### Built-Up Index (BUI)

The Built-Up Index was developed using the red and both of the SWIR bands, it was tested on Landsat ETM+ products of Thessaloniki city in Greece. Accuracy assessment of BUI showed an Overall Accuracy (OA) of 90% [30]. It can be computed using the following equation:

$$BUI = \frac{2 \times (RED \times SWIR2) - (SWIR1 \times SWIR2)}{(RED + SWIR2) \times (SWIR1 + SWIR2)} \qquad (6)$$

### Combinational Biophysical Composition Index (CBCI)

CBCI (Equation (9)) was proposed by Zhang et al. [27] for the purpose of effectively highlighting the major biophysical composition of urban areas and improving the separability between ISAs and bare soil. It was tested on multiple types of satellite imagery, including Landsat 8 OLI products, in Guangzhou city in China. It is computed as follows:

$$MBSI = \frac{2 \times (RED - GREEN)}{RED - GREEN - 2} \qquad (7)$$

$$OSAVI = \frac{2 \times (NIR - RED)}{NIR - RED + 0.16} \qquad (8)$$

$$CBCI = (A + 1) \times MBSI - OSAVI + A \qquad (9)$$

The coefficient A in Equation (9) is a correction coefficient that was calculated by linear regression to increase the value of MBSI and decrease OSAVI. It is set to 0.51.

### Perpendicular Impervious Surface Index (PISI)

PISI [26] uses the inference between the blue and NIR bands. Perpendicular indices, in general, are derived by utilizing the geometry of features by plotting reflectance values of
spectral bands against each other and defining a line of reference with maximum variations. The perpendicular distance from points in the feature space to the reference line is used as the classification basis [58]. It was computed using Equation (10).

$$PISI = 0.8192 \times \text{BLUE} - 0.5735 \times \text{NIR} + 0.0750 \quad (10)$$

## Enhanced Normalized Difference Impervious Surfaces Index (ENDISI)

Considering the spectral profiles of the Dianchi Basin in China, Chen et al. [31] analyzed the differences between the reflectance of natural surfaces and dark and bright ISAs. The ratio of SWIR bands was selected to differentiate between natural surfaces and ISAs and MNDWI as an inhibitor. $\alpha$ (Equation (13)) is a normalizing coefficient, included to obtain an index ranging between 31 and 1. ENDISI can be computed as follows:

$$\text{ENDISI} = \frac{\text{BLUE} - \alpha \times \left( \frac{\text{SWIR1}}{\text{SWIR2}} + (\text{MNDWI})^2 \right)}{\text{BLUE} + \alpha \times \left( \frac{\text{SWIR1}}{\text{SWIR2}} + (\text{MNDWI})^2 \right)} \quad (11)$$

$$\text{MNDWI} = \frac{\text{GREEN} - \text{SWIR1}}{\text{GREEN} + \text{SWIR1}} \quad (12)$$

$$\alpha = \frac{2 \times (\text{BLUE})_{\text{Mean}}}{\left(\frac{\text{SWIR1}}{\text{SWIR2}}\right)_{\text{Mean}} + ((\text{MNDWI})^2)_{\text{Mean}}} \quad (13)$$

### 2.2.3. Thresholding Methods

In this section, we provide details about the thresholding process and an overview of the methods we utilized.

In the context of SIs, higher pixel values are considered to be the foreground and are, in our case, classified as BUAs as long as they exceed the computed threshold $T$. As suggested by Sezgin and Sankur [59], thresholding methods can be better represented under unified notations, knowing that they usually employ the statistical properties of the distribution of pixel values in an image. Pixel values of index-based maps are noted here as $g$ such as $g \in [l, L]$, with $l$ and $L$ being, respectively, the minimum and maximum values of each index-based map.

Following the 8-bit shading resolution adopted by the ImageJ software, algorithms were based on 256 bin histograms to compute the thresholds. However, these algorithms were modified to use bin centers instead of gray pixel intensities in order to accommodate unique ranges of each spectral index. Probability mass functions of the foreground and background are expressed as: $p_f(g)$ with $l \le g \le T$ and $p_b(g)$ with $T' \le g \le L$, $T'$ represents the following bin center to the value of the threshold $T$. Consequently, cumulative probability functions $P_f$ and $P_b$ (Equation (14)), means $M_f$ and $M_b$ (Equation (15)), variances $\sigma_f^2$ and $\sigma_b^2$ (Equation (16)) and Shannon entropies $H_f$ and $H_b$ (Equation (17)) of the foreground and background, respectively, all as functions of $T$ are then derived:

$$P_f = \sum_{g=l}^{T} p(g), \quad P_b = \sum_{g=T'}^{L} p(g) \quad (14)$$

$$M_f = \sum_{g=l}^{T} gp(g), \quad M_b = \sum_{g=T'}^{L} gp(g) \quad (15)$$

$$\sigma_f^2 = \sum_{g=l}^{T} [g - M_f]^2 p(g), \quad \sigma_b^2 = \sum_{g=T'}^{L} [g - M_b]^2 p(g) \quad (16)$$
$$H_f = - \sum_{g=l}^{T} p_f(g) \log p_f(g), H_b = - \sum_{g=T'}^{L} p_b(g) \log p_b(g) \quad (17)$$

An important concept of thresholding used in some of these techniques is the fuzzy set theory [60]. In image processing, fuzziness is defined as the measure of a pixel belonging to both the foreground and the background at a certain degree. It can be obtained by applying membership functions [61]. These measures are denoted as $\mu_f$ for the degree to which a pixel value $g$ belongs to the foreground and conversely $\mu_b$ to the background.

### Mean thresholding method

The mean method [62] simply takes the mean of gray pixels as a threshold (18).

$$T = \frac{l+L}{2} \quad (18)$$

### Percentile thresholding method

Percentile thresholding method [63] is a simple but effective thresholding method that involves the selection of a percentile *ptile* of pixel values that are considered to be background. This information is extracted by constructing a cumulative probability function of the image. The value was set to 70% for all images in this study. This can be represented as follows:

$$T = \arg \min_{g} \left\{ \left| \sum_{g=l}^{g} p(g) - ptile \times \sum_{g=l}^{L} p(g) \right| \right\} \quad (19)$$

### Otsu's thresholding method

Otsu's method was introduced by Otsu et al. [64] and is a popular histogram-based method in image segmentation. It aims to find the threshold that maximizes the inter-class variance, identified as the variance between the foreground and the background. It produces better results when two distinct peaks in the histogram are apparent. It can be computed using the following function:

$$T = \arg \max_{T} \left\{ \frac{P(T)[1-P(T)][M_f - M_b(T)]^2}{P(T)\sigma_f^2(T) + [1-P(T)]\sigma_b^2} \right\} \quad (20)$$

### Minimum Error thresholding method

The Minimum Error Thresholding method [65] is based on the minimization of the probability of misclassifying pixels into the foreground and the background. It assumes that pixel values of each of the two classes follow a Gaussian distribution, and it calculates the threshold that minimizes the total classification error (Equation (21)).

$$T = \arg \min_{T} \left\{ P(T) \log \sigma_f(T) + [1-P(T)] \log \sigma_b(T) - P(T) \log P(T) - [1-P(T)] \log [1-P(T)] \right\} \quad (21)$$

$\sigma_f(T)$ and $\sigma_b(T)$ are foreground and background standard deviations.

### K-means clustering method

The K-means clustering algorithm is a popular technique used in multiple fields, including image processing [66,67]. It consists of clustering pixel values into distinct groups and selecting a threshold based on the properties of these groups. It generally works by first choosing a number of clusters $k$, in our case $k = 2$, and then centroids representing each cluster are initialized. Each pixel is assigned to the nearest centroid based on the distance metric (Euclidean distance). The centroids are updated by computing the
mean of pixel values assigned to each cluster in an iterative process until the centroids exhibit no significant change. In the case of thresholding, the threshold is chosen as the mean of the two centroids.

### IsoData clustering method

IsoData [68] is an iterative procedure that starts with an initial threshold and iteratively computes the empirical averages of pixel values below and above the current threshold $T_n$, until a convergence check is reached. This check ensures that the threshold value is stabilized and no longer undergoes significant change. It can be expressed as follows:

$$T = \lim_{x \to \infty} \frac{\mu_f(T_n) + \mu_b(T_n)}{2} \qquad (22)$$

$\mu_f(T_n)$ and $\mu_b(T_n)$ are foreground and background empirical averages for the new threshold in each iteration.

### Intermodes thresholding method

Introduced by Prewitt and Mendelsohn [69], Intermodes is a method that assumes a smoothed bimodal histogram. A mean filter of size three was used iteratively until only two local maxima $j$ and $k$ were found. The threshold is then chosen as $T = j + k/2$.

### Minimum thresholding method

Another thresholding method developed by Prewitt and Mendelsohn [69] that also assumes iteratively smoothed histogram is Prewitt and Mendelsohn's minimum thresholding. The difference between Intermodes and Minimum is that the minimum point between the two resulting peaks is taken as the threshold instead of the midpoint.

### Moments thresholding method

Moment-preserving thresholding is a method introduced by Tsai [70]. It uses the statistical moments to preserve the shape characteristics of the histogram. This preservation ensures that key characteristics of the image are maintained. Up to the third moment (skewness) is considered. Gray level moments and binary image moments of order $k$ are denoted $m_k$ and $b_k$ and are calculated using Equation (24). The value of $T$ is determined by solving the system of moment preservation equations. They can be represented by Equation (23).

$$m_k = \sum_{g=l}^{L} p(g)g^k, b_k = P_f m_k^k + P_b m_b^k \qquad (23)$$

$$T = \arg\arg\min\{m_1 = b_1(T), m_2 = b_2(T), m_3 = b_3(T)\} \qquad (24)$$

$m_k^k$ and $m_b^k$ are the $k$-th moments of the foreground and the background regions, respectively.

### Triangle thresholding method

The triangle thresholding method is a histogram-based method suitable for images with one dominant peak histogram [71]. The basic idea of this method is that a line is drawn starting from the peak of the histogram to the last bin that contains pixel values, then a perpendicular line towards the histogram is considered. The chosen threshold is the value of pixel intensity that maximizes the distance between the line and the histogram. This process can be expressed as follows:

$$T = \arg\max_{n} \left\{ n \in [x_{\min}, x_{\max}], \frac{|an + bh(n) + c|}{\sqrt{n_x^2 + n_y^2}} \right\} \qquad (25)$$
with

$$ 
\begin{aligned}
x_{min} &= \arg\min_{g} \left\{ \sum_{g=l}^{L} p(g) \right\}, y_{min} = h(x_{min}) \\
x_{max} &= \arg\max_{g} \left\{ \sum_{g=l}^{L} p(g) \right\}, y_{max} = h(x_{max})
\end{aligned} 
\qquad (26) $$

## Rényi's entropy and Shaboo's thresholding method

Rényi's entropy, introduced by Rényi [72], is a generalization of Shannon's entropy and is used in many fields to measure diversity, uncertainty, or randomness in a system. For a Rényi power $\rho$, it is defined in Equation (27). The Maximum entropy maximizes the sum of Rényi's entropies for $\rho \to 1$ and Yen's thresholding method for $\rho = 2$.

$$ 
\begin{aligned}
H_f^\rho(T) &= \frac{1}{1-\rho} \ln \left( \sum_{g=l}^{T} \left[ \frac{p(g)}{P(T)} \right]^\rho \right) \\
H_b^\rho(T) &= \frac{1}{1-\rho} \ln \left( \sum_{g=T'}^{L} \left[ \frac{p(g)}{1-P(T)} \right]^\rho \right)
\end{aligned} 
\qquad (27) $$

Sahoo et al. [73] used a combination of three Rényi's entropy-based thresholds, namely $T_{[1]}$, $T_{[2]}$ and $T_{[3]}$ for $\rho \to 1$ (Maximum), $0 < \rho < 1$ and $\rho > 1$ (Yen), respectively. Values considered in our case are $\rho = 0.5$, $\rho \to 1$ and $\rho = 2$. The threshold value proposed in this method is as follows:

$$ T = T_{[1]} \cdot [P_{T_{[1]}} + \frac{1}{4} \cdot w \cdot B_1] + \frac{1}{4} \cdot T_{[2]} \cdot w \cdot B_2 + T_{[3]} \cdot [1 - P_{T_{[3]}} + \frac{1}{4} \cdot w \cdot B_3] \quad (28) $$

with

$$ P[T_{[k]}] = \sum_{g=l}^{T_{[k]}} p(g), k = 1, 2, 3, w = P[T_{[3]}] - P[T_{[1]}] \qquad (29) $$

$$ B_1, B_2, B_3 = \begin{cases} (1, 2, 1) & \text{if } |T_{[1]} - T_{[2]}| \le 5 \text{ and } |T_{[2]} - T_{[3]}| \le 5 \\ & \text{or } |T_{[1]} - T_{[2]}| > 5 \text{ and } |T_{[2]} - T_{[3]}| > 5 \\ (0, 1, 3) & \text{if } |T_{[1]} - T_{[2]}| \le 5 \text{ and } |T_{[2]} - T_{[3]}| > 5 \\ (3, 1, 0) & \text{if } |T_{[1]} - T_{[2]}| > 5 \text{ and } |T_{[2]} - T_{[3]}| \le 5 \end{cases} \qquad (30) $$

## Maximum entropy

The Maximum entropy for thresholding was introduced by Kapur et al. [74]. It is based on the consideration of an image as two different signals representing the foreground and background. The optimal threshold is considered to be the threshold that maximizes the sum of the two entropies. This can be represented as follows:

$$ T = \arg\max_{T} \{H_{f}(T) + H_{b}(T)\} \qquad (31) $$

with

$$ H_{f}(T) = - \sum_{g=l}^{T} \frac{p(g)}{P(T)} \log \frac{p(g)}{P(T)} \qquad (32) $$

$$ H_{b}(T) = - \sum_{g=T'}^{L} \frac{p(g)}{P(T)} \log \frac{p(g)}{P(T)} $$
Yen's thresholding method

Following the ideas of Pun [75] and Kapur et al. [74], Yen et al. [76] defined the entropic correlation as the sum of correlations $C_f(T)$ and $C_b(T)$ defined in Equation (34). The threshold is the function that maximizes the total entropic correlation.

$$T = \arg \max_{T} \{C_{f}(T) + C_{b}(T)\} \qquad (33)$$

with

$$C_f(T) = -\log \left( \sum_{g=l}^{T} \left[ \frac{p(g)}{P(T)} \right]^2 \right) \qquad (34)$$

$$C_b(T) = -\log \left( \sum_{g=T'}^{L} \left[ \frac{p(g)}{1-P(T)} \right]^2 \right)$$

Huang's thresholding method

This method was proposed by Huang and Wang [77] and leverages the measure of fuzziness $\mu_f[i, j]$ (Equation (36)) representing the measure for a pixel in position $(i, j)$ to belong to the foreground, where $0 \le \mu_f[i, j] \le 1$. An index of fuzziness is derived by applying the Shannon entropy, given the fuzzy measure for each pixel in the image. The threshold is then defined by minimizing the fuzziness index (Equation (35)).

$$T = \arg \min \left\{ -\frac{1}{N^2 \log 2} \sum_{g=l}^{L} \left[ \mu_f(g, T) \log(\mu_f(g, T)) + (1 - \mu_f(g, T)) \log(1 - \mu_f(g, T)) \right] p(g) \right\} \quad (35)$$

with

$$\mu_f(I[i, j], T) = \frac{L}{L + |I[i, j] - M_f(T)|} \qquad (36)$$

Shanbhag's thresholding method

This method is a modification of the Maximum entropy thresholding method, where Shanbhag [78] suggested the use of measures of fuzziness, such as the farther a gray value is from the threshold, the higher potential it obtains to belong to a specific class. The optimal threshold is found by minimizing the sum of the fuzzy entropies (Equation (38))

$$T = \arg \min_{T} \{|H_{f}(T) + H_{b}(T)|\} \qquad (37)$$

with

$$H_f(T) = - \sum_{g=l}^{T} \frac{p(g)}{P(T)} \log(\mu_f) \qquad (38)$$

$$H_b(T) = - \sum_{g=T'}^{L} \frac{p(g)}{P(T)} \log(\mu_b)$$

Li's thresholding method

Li and Tam [79] proposed a method based on the minimization of the measure of Kullback-Leibler distance between the observed image and the reconstructed image. This measure is minimized, such as the observed image and reconstructed image have the same average values in the foreground and background. It can be computed as follows:

$$T = \arg \min_{T} \left\{ \sum_{g=l}^{T} g p(g) \log \frac{g}{M_f(T)} + \sum_{g=T'}^{L} g p(g) \log \frac{g}{M_b(T)} \right\} \qquad (39)$$
with

$$ \sum_{g \le T} g = \sum_{g \le T} M_f(T), \quad \sum_{g \ge T} g = \sum_{g \ge T} M_f(T) \qquad (40) $$

### 2.2.4. Separability Analysis

The Spectral Discrimination Index (SDI) was chosen to assess the separability between built-up and other land covers. SDI is a quantitative measure used to evaluate the effectiveness of spectral bands or indices in separating LULC classes. It can be calculated using the following formula:

$$ \text{SDI} = \frac{|\mu_1 - \mu_2|}{\sigma_1 + \sigma_2} \qquad (41) $$

$\mu_1$ and $\mu_2$ are the mean spectral values of the two classes being compared, and $\sigma_1$ and $\sigma_2$ are the standard deviations of spectral values for the two classes. We also constructed histograms of classified pixel values for the visual interpretation of overlapping.

### 2.2.5. Accuracy Assessment

The assessment of SIs after the thresholding process was conducted on the basis of the Overall Accuracy (OA) and the Kappa coefficient [80]. LULC reference points were reclassified as built-up and non-built-up. All these operations were performed by flattening the resulting masks and constructing a DataFrame structure to which the mentioned metrics were applied.

## 3. Results

### 3.1. Index-Based Maps

A total of 18 index-based maps were constructed (Figures 5 and 6) as a result of the computation process of SIs. At first glance, it is noticed that NDBI in summer and SWIRED in both spring and summer resulted in altered pixel values for BUAs, meaning that they are indicated by lower values. This can be explained by analyzing the spectral profiles (Figure 4) alongside the formulation of these indices. In comparison with agriculture and bare land, the band difference of the numerators in these indices yielded lower values for BUAs due to the contrast between the bands. Therefore, to preserve the integrity of all maps for comparison purposes, they were subjected to inversion (negative transformation) by subtracting each pixel value from the corresponding maximum value of each of the maps. These maps also visually illustrate how BUAs are spatially distributed and to what extent an index is sensitive to the spectral properties of these areas.

In spring, most maps show greater contrast between BUAs and vegetation (Taking the SVM classification map in Figure 2 as a visual reference to the spatial distribution of land covers) due to the enhanced spectral distinctions from increased vegetation and moisture, at the same time, some indices show similar pixel values of bare land, trees and/or water to BUAs, such as NDBI, ENDISI, and PISI. The contrast between land covers in summer becomes generally less pronounced in the majority of indices, and more confusion between land covers is noticed in NBAI, VRNIR, PISI, and ENDISI.

### 3.2. Separability Potential of SIs and Sensitivity to Seasonal Variations

The chosen SIs leverage different design principles to capture BUAs. Their ability to distinguish BUAs from other land cover types, evaluated via SDI (Table 3) and visual histogram overlap (Figure 7), demonstrated differences in behavior across land covers and seasons. Even though it can be seen that all SIs performed well for distinguishing BUAs from agriculture during spring, NDBI, VRNIR, ENDISI, PISI, and NBAI showed a decline in SDI values and, therefore, moderate to poor separability during summer, showcasing their sensitivity to seasonal variation. Similarly, the same indices also show
poor to moderate separability regarding bare land in both seasons. This decline of separability potential introduces the need for seasonally adaptive indices for BUAs mapping in semi-arid climates.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Index-based maps in spring, plotted with a 2% standard deviation stretching to enhance more subtlety differences.</figcaption></figure>

SWIRED, BUI, and BLFEI, on the other hand, demonstrated greater efficiency, maintaining SDI values above 1 in both seasons for the different land covers except water. The lowest SDI values and overlap in these indices were yielded for BUAs vs. bare land in summer (1.08 and 1.05). They also show the highest SDI values (1.23, 1.22, and 1.42) for BUAs vs. agriculture in the same season compared to other indices, making them the most reliable for summer scenes. With respect to water, a final separability assessment cannot be concluded for most SIs, as the number of water pixels is very small. Thus, a water masking process is recommended in areas where water bodies are heavily featured.

<table><caption>**Table 3.** SDI values of BUAs and each of the other land covers and for each spectral index in summer and spring. Here, $SDI > 1$ indicates good separability, $0.5 < SDI < 1$ indicates moderate separability, $SDI < 0.5$ indicates poor separability.</caption><thead><tr><th rowspan="2">Index</th><th colspan="2">Built-Up/Agriculture</th><th colspan="2">Built-Up/Bare Land</th><th colspan="2">Built-Up/Trees</th><th colspan="2">Built-Up/Water</th></tr><tr><th>Spring</th><th>Summer</th><th>Spring</th><th>Summer</th><th>Spring</th><th>Summer</th><th>Spring</th><th>Summer</th></tr></thead><tbody><tr><td>NDBI</td><td>1.57</td><td>0.67</td><td>0.34</td><td>0.83</td><td>0.65</td><td>0.00</td><td>0.62</td><td>0.57</td></tr><tr><td>VRNIR</td><td>2.15</td><td>0.56</td><td>0.89</td><td>0.51</td><td>1.67</td><td>1.16</td><td>0.51</td><td>0.01</td></tr><tr><td>SWIRED</td><td>2.14</td><td>1.23</td><td>1.22</td><td>1.08</td><td>2.03</td><td>1.67</td><td>0.19</td><td>0.28</td></tr><tr><td>BUI</td><td>2.10</td><td>1.22</td><td>1.21</td><td>1.08</td><td>2.02</td><td>1.66</td><td>0.19</td><td>0.27</td></tr><tr><td>ENDISI</td><td>1.72</td><td>0.52</td><td>0.45</td><td>0.31</td><td>1.25</td><td>1.25</td><td>1.41</td><td>2.14</td></tr><tr><td>PISI</td><td>1.34</td><td>0.79</td><td>0.93</td><td>0.76</td><td>0.36</td><td>0.54</td><td>1.00</td><td>1.01</td></tr><tr><td>CBCI</td><td>2.02</td><td>0.92</td><td>1.01</td><td>0.89</td><td>1.19</td><td>0.90</td><td>0.38</td><td>0.75</td></tr><tr><td>BLFEI</td><td>2.28</td><td>1.42</td><td>1.21</td><td>1.05</td><td>1.97</td><td>1.67</td><td>0.06</td><td>0.34</td></tr><tr><td>NBAI</td><td>1.76</td><td>0.59</td><td>0.46</td><td>0.29</td><td>1.66</td><td>1.30</td><td>2.44</td><td>2.39</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Index-based maps in summer, plotted with a 2% standard deviation stretching to enhance more subtlety differences.</figcaption></figure>

### 3.3. Characteristics of Pixel Intensity Distributions

Histogram shapes and the distribution of pixel values are a tipping point in the performance of index-based maps and are critical in the thresholding process. The resulting overall intensity distributions in Figure 7a reveal that spring maps generally show bimodal or slightly multimodal characteristics with sufficient range between the peaks, as well as most distributions being right-skewed and asymmetric, due to agriculture and bare land combined having a larger proportion than BUAs in the study area. In summer, however (Figure 7b), distributions become more compressed with high peaks and small ranges. Standard deviations, skewness, and kurtosis (flatness) can serve as valuable metrics to evaluate the variations in histogram properties [81]. They are shown in Table 4.

<table><caption>**Table 4.** Standard deviation, skewness, and kurtosis of pixel values in index-based maps.</caption><thead><tr><th rowspan="2">Index</th><th colspan="2">Standard Deviation</th><th colspan="2">Skewness</th><th colspan="2">Kurtosis</th></tr><tr><th>Spring</th><th>Summer</th><th>Spring</th><th>Summer</th><th>Spring</th><th>Summer</th></tr></thead><tbody><tr><td>NDBI</td><td>0.18</td><td>0.06</td><td>0.11</td><td>1.40</td><td>-1.25</td><td>6.05</td></tr><tr><td>VRNIR</td><td>0.23</td><td>0.07</td><td>0.57</td><td>-2.18</td><td>-0.93</td><td>12.58</td></tr><tr><td>SWIRED</td><td>0.15</td><td>0.07</td><td>0.51</td><td>0.48</td><td>-0.89</td><td>0.54</td></tr><tr><td>BUI</td><td>0.17</td><td>0.08</td><td>0.36</td><td>0.38</td><td>-0.99</td><td>0.71</td></tr><tr><td>ENDISI</td><td>0.27</td><td>0.12</td><td>0.54</td><td>0.44</td><td>-0.66</td><td>2.38</td></tr><tr><td>PISI</td><td>0.05</td><td>0.02</td><td>0.31</td><td>-0.01</td><td>-0.41</td><td>5.25</td></tr><tr><td>CBCI</td><td>0.13</td><td>0.04</td><td>0.55</td><td>-0.02</td><td>-0.80</td><td>3.84</td></tr><tr><td>BLFEI</td><td>0.10</td><td>0.05</td><td>0.53</td><td>0.87</td><td>-0.84</td><td>1.32</td></tr><tr><td>NBAI</td><td>0.07</td><td>0.05</td><td>1.09</td><td>1.09</td><td>1.46</td><td>5.32</td></tr></tbody></table>
As a consequence, the adaptability of thresholding methods can be seasonally driven, as a method may not be suitable for the same index in both seasons.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Histograms of pixel values and their respective composition of land cover distributions according to the SVM classification, grouped in 30 bins in (a) Spring; (b) Summer.</figcaption></figure>

## 3.4. Accuracy Assessment

### 3.4.1. Thresholding Consistency

Consequent to the thresholding process, a total of 288 binary maps were generated and assessed by calculating their OAs and Kappa coefficients. The radar plot displayed in Figure 8 discloses that some thresholding methods can be consistent in both seasons and can extract reliable thresholds for indices with stable or relatively stable seasonal behavior, being SWIRED, BUI, BLFEI, CBCI, and PISI. These methods are K-means, Huang, IsoData, Li, Otsu, Percentile, Shanbhag, and Moments, resulting in mean accuracies between 78% and 85% and mean Kappa values between 55% and 65% across all indices. Intermodes and entropy-based methods (Max Entropy, Renyi's entropy, and Yen) generally struggled with spectral indices producing low contrast (standard deviations between 0.02 and 0.10) such as NBAI and PISI in spring and most if not all indices during summer.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Radar plots of OAs and Kappa coefficients corresponding to thresholding methods and each of the SIs in: (a) Spring, (b) Summer.</figcaption></figure>

### 3.4.2. Overview of Highest Performing Binary Maps and Manual Thresholding

In Table 5, a detailed summary of thresholding performance metrics for the best-performing binary masks of each SI in spring and summer is presented. Metrics include OAs and Kappa coefficients, comparing outcomes from the most accurate automatic thresholding methods (ATM) with results from the manual thresholding process for each index and season. Thresholds considered in automatic and manual methods are the ones that maximize Kappa values.

In spring, BLFEI, SWIRED, and BUI achieved the highest accuracies with OAs of 93.97%, 93.39%, and 92.81%, respectively, with the Minimum method defining close thresholds to the thresholds found by manual thresholding. VRNIR and CBCI also showed strong performance by yielding OAs of 88.35% and 88.99% and Kappa values exceeding
70%. Percentile method (ptile value of 0.7), in this case, extracted thresholds with more offset to manual thresholding but without significantly affecting accuracy (1.24% and 1.57% OA differences in favor of the manual thresholding). This method was the most effective for NDBI and NBAI as well. It is worth pointing out that Percentile can be considered a semi-automatic method. It can compute accurate thresholds as long as an estimated percentage of the proportion of non-built-up areas in the data is passed.

In summer, however, results vary due to changes in histogram characteristics. BLFEI, SWIRED, and BUI maintain OAs above or little below 90% and Kappa scores between 73% and 77%, by means of Huang, Shanbhag, and Li, which correlates with their adaptability to bimodal, narrower peaked distributions. In the case of VRNIR, which shows a very high kurtosis (high peak) but has a shallow bimodal histogram, manual thresholding resulted in a 9% higher Kappa value over Percentile, indicating that global thresholding can be limited in very high peaked distributions. On the other hand, indices with unimodal distributions like NBAI and ENDISI showed poor performances for both manual and automatic thresholding, suggesting that these indices are not applicable for separating BUAs from other land covers in summer.

<table><caption>**Table 5.** Thresholds, Overall Accuracies, Kappa coefficients and used thresholding methods resulting in the best-performing binary masks in terms of maximum Kappa values for spring and summer, in comparison to manual thresholding.</caption><thead><tr><th rowspan="2">Index</th><th rowspan="2">Season</th><th colspan="2">Threshold</th><th colspan="2">OA (%)</th><th colspan="2">Kappa (%)</th><th rowspan="2">ATM</th></tr><tr><th>ATM</th><th>Manual</th><th>ATM</th><th>Manual</th><th>ATM</th><th>Manual</th></tr></thead><tbody><tr><td rowspan="2">NDBI</td><td>Spring</td><td>-0.11215</td><td>-0.16703</td><td>81.49</td><td>80.57</td><td>53.91</td><td>56.15</td><td>Percentile</td></tr><tr><td>Summer</td><td>0.20855</td><td>0.21272</td><td>83.30</td><td>84.13</td><td>60.00</td><td>60.10</td><td>Percentile</td></tr><tr><td rowspan="2">VRNIR</td><td>Spring</td><td>-0.47940</td><td>-0.38573</td><td>88.01</td><td>89.50</td><td>71.09</td><td>71.91</td><td>Shanbhag</td></tr><tr><td>Summer</td><td>-0.17484</td><td>-0.16076</td><td>80.90</td><td>86.19</td><td>52.32</td><td>61.44</td><td>Percentile</td></tr><tr><td rowspan="2">SWIRED</td><td>Spring</td><td>0.49814</td><td>0.49438</td><td>93.39</td><td>93.47</td><td>82.80</td><td>83.09</td><td>Minimum</td></tr><tr><td>Summer</td><td>0.44645</td><td>0.45980</td><td>89.17</td><td>90.41</td><td>73.02</td><td>75.01</td><td>Huang & Shanbhag</td></tr><tr><td rowspan="2">BUI</td><td>Spring</td><td>-0.34257</td><td>-0.32459</td><td>92.81</td><td>93.55</td><td>81.69</td><td>83.23</td><td>Minimum</td></tr><tr><td>Summer</td><td>-0.26458</td><td>-0.25523</td><td>89.67</td><td>90.58</td><td>73.90</td><td>75.49</td><td>Huang</td></tr><tr><td rowspan="2">ENDISI</td><td>Spring</td><td>-0.19280</td><td>-0.19998</td><td>85.78</td><td>85.78</td><td>65.42</td><td>65.68</td><td>Max Entropy</td></tr><tr><td>Summer</td><td>-0.28870</td><td>-0.25905</td><td>75.70</td><td>78.92</td><td>39.26</td><td>41.00</td><td>Percentile</td></tr><tr><td rowspan="2">PISI</td><td>Spring</td><td>-0.05830</td><td>-0.05483</td><td>87.77</td><td>88.18</td><td>68.88</td><td>69.28</td><td>Renyi's Entropy</td></tr><tr><td>Summer</td><td>-0.01663</td><td>-0.01462</td><td>83.88</td><td>84.38</td><td>59.03</td><td>59.40</td><td>Huang & Shanbhag</td></tr><tr><td rowspan="2">CBCI</td><td>Spring</td><td>0.13308</td><td>0.15883</td><td>88.68</td><td>89.91</td><td>72.27</td><td>73.71</td><td>Percentile</td></tr><tr><td>Summer</td><td>0.25985</td><td>0.26006</td><td>86.77</td><td>86.69</td><td>65.41</td><td>65.66</td><td>IsoData</td></tr><tr><td rowspan="2">BLFEI</td><td>Spring</td><td>-0.26332</td><td>-0.26856</td><td>93.97</td><td>94.13</td><td>84.29</td><td>84.84</td><td>Minimum</td></tr><tr><td>Summer</td><td>-0.23654</td><td>-0.22596</td><td>91.57</td><td>91.74</td><td>77.77</td><td>78.06</td><td>Li</td></tr><tr><td rowspan="2">NBAI</td><td>Spring</td><td>-0.84534</td><td>-0.84341</td><td>86.94</td><td>87.60</td><td>67.74</td><td>68.95</td><td>Percentile</td></tr><tr><td>Summer</td><td>-0.75589</td><td>-0.75377</td><td>79.25</td><td>80.50</td><td>44.16</td><td>45.18</td><td>Moments</td></tr></tbody></table>

A visual representation of a sample area of best-performing binary maps is shown in Figures 9 and 10. All SIs during spring had a high identification performance of BUAs and generally had a low misclassification output for other land covers except NDBI, NBAI, and ENDISI, misclassifying bare land and hence overestimating BUAs. With regard to trees, both NDBI and PISI misclassified a small portion of the chosen woodlot (area n°4) and other tree pixels from the northern-west part of the sample map. Additionally, an underestimation of BUAs is the most apparent among indices like VRNIR, NBAI, and CBCI
during summer, and an overestimation of BUAs was observed for trees in NDBI and bare land in NBAI and ENDISI.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Best performing binary masks for each index in a sample area during spring and Google Earth imagery from February 2021. Refer to Table 5 for accuracies and thresholding methods. Red marks represent a sample of: 1. BUAs, 2. Bare land, 3. Agriculture, 4. Trees and 5. Water.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Best performing binary masks for each index in a sample area during summer and Google Earth imagery from November 2021. Refer to Table 5 for accuracies and thresholding methods. Red marks represent a sample of: 1. BUAs, 2. Bare land, 3. Agriculture, 4. Trees and 5. Water.</figcaption></figure>
4. Discussion

This paper suggests an evaluation of the performance of nine SIs (NDBI, VRNIR, SWIRED, NBAI, BLFEI, BUI, CBCI, PISI, and ENDISI) and offers an examination of their sensitivity to variations between spring and summer in a semi-arid climate, based on Landsat 8 imagery. Additionally, sixteen thresholding methods were applied to test which are the most suitable and consistent in accordance with the numerical distribution of pixel values in each of the resulting index-based maps. BLFEI, SWIRED, and BUI resulted in the highest consistency of separation between land covers in both seasons and, therefore, the highest accuracy, while other indices also showed good accuracies but were more sensitive to seasonal variations. Moreover, multiple thresholding methods proved consistency in setting appropriate thresholds when bimodal distributions are manifested in both seasons, including K-means, Huang, IsoData, Li, Moments, Otsu, Percentile, and Shanbhag.

The quantitative and visual separability analysis reveals that indices distinguish
between built-up and other land covers differently. Accuracy assessment indicates that
the most affected indices by seasonal variation are VRNIR, ENDISI, and NBAI. VRNIR
uses the NIR and Red bands to promote the contrast between agriculture/trees, which
are characterized by high phenology and moisture in spring, and BUAs, making VRNIR
effective during this season. The same can be said for bare land, but with less contrast in
the visible bands. This contrast is significantly reduced during summer in these bands due
to the dryness of the soil, as the three land covers exhibit similar reflectance, which explains
its poor performance in summer. Although SWIR bands mark a sharp increase in both
reflectance and contrast during summer, the reliance on ratios involving these bands as
normalizers in NBAI and ENDISI is susceptible of leading to BUAs having similar values
to agriculture and bare land in summer, considering the decrease of reflectance in SWIR2
(Figure 3b). Under these circumstances, some indices were less affected by variations but
still show some degree of confusion between land covers, such as PISI and CBCI in summer
and NDBI in spring.

Furthermore, it is revealed by looking into the literature that statements with regard
to the optimal season for extracting BUAs from satellite imagery differ greatly [40–42,44].
However, our findings correlate with the ideas confirmed by Xi et al. [44], Capolupo et al. [57]
and Bouhennache et al. [32], stating that the use of a combination of visible and SWIR
bands, such as in the form of BLFEI, SWIRED and BUI can be more effective. It is further
proved by the present study that they can be used in both wet and dry conditions to
obtain high accuracies in regions with similar geographical and reflectance characteristics
to those of our study area. With respect to thresholding methods, the most used methods
in the literature are the manual, Otsu, and IsoData. It has been proved in multiple studies
that Otsu's method is robust when the distribution of pixel values shows a clear bimodal
histogram [26,32]. Li et al. [40] showed that IsoData was better than Otsu at extracting good
thresholds when histograms lean towards multimodal or unimodal distributions. This
study, however, suggests that Otsu's method can also yield good results for multimodal
distributions such as in NDBI, CBCI, and ENDISI in spring. Nonetheless, such distributions
lead to considerable confusion with bare land. This study also adds to the conversation
that the Minimum method can be powerful in finding the closest thresholds to those
found by manual thresholding for indices with high separation potential in spring (BLFEI,
SWIRED and BUI). Günen and Atasever [48] also found this particular method to be the
most effective for water bodies mapping. This can be explained by the fact that BUAs on a
city scale can share, in most cases, bimodal characteristics in histograms with lakes at their
extent scale. On the other hand, Huang, Shanbhag, Li, and Percentile were found to yield
the highest accuracies in more peaked bimodal distributions.
The use of adaptive SIs combined with reliable thresholding methods can offer an alternative or complementary approach to classification algorithms and machine learning-based methods by providing a practical, scalable, data-efficient, and cost-effective outcome that reduces the reliance on labeled datasets and field surveys. Incorporating adaptive SIs into automated remote sensing workflows can serve urban planners and policymakers to track urban expansion more efficiently and dynamically, particularly under more and more intensified effects of climate change, as per the alteration of rainfall patterns that leads to more variability in soil moisture and vegetation phenology. Alternatively, this study can serve as a reference for supporting the choice of SIs in regions with similar reflectance and climate properties to those of our study area for the purpose of mapping BUAs. It can also be informative in other frameworks, such as in other LULC applications or in the context of automatic image thresholding processes.

**5. Limitations and Future Research**

In terms of consistency, it is safe to say that SIs are generally better used when moisture and phenology of soils are high in semi-arid regions or regions with similar reflectance responses to those of our study area. Considering the latter statement, limitations of this study arise from the inability to project these findings on other geographical properties, such as topography, soil composition, or coastal proximity. On the other hand, the influence of land cover proportions with reference to the separability potential of SIs was not explicitly analyzed in this study, as histogram shapes and thresholding performance can be affected by these proportions and the size of the area of interest. Hence, these elements should be considered in future research for both existing and upcoming design processes of SIs. Furthermore, the use of Landsat imagery to map BUAs can be beneficial for its long-term availability and suitability to most SIs in the literature. Nevertheless, similar to Xi et al. [44], it has been observed that its moderate spatial resolution resulted in finer urban features such as peri-urban areas and thinner roads not being detected. The use of higher-resolution multi-spectral imagery, such as Sentinel-2 imagery, can be more practical. Thus, the design of new SIs and more testing of existing SIs on this type of imagery should be explored in more depth. Lastly, this study focused primarily on spectral separability and did not assess the impact of urban morphology, such as building density and material composition, on the performance of SIs. Thus, future research could analyze how different urban structures influence spectral reflectance patterns in medium-resolution imagery to further refine index-based methods.

The findings of this study are meant to be used to conduct a long-term, detailed remote sensing-based assessment of the Surface Urban Heat Island (SUHI) effect at the city scale in the different urban settings of the province of Settat, where automatically mapping BUAs in different seasons is required to extract Land Surface Temperatures (LSTs) in these areas, thus gaining more insights on this phenomenon in such regions.

**6. Conclusions**

The need for seasonally robust SIs is a cornerstone for ensuring a reliable and automated BUAs mapping process. It not only enhances the resilience of these processes to environmental changes but also supports more accurate and detailed monitoring of urban expansion and resource allocation in the face of a rapidly changing climate. On account of the assessment of nine SIs (NDBI, VRNIR, SWIRED, NBAI, BLFEI, BUI, CBCI, PISI and ENDISI) and sixteen thresholding methods (K-means, Huang, Intermodes, IsoData, Li, Maximum entropy, Mean, Minimum error, Minimum, Moments, Otsu, Percentile, Renyi's entropy, Shanbhag, Triangle, and Yen) performed in this study, the following conclusions were reached:
* The performance of an SI depends on the contrast it produces between BUAs and each of the other land covers. As a consequence, it can be deduced that the proportions of these land covers in an area of interest are also to be taken into consideration when adopting an index, as some SIs can distinguish BUAs from one or more particular land cover types better than others.

* Overall, SIs are better applied in wet conditions for semi-arid climate or regions with similar reflectance properties to the study area. Nonetheless, BLFEI, SWIRED, and BUI can achieve high accuracies in both wet and dry conditions, as they provide a balanced separability between BUAs and other land covers except water in both conditions. In contrast, ENDISI and NBAI resulted in poor separability and unimodal distributions in summer and, therefore, are deemed unsuitable in dry conditions.

* In terms of consistency, multiple thresholding methods are applicable in both dry and wet conditions when bimodal distributions are observed and can be reliable for automation processes. These methods are K-means, Huang, IsoData, Li, Moments, Otsu, Percentile, and Shanbhag. Other methods, such as Intermodes, Maximum entropy, Renyi's entropy, Yen, Minimum error, Minimum, and Triangle, are not recommended when no evidence of the distribution of pixel values is available.

* The use of BLFEI in combination with the Minimum method was found to be the most effective approach under wet conditions in semi-arid climates. Likewise, the use of BLFEI in combination with either Li, Huang, Triangle, Otsu, K-means, or IsoData is the most effective approach under dry conditions.

**Author Contributions:** Conceptualization, Y.H., A.R. and R.A.; methodology, Y.H., A.R. and R.A.; software, Y.H.; validation, Y.H., A.R. and R.A.; formal analysis, Y.H.; investigation, Y.H., A.R. and R.A.; resources, Y.H.; data curation, Y.H.; writing—original draft preparation, Y.H.; writing—review and editing, Y.H., A.R. and R.A.; visualization, Y.H.; supervision, A.R. and R.A.; project administration, Y.H., A.R. and R.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Landsat 8 satellite imagery used in this study is publicly available from the United States Geological Survey (USGS) through the EarthExplorer platform (https://earthexplorer.usgs.gov (accessed on 10 May 2024)). It was derived from resources available in the public domain under the USGS open data policy. The OpenStreetMap (OSM) data utilized in Figure 1 was obtained from the official OpenStreetMap database (https://www.openstreetmap.org (accessed on 10 April 2024)). It was also accessed from public domain resources and was provided under the Open Database License (ODbL).

**Conflicts of Interest:** The authors declare no conflicts of interest.

## References

1. Li, X.; Zhou, Y.; Hejazi, M.; Wise, M.; Vernon, C.; Iyer, G.; Chen, W. Global urban growth between 1870 and 2100 from integrated high resolution mapped data and urban dynamic modeling. *Commun. Earth Environ.* **2021**, *2*, 201. [CrossRef]

2. Schiavina, M.; Melchiorri, M.; Corbane, C.; Freire, S.; Batista e Silva, F. Built-up areas are expanding faster than population growth: Regional patterns and trajectories in Europe. *J. Land Use Sci.* **2022**, *17*, 591–608. [CrossRef]

3. Abd-Elhamid, H.F.; Zeleňáková, M.; Vranayová, Z.; Fathy, I. Evaluating the impact of urban growth on the design of storm water drainage systems. *Water* **2020**, *12*, 1572. [CrossRef]

4. Fattah, M.A.; Morshed, S.R.; Morshed, S.Y. Impacts of land use-based carbon emission pattern on surface temperature dynamics: Experience from the urban and suburban areas of Khulna, Bangladesh. *Remote Sens. Appl. Soc. Environ.* **2021**, *22*, 100508. [CrossRef]

5. Rahman, M.N.; Akter, K.S.; Faridatul, M.I. Assessing the Impact of Urban Expansion on Carbon Emission. *Environ. Sustain. Indic.* **2024**, *23*, 100416. [CrossRef]

6. Huang, X.; Liu, Y.; Stouffs, R. Human-earth system dynamics in China's land use pattern transformation amidst climate fluctuations and human activities. *Sci. Total Environ.* **2024**, *954*, 176013. [CrossRef]
7. Oyadeyi, O.A.; Oyadeyi, O.O. Towards inclusive and sustainable strategies in smart cities: A comparative analysis of Zurich, Oslo, and Copenhagen. *Res. Glob.* **2025**, *10*, 100271. [CrossRef]

8. Belaïd, F.; Apeaning, R.; Arora, A. Advancing Cities' Resilience: A Comprehensive Review and Mapping of Research. In *Climate-Resilient Cities*; Springer: Cham, Switzerland, 2025; pp. 1–15.

9. Zhao, Y.; Zhu, Z. ASI: An artificial surface Index for Landsat 8 imagery. *Int. J. Appl. Earth Obs. Geoinf.* **2022**, *107*, 102703. [CrossRef]

10. Nguyen, T.; Nguyen, N.P.; Savaglio, C.; Zhang, Y.; Dumba, B. Future of urban remote sensing and new sensors. *Eur. J. Remote Sens.* **2023**, *56*, 2281073. [CrossRef]

11. More, N.; Kanthe, A.; Gangane, J.; Padiya, P.; Sumithra, T. A Complete Study of Remote Sensing-Sentinel-2 Satellite Data for Land Use/Land Cover (LULC) Analysis. *Panam. Math. J.* **2025**, *35*, 231–249. [CrossRef]

12. Kebede, T.A.; Hailu, B.T.; Suryabhagavan, K.V. Evaluation of spectral built-up indices for impervious surface extraction using Sentinel-2A MSI imageries: A case of Addis Ababa city, Ethiopia. *Environ. Chall.* **2022**, *8*, 100568. [CrossRef]

13. Sun, G.; Chen, X.; Jia, X.; Yao, Y.; Wang, Z. Combinational build-up index (CBI) for effective impervious surface mapping in urban areas. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2015**, *9*, 2081–2092. [CrossRef]

14. Lillesand, T.; Kiefer, R.W.; Chipman, J. *Remote Sensing and Image Interpretation*; John Wiley & Sons: Hoboken, NJ, USA, 2015.

15. Soni, P.K.; Rajpal, N.; Mehta, R.; Mishra, V. Urban land cover and land use classification using multispectral sentinal-2 imagery. *Multimed. Tools Appl.* **2021**, *81*, 36853–36867. [CrossRef]

16. Blaschke, T. Object based image analysis for remote sensing. *Isprs J. Photogramm. Remote Sens.* **2010**, *65*, 2–16. [CrossRef]

17. Zhang, C.; Wang, Q.; Atkinson, P.M. Unsupervised object-based spectral unmixing for subpixel mapping. *Remote Sens. Environ.* **2025**, *318*, 114514. [CrossRef]

18. Javed, A.; Cheng, Q.; Peng, H.; Altan, O.; Li, Y.; Ara, I.; Huq, E.; Ali, Y.; Saleem, N. Review of spectral indices for urban remote sensing. *Photogramm. Eng. Remote Sens.* **2021**, *87*, 513–524. [CrossRef]

19. Kaur, R.; Pandey, P. A review on spectral indices for built-up area extraction using remote sensing technology. *Arab. J. Geosci.* **2022**, *15*, 391. [CrossRef]

20. Jozdani, S.E.; Johnson, B.; Chen, D. Comparing Deep Neural Networks, Ensemble Classifiers, and Support Vector Machine Algorithms for Object-Based Urban Land Use/Land Cover Classification. *Remote Sens.* **2019**, *11*, 1713. [CrossRef]

21. Maxwell, A.E.; Warner, T.A.; Fang, F. Implementation of machine-learning classification in remote sensing: An applied review. *Int. J. Remote Sens.* **2018**, *39*, 2784–2817. [CrossRef]

22. Youssef, R.; Aniss, M.; Jamal, C. Machine Learning and Deep Learning in Remote Sensing and Urban Application: A Systematic Review and Meta-Analysis. In *Proceedings of the 4th Edition of International Conference on Geo-IT and Water Resources 2020, Geo-IT and Water Resources 2020, Hoceima, Morocco, 11–12 March 2020*. [CrossRef]

23. Thunig, H.; Wolf, N.; Naumann, S.; Siegmund, A.; Jürgens, C.; Uysal, C.; Maktav, D. Land use/land cover classification for applied urban planning-the challenge of automation. In *Proceedings of the 2011 Joint Urban Remote Sensing Event, Munich, Germany, 11–13 April 2011*; pp. 229–232.

24. Polykretis, C.; Grillakis, M.G.; Alexakis, D.D. Exploring the impact of various spectral indices on land cover change detection using change vector analysis: A case study of Crete Island, Greece. *Remote Sens.* **2020**, *12*, 319. [CrossRef]

25. Estoque, R.C.; Murayama, Y. Classification and change detection of built-up lands from Landsat-7 ETM+ and Landsat-8 OLI/TIRS imageries: A comparative assessment of various spectral indices. *Ecol. Indic.* **2015**, *56*, 205–217. [CrossRef]

26. Tian, Y.; Chen, H.; Song, Q.; Zheng, K. A novel index for impervious surface area mapping: Development and validation. *Remote Sens.* **2018**, *10*, 1521. [CrossRef]

27. Zhang, S.; Yang, K.; Li, M.; Ma, Y.; Sun, M. Combinational biophysical composition index (CBCI) for effective mapping biophysical composition in urban areas. *IEEE Access* **2018**, *6*, 41224–41237. [CrossRef]

28. Zha, Y.; Gao, J.; Ni, S. Use of normalized difference built-up index in automatically mapping urban areas from TM imagery. *Int. J. Remote Sens.* **2003**, *24*, 583–594. [CrossRef]

29. Waqar, M.M.; Mirza, J.F.; Mumtaz, R.; Hussain, E. Development of New Indices for Extraction of Built-Up Area and Bare Soil from Landsat Data. *Open Access Sci. Rep.* **2012**, *1*, 4.

30. Kaimaris, D.; Patias, P. Identification and area measurement of the built-up area with the Built-up Index (BUI). *Int. J. Adv. Remote Sens.* **GIS** **2016**, *5*, 1844–1858. [CrossRef]

31. Chen, J.; Yang, K.; Chen, S.; Yang, C.; Zhang, S.; He, L. Enhanced normalized difference index for impervious surface area estimation at the plateau basin scale. *J. Appl. Remote Sens.* **2019**, *13*, 016502. [CrossRef]

32. Bouhennache, R.; Bouden, T.; Taleb-Ahmed, A.; Cheddad, A. A new spectral index for the extraction of built-up land features from Landsat 8 satellite imagery. *Geocarto Int.* **2019**, *34*, 1531–1551. [CrossRef]

33. As-Syakur, A.R.; Adnyana, I.W.S.; Arthana, I.W.; Nuarsa, I.W. Enhanced built-up and bareness index (EBBI) for mapping built-up and bare land in an urban area. *Remote Sens.* **2012**, *4*, 2957–2970. [CrossRef]
34. Piyoosh, A.K.; Ghosh, S.K. Development of a modified bare soil and urban index for Landsat 8 satellite data. *Geocarto Int.* **2018**, *33*, 423–442. [CrossRef]

35. Chen, Z.; Yu, B.; Zhou, Y.; Liu, H.; Yang, C.; Shi, K.; Wu, J. Mapping global urban areas from 2000 to 2012 using time-series nighttime light data and MODIS products. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2019**, *12*, 1143–1153. [CrossRef]

36. Senthilkumaran, N.; Vaithegi, S. Image segmentation by using thresholding techniques for medical images. *Comput. Sci. Eng. Int. J.* **2016**, *6*, 1–13.

37. Sekertekin, A. A survey on global thresholding methods for mapping open water body using Sentinel-2 satellite imagery and normalized difference water index. *Arch. Comput. Methods Eng.* **2021**, *28*, 1335–1347. [CrossRef]

38. Li, Y.; Xu, X.; Wu, Z.; Fan, H.; Tong, X.; Liu, J. A forest type-specific threshold method for improving forest disturbance and agent attribution mapping. *Giscience Remote Sens.* **2022**, *59*, 1624–1642. [CrossRef]

39. Heiden, U.; d'Angelo, P.; Schwind, P.; Karlshöfer, P.; Müller, R.; Zepp, S.; Wiesmeier, M.; Reinartz, P. Soil reflectance compos- ites—improved thresholding and performance evaluation. *Remote Sens.* **2022**, *14*, 4526. [CrossRef]

40. Li, C.; Shao, Z.; Zhang, L.; Huang, X.; Zhang, M. A comparative analysis of index-based methods for impervious surface mapping using multiseasonal Sentinel-2 satellite data. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2021**, *14*, 3682–3694. [CrossRef]

41. Almohamad, H.; Alshwesh, I.O. Evaluation of Index-Based Methods for Impervious Surface Mapping from Landsat-8 to Cities in Dry Climates: A Case Study of Buraydah City, KSA. *Sustainability* **2023**, *15*, 9704. [CrossRef]

42. Chen, J.; Chen, S.; Yang, C.; He, L.; Hou, M.; Shi, T. A comparative study of impervious surface extraction using Sentinel-2 imagery. *Eur. J. Remote Sens.* **2020**, *53*, 274–292. [CrossRef]

43. Sekertekin, A.; Abdikan, S.; Marangoz, A.M. The acquisition of impervious surface area from LANDSAT 8 satellite sensor data using urban indices: A comparative analysis. *Environ. Monit. Assess.* **2018**, *190*, 381. [CrossRef]

44. Xi, Y.; Thinh, N.X.; Li, C. Preliminary comparative assessment of various spectral indices for built-up land derived from Landsat-8 OLI and Sentinel-2A MSI imageries. *Eur. J. Remote Sens.* **2019**, *52*, 240–252. [CrossRef]

45. Zhang, H.; Zhang, Y.; Lin, H. Seasonal effects of impervious surface estimation in subtropical monsoon regions. *Int. J. Digit. Earth* **2014**, *7*, 746–760. [CrossRef]

46. Su, S.; Tian, J.; Dong, X.; Tian, Q.; Wang, N.; Xi, Y. An impervious surface spectral index on multispectral imagery using visible and near-infrared bands. *Remote Sens.* **2022**, *14*, 3391. [CrossRef]

47. Shao, Z.; Cheng, T.; Fu, H.; Li, D.; Huang, X. Emerging issues in mapping urban impervious surfaces using high-resolution remote sensing images. *Remote Sens.* **2023**, *15*, 2562. [CrossRef]

48. Günen, M.A.; Atasever, U.H. Remote sensing and monitoring of water resources: A comparative study of different indices and thresholding methods. *Sci. Total Environ.* **2024**, *926*, 172117. [CrossRef] [PubMed]

49. Beck, H.E.; Zimmermann, N.E.; McVicar, T.R.; Vergopolan, N.; Berg, A.; Wood, E.F. Present and future Köppen-Geiger climate classification maps at 1-km resolution. *Sci. Data* **2018**, *5*, 180214. [CrossRef] [PubMed]

50. HCP. *Monographie de la Province de Settat*; Technical Report; Haut Commissariat du Plan du Royaume du Maroc: Casablanca, Morocco, 2020.

51. Driouech, F. Distribution des précipitations hivernales sur le Maroc dans le cadre d'un changement climatique: Descente d'échelle et incertitudes. Ph.D. Thesis, Institut National Polytechnique de Toulouse-INPT, Toulouse, France, 2010.

52. DGM. *Maroc-Etat du climat en 2023*; Technical Report; Direction Général de Météorologie-Ministér de l'Equipement et de l'Eau: Casablanca, Morocco, 2024.

53. HCP. *Recessement Général de la Population et de l'Habitat 2024*; Technical Report; Haut Commissariat du Plan du Royaume du Maroc: Casablanca, Morocco, 2024.

54. Lewis, A.; Lacey, J.; Mecklenburg, S.; Ross, J.; Siqueira, A.; Killough, B.; Szantoi, Z.; Tadono, T.; Rosenavist, A.; Goryl, P.; et al. CEOS analysis ready data for Land (CARD4L) overview. In Proceedings of the IGARSS 2018—2018 IEEE International Geoscience and Remote Sensing Symposium, Valencia, Spain, 22–27 July 2018; pp. 7407–7410.

55. John, S.; Varghese, A. Analysis of support vector machine and maximum likelihood classifiers in land cover classification using Sentinel-2 images. *Proc. Indian Natl. Sci. Acad.* **2022**, *88*, 213–227. [CrossRef]

56. Ahmed, S.; N, H. Land use and land cover classification using machine learning algorithms in google earth engine. *Earth Sci. Inform.* **2023**, *16*, 3057–3073.

57. Capolupo, A.; Monterisi, C.; Tarantino, E. Landsat images classification algorithm (LICA) to automatically extract land cover information in Google Earth Engine environment. *Remote Sens.* **2020**, *12*, 1201. [CrossRef]

58. Richardson, A.J.; Wiegand, C. Distinguishing vegetation from soil background information. *Photogramm. Eng. Remote Sens.* **1977**, *43*, 1541–1552.

59. Sezgin, M.; Sankur, B.L. Survey over image thresholding techniques and quantitative performance evaluation. *J. Electron. Imaging* **2004**, *13*, 146–168.

60. Zadeh, L.A. Fuzzy sets. *Inf. Control* **1965**, *8*, 338–353. [CrossRef]

61. Kosko, B.; Toms, M. *Fuzzy Thinking: The New Science of Fuzzy Logic*; Hyperion New York: New York, NY, USA, 1993; Volume 288.
62. Glasbey, C.A. An analysis of histogram-based thresholding algorithms. CVGIP Graph. Model. Image Process. **1993**, *55*, 532–537.
[CrossRef]

63. Doyle, W. Operations useful for similarity-invariant pattern recognition. J. ACM (JACM) **1962**, *9*, 259–267. [CrossRef]

64. Otsu, N. A threshold selection method from gray-level histograms. Automatica **1975**, *11*, 23–27. [CrossRef]

65. Kittler, J.; Illingworth, J. Minimum error thresholding. Pattern Recognit. **1986**, *19*, 41–47. [CrossRef]

66. Burney, S.A.; Tariq, H. K-means cluster analysis for image segmentation. Int. J. Comput. Appl. **2014**, *96*.

67. Panwar, P.; Gopal, G.; Kumar, R. Image Segmentation using K-means clustering and Thresholding. Image **2016**, *3*, 1787–1793.

68. Ridler, T.; Calvard, S. Picture thresholding using an iterative selection method. IEEE Trans. Syst. Man Cybern **1978**, *8*, 630–632.

69. Prewitt, J.M.; Mendelsohn, M.L. The analysis of cell images. Ann. N. Y. Acad. Sci. **1966**, *128*, 1035–1053. [CrossRef]

70. Tsai, W.H. Moment-preserving thresolding: A new approach. Comput. Vis. Graph. Image Process. **1985**, *29*, 377–393. [CrossRef]

71. Zack, G.W.; Rogers, W.E.; Latt, S.A. Automatic measurement of sister chromatid exchange frequency. J. Histochem. Cytochem. **1977**, *25*, 741–753. [CrossRef]

72. Rényi, A. On the foundations of information theory. Rev. L'Institut Int. Stat. **1965**, *33*, 1–14. [CrossRef]

73. Sahoo, P.; Wilkins, C.; Yeager, J. Threshold selection using Renyi's entropy. Pattern Recognit. **1997**, *30*, 71–84. [CrossRef]

74. Kapur, J.N.; Sahoo, P.K.; Wong, A.K. A new method for gray-level picture thresholding using the entropy of the histogram. Comput. Vis. Graph. Image Process. **1985**, *29*, 273–285. [CrossRef]

75. Pun, T. A new method for grey-level picture thresholding using the entropy of the histogram. Signal Process. **1980**, *2*, 223–237.
[CrossRef]

76. Yen, J.C.; Chang, F.J.; Chang, S. A new criterion for automatic multilevel thresholding. IEEE Trans. Image Process. **1995**, *4*, 370–378.

77. Huang, L.K.; Wang, M.J.J. Image thresholding by minimizing the measures of fuzziness. Pattern Recognit. **1995**, *28*, 41–51.
[CrossRef]

78. Shanbhag, A.G. Utilization of information measure as a means of image thresholding. CVGIP Graph. Model. Image Process. **1994**, *56*, 414–419. [CrossRef]

79. Li, C.; Tam, P.K.S. An iterative algorithm for minimum cross entropy thresholding. Pattern Recognit. Lett. **1998**, *19*, 771–776.
[CrossRef]

80. Rwanga, S.S.; Ndambuki, J.M. Accuracy assessment of land use/land cover classification using remote sensing and GIS. Int. J. Geosci. **2017**, *8*, 611. [CrossRef]

81. Yamashiro, T.; Matsuoka, S.; Estépar, R.S.J.; Bartholmai, B.J.; Diaz, A.; Ross, J.C.; Murayama, S.; Silverman, E.K.; Hatabu, H.; Washko, G.R. Kurtosis and skewness of density histograms on inspiratory and expiratory CT scans in smokers. COPD J. Chronic Obstr. Pulm. Dis. **2011**, *8*, 13–20. [CrossRef] [PubMed]

**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.