Automatic built-up area extraction by feature-level fusion of Luojia 1–01 nighttime light and Sentinel satellite imageries in Google Earth Engine

Today, with the non-stop expansion of urbanization, mapping urban areas and monitoring their dynamic changes have become challenges for governments and also a hot topic for researchers. Remote sensing imageries play a key role in urban studies, the extraction of urban built-up areas, and monitoring their changes. A variety of studies have proposed methods for the extraction of regional, national, and global built-up areas. However, the majority of them used limited features and applied a manual sample selection strategy for clas-siﬁcation, leading to time-consuming and low-eﬃcient algorithms. This paper proposes a fully automatic procedure to real-time extract built-up areas by integrating the Luojia 1–01 nighttime lights (NTL) images, Sentinel-2 multispectral data, Sentinel-1 Radar images, and SRTM elevation data in cloud-computing Google Earth Engine. Firstly, potential built-up areas (PBA) and non-built-up areas (NBA) are obtained by applying Otsu and multi-level thresholding to some of the extracted spectral-textural-spatial (STS) features and by applying logical rules. Secondly, built-up and non-built-up samples are automatically selected and are used to train a Support Vector Machine (SVM) supervised classiﬁer and to classify the hybrid feature set so that a preliminary classiﬁed map (PCM) can be obtained. Thirdly, the PCMs are automatically corrected using the non-built-up area, and morphological operations in the so-called post-classiﬁcation to provide a reﬁned classiﬁed map (RCM) and ﬁnal built-up map. Four study areas in Northern America, Europe (Scandinavia), the Middle East, and Eastern Asia were selected to test the proposed method. Also, ﬁve state-of-the-art built-up products, accompanied by Google Earth images, were used as the reference data. The results indicate that the proposed method can accurately and automatically select samples and map built-up areas with a spatial resolution of 10 m. Its performance is validated with an average overall accuracy of 94.4% and an average Kappa coeﬃcient of 0.89 and by visual comparison of our method results with other reference data. The proposed method has signiﬁcant potential to be used in real-time extracting built-up areas and in monitoring their dynamic changes on national and global scales. (cid:1) 2023 COSPAR. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Urban built-up areas, as living spaces for the majority of humans and a place for their activities, have rapidly expanded in recent decades (Jafari and Attarchi, 2021;Liu et al., 2019).Due to the fast rate of urbanization expansion, it is one of the most challenging issues for governments and policy-makers (Esch et al., 2017).With this regard, the accurate mapping of built-up areas and monitoring their dynamic changes is crucial for decision makers to manage better human lives (Liu et al., 2019).Compared to mapping built-up areas on the ground, which is so timeconsuming and expensive, remote sensing (RS) science and technology provide experts with the opportunity to accurately extract built-up areas and monitor their dynamic changes (Liu et al., 2019;Zhou et al., 2017).Therefore, accurate and fast extraction of urban built-up areas has always been a hot topic for researchers in RS (Djerriri et al., 2019).
The expansion of urban environments is taking place rapidly, and to monitor and manage these urban areas, many up-to-date and accurate maps are required.Since the methods mentioned above do not allow the user to effectively produce up-to-date maps for any parts of the Earth in a short time we developed a fully automatic framework that uses open-access RS data, real-time cloud-based image processing tools, and efficient artificial intelligence (AI) algorithms to produce urban built-up maps for any given study area and any time.For this purpose, inspired by Liu et al. (2019) and in continuation of our previous research on automatic LULC mapping (Toosi et al., 2022), in this paper, an improved, automatic, and real-time built-up area extraction method was proposed, which integrates LJ 1-01 NTL data, Sentinel-1 and Sentinel-2 imageries, and SRTM DEM.In the developed method, first, a hybrid feature space is established.Second, a variety of builtup and non-built-up training and validation samples are collected automatically using multi-level/Otsu thresholding and logical rules.Third, the training samples are used for machine learning (ML)-based supervised classification to obtain a preliminary classified map.Then, it is corrected within a two-level post-classification phase to obtain a fine classified thematic map and built-up map.By using numerous powerful features, the fully automatic proposed algorithm can reduce the confusion between built-up areas and other similar land coverages, thereby improving the accuracy of built-up area mapping.
The main contribution of this paper is to establish a fast and simple method that enables the users to automatically extract built-up areas at no charge for any desired region of interest at any time using well-established powerful ML techniques and open-access satellite data.

Study area
The following four areas from North America, Europe (Scandinavia), the Middle East, and Eastern Asia were selected as study areas: Mexico City, Mexico; Stockholm, Sweden; Tehran, Iran; and Seoul, South Korea (Fig. 1).The above-mentioned study areas which are chosen randomly are from continents with different climates, urban structures, and levels of economic condition.The choice of these four areas allows us to validate the proposed method with different built-up area distribution characteristics and to show the ability of the method for the efficient and automatic extraction of urban built-up areas.

Dataset
We used 10 m Sentinel-2 level-2A multispectral imageries, 10 m Sentinel-1 SAR images, and 130 m LJ 1-01 NTL data from the year 2018 for the selected study areas as input data.The dataset is described in Table 1 and is shown in Fig. 2. We used cloud-free Sentinel-2 images for a short period of time within a year to avoid any issues caused by any noticeable changes to built-up areas over time.Sentinel-2 images are bottom of atmosphere (BOA) or surface reflectance (SR) products which are not affected by the atmosphere.The preprocessed Sentinel-1, Sentinel-2, and SRTM DEM data were called back from the Google Earth Engine (GEE) dataset repository (https://developers.google.com/earth-engine/datasets/),while the LJ 1-01 images were downloaded from the official website of Hubei data and application network (https://www.hbeos.org.cn;https://59.175.109.173:8888/).

Processing tool
The GEE cloud computing platform is used to access a variety of open-source geospatial data without the need to download them.On the other hand, some levels of preprocessing (i.e.geometric and atmospheric corrections) were applied by the developers to the data which can facilitate working with big satellite data.

Methodology
According to the block diagram in Fig. 3, the proposed method of this study is made up of three steps including ''data preparation and feature extraction," ''classification and post-classification," and ''validation" which are separately described as follows:

Data preparation and feature extraction
In this initial step, the Sentinel-1 SAR images are called back from the GEE data repository, filtered by the boundary (for each study area) and date (for the year 2018), and all the available polarization/bands (in these case studies VV and VH are available) are selected.Then, the subtraction of VV and VH (VV-VH) and their average ( VVþVH 2 ) are calculated as features since they can help to discriminate built-up areas from other land coverages (Lino et al.,  2018).For Sentinel-2 data, the images are called back in GEE and filtered by the boundary of study areas.Then, the remaining images are filtered by date (the summer season was set for the date duration to have the least cloudy or foggy pixels) and ''percentage of cloudy pixels" where only the images with cloudy pixels less than 1% are selected.All the spectral bands having a spatial resolution of more than 10 m are manipulated (resampled) so that all the channels can reach a spatial resolution of 10 m (since the Sentinel-2 dataset lacks a high-resolution panchromatic band that has complete spectral coverage with other bands, no pansharpening operation has been carried out).Since built-up areas have different textures compared to other land coverages, the first group of features that are obtained from the average of three original 10 m-visible (RGB) bands of the Sentinel-2 dataset is textural indices which are used for discriminating built-up from non-built-up areas.We calculated 11 commonly used textural descriptors including mean, variance, homogeneity, contrast, dissimilarity, second moment, correlation, prominence, inertia, shade, and entropy which are calculated from Gray-Level Cooccurrence Matrices (GLCM) (Haralick et al., 1973).The formulas related to the textural indices are provided in Table 2.
The SRTM slope maps and the features derived from the arithmetic combination of Sentinel-1 Radar bands form the spatial feature set.Due to their relatively coarse spatial resolution, the SRTM slope maps are not included in the hybrid feature set (which will be imported into the classification), and they are only used for masking mountainous areas.It is noteworthy that as the radar backscatterings are sensitive to the spatial pattern of the phenomenon, we categorized them as spatial features, and as their main duty, they can help to identify built-up areas (Semenzato et al., 2020).Furthermore, as their sec-ondary task, radars can effectively help us in discriminating water and vegetation from other objects (Bioresita et al., 2018).
The Otsu thresholding is used to provide a threshold to obtain water, no-water, vegetation, and no-vegetation masks.It segments the indices by exhaustively searching for the optimum threshold (T) that minimizes the intraclass variance (Eq.( 22)).
ð Þ are variances of the two desired classes and w 1 ¼ P TÀ1 i¼0 pðiÞ and w 2 ¼ P LÀ1 i¼0 pðiÞ are the probabilities of the two classes separated by T (Sezgin et al., 2004;Toosi et al., 2022).
As the LJ 1-01 NTL images are not covered by the GEE data repository, it is necessary to freely download them from the official website of the Hubei data and application network and then upload them to the GEE platform.Regarding the breadth of products, image mosaicking is necessary for case studies extending more than 250 km Â 250 km.For some case studies, there may be slight spatial misregistration between the LJ 1-01 and Sentinel data that in such cases LJ 1-01 should be co-registered to Sentinel-2 10 m bands, for instance.The basic format of data is the digital numbers (DN) which are converted to radiance (q) according to the exponential relation of Eq. ( 23).
where q is the input radiance (W=ðm 2 srlmÞ) and DN is an abbreviation for the digital number of the LJ 1-01 images.Afterward, multilevel thresholding is applied to q in order to obtain q = 0 (completely dark areas on the Earth) and q !T (T is a threshold that presents highly lit-up areas on the Earth1 ).The reason is that the NTL images are correlated with human activities and that built-up areas have higher DN values, compared to natural non-built-up areas with very low DN values (and sometimes the values are equal to zero) (Liu et al., 2019).It is noteworthy that the LJ data are not significantly affected by different types of noises (Li et al., 2019).
SRTM Digital Elevation Data (DEM) Version-4 is called back in GEE, and its elevation band is filtered by boundary and date.Then, the slope map is calculated from the DEM layer.T = 15 is used as a threshold to obtain mountain and no-mountain masks.T = 15 is the most common threshold, which is used in studies conducted for built-up mapping (Ban et al., 2015).The AND and OR logical rules are applied to all the masks obtained from the SRTM, LJ 1-01 and Sentinel-2 data to make two base regions called potential built-up area (PBA) and non-builtup area (NBA).The AND operator states that PBA is an area that has high light-up and at the same time is neither covered by vegetation nor located on steep slopes.In addition, the OR operator indicates that NBA is the summation of completely dark (with no nighttime light) areas or areas with steep slopes or covered with water/vegetation (at least one of the conditions is needed so that it can true).On the other hand, all the spectral, textural, and spatial features are stacked together to obtain a hybrid feature space.

Classification and Post-classification
A variety of sample points are randomly and automatically selected in the PBA and NBA regions.The samples are generated based on the GEE ''RandomPoints" function which generates points that are uniformly random on the sphere and are within the given region (GEE Ref., 2022).For random sample generation, the given polygon, which are our study area, is decomposed into n different triangles (i ¼ f1; 2; 3; Á Á Á ; ng) with vertices V 1i , V 2i , and V 3i .Therefore, a random point (P j ) can be generated uniformly within triangle V 1i V 2i V 3i according to the convex combination of the vertices (Eq.( 24)).
where w 1i and w 2i , which are called convex weights, are uniformly and randomly drawn from the range of [0,1] for the ith triangle.It is also noteworthy that always , 2011;Osada et al., 2002;SOF, 2011).
A preliminary sample collection is obtained, which can also be improved by human supervision/revision.It is noteworthy that a manual sample selection module is also set up, which allows the human expert to manually add some samples or delete some of them.After the final sample collection is prepared, it is then split into training and validation samples in an 80%-20% ratio.The training samples are used to train the Support Vector Machine (SVM) as one of the most fashionable ML supervised classification methods (Yang et al., 2018), and the validation samples are used to evaluate the classification performance.It is worth noting that here, traditional ML methods are used instead of DL-based methods so that they can be suitable for real-time, fast (as one of the goals of our proposed method), and precise processing.This is because despite the relatively higher accuracy of DL methods (Karaca et al., 2017), training of DL models is so time-consuming in comparison to traditional methods and the time will increase drastically with the increase in the size of study areas (Yilmaz et al., 2020).It is also noteworthy that in many cases, the well-established traditional ML methods, such as SVM, can provide comparable results to the DL methods (Lai, 2019).The classification process provides a preliminary classified map (PCM).Post-classification operations are set up to improve PCM.Post-classification is conducted in two steps including (1) overlaying NBA to the PCM, and then (2) applying morphological operations, i.e. opening and closing.In the first step, the NBA is overlaid to the PCM (Eq.( 25)), i.e. an AND logical operation is performed between the PCM and NBA binary images.Assume the pixels in the PCM that are related to builtup and non-built-up classes be shown by 1 and 0, respectively.Also, the pixels with 1 in the NBA indicate the pixels that are definitely non-built-up.Thus, the refined classified map (RCM) obtained by overlaying (RCM OL ) is a binary image in which the misclassified pixels are corrected.
where | denotes conditional equality, which means that RCM is the same as PCM with the condition that the statement in the braces is exerted on the PCM.In the second level, RCM OL is refined by morphological opening and closing (Eq.( 26)) to remove small misclassified pixels (speckles) and to close (fill) the gaps and holes in the classified map.
where RCM F is the final version of RCM.The opening, closing, erosion, and dilation are presented by ''o, _ s, É, and È", respectively.Also, H represents the structural element.After post-classification, the final built-up map is obtained, and the area of built-up area can be calculated using the obtained maps, if necessary.

Validation
Validation is done in both quantitative and qualitative manners.In the former method, the evaluation is conducted using the elements of confusion (error) matrix, while in the latter, the assessment is done visually by comparing the obtained results with base high-resolution Google Earth images, state-of-the-art reference maps, and the result of other studies.In our error matrix, rows and columns correspond to reference and predicted classes (BU: built-up and NBU: non-built-up), respectively.The elements of the error matrix are the User's Accuracy (UA), Producer's Accuracy (PA), Overall Accuracy (OA), and Kappa Coefficient (KC).UA indicates the probability that a classified built-up pixel is truly a built-up area.On the other hand, PA calculates the probability that a built-up pixel is correctly identified.OA indicates the proportion of pixels that are correctly identified.Finally, KC is an indicator that indicates the agreement between classified results and ground truth.UA, PA, OA, and KC are defined as Eqs.( 27)-(30), respectively.
where X ii represents the number of pixels that are correctly identified, X iþ is the number of pixels that belong to class i, X þi is the number of pixels identified as class i, and N is the total number of image pixels (Li et al., 2018;Pearson and Blakeman, 1904).

Results and discussions
Fig. 4 presents the overall hybrid spectral, textural, and spatial (STS) feature space including 65 feature images in four study areas.In the case of textural descriptors, the descriptors are obtained by setting the parameters as follows (based on the GEE's default): The size of the neighborhood is set to 1, and a 3 Â 3 square window is utilized, leading to four GLCMs with the offsets (À1,À1), (0,À1), (1,À1), and (À1,0).It is noteworthy that to avoid the low-quality STS features, i.e. features containing high noise or with relatively low contrast, the quality of features can be assessed both visually and numerically using the simple no-reference quality assessment metrics such as Standard Deviation, BRISQUE, DIIVINE, BLIINDS-II, SSEQ, BIQI, etc. (here we used the well-established Standard Deviation metric) (Al-Wassai and Kalyankar, 2012; Javan et al., 2019;Mhangara et al., 2020) This led to keeping only the high-quality STS features (that properly discriminate between the phenomena) and deleting useless features from the feature set.
As mentioned before, built-up indices form a significant portion of spectral index groups.According to Fig. 5, which shows the correlation between the corresponding pixels of built-up indices for SA1 (as an example), it can be seen that except in a few cases (with a weak correlation or R%0) there is a significant positive correlation (R) between different pairs of indices at the level of a = 0.01 in a two-tailed statistical test.What is more important than the correlation of indices is the ability of indices to discriminate between different land covers (the high power of the indices to perform such a task is indicated in Fig. 4).
Fig. 6 shows the obtained PBA and NBA accompanied by the Otsu thresholding outputs, i.e., water mask and vegetation mask.The thresholds of the histograms, i.e.T1-T8 are the Otsu thresholds that segment the water-covered and vegetated areas.The sample points including training feature collection and validation feature collection for both built-up and non-built-up classes are presented in Fig. 7. Table 4 contains information on the number of sample points.It is noteworthy that the number of samples is set by trial and error, which provides the best outcomes.Due to the computational load issue and the more accurate selection of samples and avoiding the selection of mixed pixels instead of pure pixels, a point-wise selection method was used instead of the polygon-based strategy.Hence this strategy can lead to a relatively faster classification process and more accurate and reliable results (Toosi et al., 2022).Furthermore, the manual sample selection method for large-scale case studies such as provinces, countries, and continents is so time-consuming; in some cases, humanbased errors may emerge.Thus, we believe the automatic sample selection strategy of our method seems to be more efficient, which decreases the number of mistakes.
In order to tune the SVM supervised classifier parameters, the default values of GEE's references are considered, i.e.Radial Basis Function (RBF) is chosen as the kernel type, and c and cost parameters are set to 0.5 and 10, respectively.SVM is trained by training features, and then the hybrid feature set is introduced to the classifier.By producing the initial classified maps, the produced maps are morphologically corrected using opening and closing operations with circle-shaped normalized Boolean structural elements (kernels) of radius 1 to obtain the RCMs.The mentioned 3 Â 3 kernel is expressed as [0, 0.2, 0; 0.2, 0.2, 0.2; 0, 0.2, 0].In this matrix, the semicolons separate the rows, and the colons separate the elements in each row.The RCM F and the final built-up map for all study areas are produced, which are shown in Fig. 8.
The confusion (error) matrix for the classification performance in all study areas is presented in Fig. 9.The average overall accuracy (OA) and Kappa coefficient (K) were calculated as 94.4% and 0.89, respectively.According to the USGS standard limit (OA !85 %), all the classification   outputs are acceptable (Chen et al., 2018c).The obtained results of our proposed method were also visually compared with some other state-of-the-art national and global built-up maps (Fig. 10).The products include 1000 m-Global Human Settlement Layers (GHSL) (Pesaresi et al., 2016), 100 m-Copernicus Global Land Cover Layers, i.e., CGLS-LC100 Collection 3 (CGLS) (Buchhorn et al., 2020), ESA WorldCover 10 m v100 (ESA WC) (Van De The results showed that all the methods similarly and effectively extracted the general pattern of urban built-up areas and that the slight differences are related to some challenging areas such as bare lands (which are similar to built-up areas from the visual/spectral similarity point of view) and narrow road network.By scrutinizing the results and comparing them with each other, it is found that our proposed method effectively extracts built-up areas and  that even it has a good ability to map narrow roads.Furthermore, the GHSL and ESA WC, with a spatial resolution of 30 m and 10 m, respectively, are two of the other powerful methods that can effectively extract urban road networks.Among the four high-resolution (10 m) builtup maps, i.e. our produced maps, CLC, and DW V1, and WSF, the first two better match the Google Earth ground truth, while in DW V1 and WSF there are some slight discontinuities which are less consistent with reality.Due to its spatial resolution (100 m) and with an OA of about 80% (which was declared in Buchhorn et al., 2020), the CGLS product can extract the overall pattern of urban built-up areas, and it can model the road structure with a width less than 100 m.Furthermore, with a resolution of about 500 m, MCD12Q1 is not able to extract small built-up areas, i.e. like CGLS only mapped the urban area in a relatively coarse scale.Our results and the WSF product have the most similarity in terms of visual characteristics, i.e. in both of them, there are the least possible discontinuities in the final map.The two-level post-classification in our method provides us with such results with the least unreal discontinuities.For the SA3, the result of our method is compared to 10 m-ILCM.According to Ghorbanian et al. (2020), ILCM was produced with an OA of 95.6%, which shows that our method had good accuracy in the extraction of built-up areas in SA3 (Note that our method has an OA of about 95.5% in SA3).It can be seen that the OA and Kappa in SA4 0 s results are smaller than other SAs.
To investigate the reason for this, it can be argued that due to the nature of our proposed method, the algorithm has been run several times and different amounts of training and validation samples have been set each time.In order to analyze the impact of the classifier algorithm in the classification process, the statistical difference of the SVM result is compared to that of the other well-known classifier, i.e.Random Forest (RF).For such a comparison, the non-parametric McNemar's statistical test is used, which defines a standardized normal indicator as Eq. ( 31) (McNemar, 1947;Samadzadegan et al., 2017).
where f 12 indicates the number of samples classified correctly by SVM and incorrectly by RF, and f 21 indicates the number of samples classified correctly by RF and incorrectly by SVM.
The test is carried out on the Seoul case study (SA4) in which more samples (compared to the other three study areas) are involved in its classification process.For the same training and validation samples (similar to that was used for SVM), the confusion matrix was obtained as Eq.(32) and Eq. ( 33) for SVM and RF, respectively.Thus Z = 1.63 indicates that the difference in accuracy between the two classifiers is not statistically significant (as |Z|less than1.96).Furthermore, as Z greater than 0, it is inferred that SVM is relatively more accurate than RF (Samadzadegan et al., 2017) One of the most prominent strength points of our realtime method, compared to its competitors, is its high flexibility.By the term ''flexibility," we mean that one can easily manipulate the hybrid feature set (increase or decrease the layers), the size of sample collections, and the key parameters as fast as possible and re-run the algorithm for the same (or even other) case study without running into any trouble.Furthermore, the joint automatic and manual modules for sample selection led to producing a sample feature collection from all kinds of built-up characteristics.This provides a diverse and rich set of training data for training the SVM classifier and prevents overfitting in the classification of new inputs (as the classifier is trained with a variety of diverse samples from all over the study area and it is somehow saturated).In order to guarantee the robustness of our method against the study area, i.e. to investigate if its performance depends on the study area and also to find out whether or not it works for any area with different built-up characteristics, it should be implemented for a large number of study areas (e.g. for tens of cases) all over the world.
Our proposed method was implemented on a laptop with a 2.50 GHz Pentium Intel Core i7 CPU, NVIDIA GeForce MX 130 GPU, and 16 GB of RAM.However, Internet speed is very more important than hardware characteristics.With a 10 Mb/s Internet speed, GEE codes consume about three minutes to be run for the study areas of this paper.As a result, the performance on national and global scales would be time-consuming while having about 1,000-3,000 samples and a hybrid feature set containing more than 15 STS bands.In the case of online cloudcomputing server-client platforms, such as GEE, where a great portion of computational affairs are conducted on the server side and only the final desired results are sent to the client side, the process gets accelerated considerably, compared to the traditional offline desktop-based image processing platforms.
As the LJ 1-01 satellite was launched very recently (in June 2018), its images do not cover the whole globe (Fig. 11 shows that by the year 2022, about 30% of the world is covered by 8,673 NTL image frames), and also they have not been provided by the GEE platform yet.In the upcoming years, by having the LJ 1-01 images in the GEE data repository, it can be a good opportunity to automatically extract the built-up areas for the whole globe with a spatial resolution of 10 m and also to monitor dynamic changes of built-up areas within two epochs in real-time mode.

Conclusions
In this study, we proposed a fully automatic method to integrate the LJ 1-01 NTL data, Sentinel-1 SAR, and Sentinel-2 multispectral imagery for real-time mapping of urban built-up areas in the GEE platform.A hybrid STS feature space was established and was accompanied by the training and validation sample collections, and then they were entered into the well-established SVM classification.The samples were automatically collected based on the LJ image and the water and vegetation masks generated by the Sentinel-2 imagery and applying thresholding.The hybrid STS feature set was classified to produce a preliminary built-up map.Then, the mentioned map was automatically improved in the two-level post-classification process.The proposed method with an average OA of 94.4% was shown to efficiently extract the built-up area.Also, compared to other state-of-the-art methods, our algorithm has promising performance.
One of the good characteristics of our method is using the cloud computing GEE platform, and this makes it suitable to be implemented in other case studies all over the globe.In this study, GEE was shown to have a brilliant capability to process multi-gigabyte big satellite data and produce desired built-up maps in a short time.Our method may face some challenges.For instance, although we filtered the input images to have the least possible cloudfree data, the cloud may cause some problem for the results by affecting the Sentinel-2 image, and the challenge may intensify when the study area drastically increases (e.g., for the whole globe).Furthermore, obviously due to the limitations in the input data, our method is not able to map the urban roads with a width lower than the spatial resolution of the Sentinel-2 images (10 m).
For future research to improve the performance of our method, one can implement the algorithm using the DLbased method (instead of traditional ML-based methods such as SVM, RF, etc.) by linking the GEE and Google Colaboratory (Colab) platforms.In order to make other improvements in future research, one can use the superresolution (SR)-based DL techniques to sharpen the initial Sentinel MS bands using their own spatial information (Zhong et al., 2016) or the high-resolution SkySat and NAIP satellite and aerial images provided by GEE2 (EEDC, 2022b).The image sharpening led to producing enhanced Sentinel data with a spatial resolution better than 10 m.It is worth noting that in Zhong et al.'s (2016) selfenhancement technique the resolution of the Sentinel images improves using the spatial information of the Sentinel multispectral bands themselves.Furthermore, a distinct module can be set up for extraction of road networks out of the LJ 1-01 NTL data by distinguishing the urban regions through a threshold-based method and by means of an unsupervised pulse coupled neural network (PCNN) (Wang et al., 2021).The output of the road extraction module can be aggregated with our base built- up mapping method using an AND logical operation to obtain the final built-up areas.
We hope that this research will open new horizons to the RS scientific community and that its promising results will be useful for them.Furthermore, it is desired that, according to these geospatial maps, the decision makers make appropriate and effective decisions toward sustainable development.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

F
Fig. 6.The results of Otsu thresholding and the obtained PBAs and NBAs.

Fig. 8 .
Fig. 8.The RCM F and final built-up maps produced by our method.

Fig. 11 .
Fig. 11.The global coverage of the LJ 1-01 NTL data (source: Hubei data and application network).

Table 1
Data descriptions.
* P(i,j): the relative frequency with which two pixels occur within a given neighborhood; i and j: pixels' intensity values; G: number of grey levels.

Table 4
Number of training and validation samples*.
. The indicator revealed that our results are not much dependent on the classifier method and that any well-established classifier with such a large number of training samples and rich datasets can provide acceptable results.