Landsat Images Classification Algorithm (LICA) to Automatically Extract Land Cover Information in Google Earth Engine Environment

: Remote sensing has been recognized as the main technique to extract land cover/land use (LC/LU) data, required to address many environmental issues. Therefore, over the years, many approaches have been introduced and explored to optimize the resultant classification maps. Particularly, index-based methods have highlighted its efficiency and effectiveness in detecting LC/LU in a multitemporal and multisensors analysis perspective. Nevertheless, the developed indices are suitable to extract a specific class but not to completely classify the whole area. In this study, a new Landsat Images Classification Algorithm (LICA) is proposed to automatically detect land cover (LC) information using satellite open data provided by different Landsat missions in order to perform a multitemporal and multisensors analysis. All the steps of the proposed method were implemented within Google Earth Engine (GEE) to automatize the procedure, manage geospatial big data, and quickly extract land cover information. The algorithm was tested on the experimental site of Siponto, a historic municipality located in Apulia Region (Southern Italy) using 12 radiometrically and atmospherically corrected satellite images collected from Landsat archive (four images, one for each season, were selected from Landsat 5, 7, and 8, respectively). Those images were initially used to assess the performance of 82 traditional spectral indices. Since their classification accuracy and the number of identified LC categories were not satisfying, an analysis of the different spectral signatures existing in the study area was also performed, generating a new algorithm based on the sequential application of two new indices (SwirTirRed (STRed) index and SwiRed index). The former was based on the integration of shortwave infrared (SWIR), thermal infrared (TIR), and red bands, whereas the latter featured a combination of SWIR and red bands. The performance of LICA was preferable to those of conventional indices both in terms of accuracy and extracted classes number (water, dense and sparse vegetation, mining areas, built-up areas versus water, and dense and sparse vegetation). GEE platform allowed us to go beyond desktop system limitations, reducing acquisition and processing times for geospatial big data.


Introduction
Accurate maps of land cover/land use (LC/LU) distribution are essential to gather information which is useful in many land management and environmental monitoring tasks. Therefore, over the last 15 years [1], several products have been generated to face the growing demands for related maps, using different approaches. Among these, the remote sensing technique has been an invaluable source of LC/LU information [2,3]. However, most of the satellite-derived maps covering the whole world have a coarse resolution, not suitable to describe the true Earth heterogeneity and urban and agricultural landscapes. For instance, the Global Land Cover product (GLC2000) carried out by the Basically, the classification algorithms are grouped into two categories: Unsupervised and supervised approaches. The former aggregates the pixels of an image in classes by analyzing the similarity of attributes, without any analyst's contribution [16]. Such methods are commonly applied when the knowledge about the land cover is scarce. By contrast, the operators' work is a key factor for the second group: they identify some training areas to coach the algorithms and to assign each pixel of the images in a specific category [22]. Although the first approaches are completely automatized, they are extremely time consuming since they require operator input to improve the accuracy of a classification map. However, supervised classification is not error-free either and the analyst has to refine the outcomes.
Among the several methodologies developed to extract LC/LU information belonging to both groups, index-based approach method [27], maximum likelihood supervised classification (ML) [28], machine learning algorithms (MLAs) [29] and object-based image analysis (OBIA) approach [30][31][32][33] are the most popular. Yet, each of them shows some strengths and weakness [34]. Although the indexbased approach allows us to reduce the amounts of components and to classify a large area in a short time, several indices must be applied to detect the different LC/LU classes since each of them is aimed at distinguishing just one category [35]; for instance, vegetation indices are intended to identify "green areas" and so on. ML is recognized as one of the simplest algorithms to implement and to interpret [36], but its results are not satisfying without introducing a large amount of training areas since, because of insufficient a priori information, it assumes an equal a priori probability for each land cover classes [29]. Completely opposite are the MLAs, which comprise different approaches, such as artificial neural networks [37], support vector analysis (SVA) [38], and random forests (RF) [39]. Nevertheless, although these algorithms are efficient [40] and show more accurate results than the other conventional methods [29,41], MLAs are difficult to be implemented since, generally, a large volume of parameters must be fixed [40]. Moreover, MLAs tend to over-fit data [40]. There are some exceptions since each MLA shows peculiar traits and reveals different performances. Conversely to the other approaches included in the MLA group, SVA requires a smaller number of data training [42] and RF does not over-fit data due to the law of large numbers and allows us to reduce training dataset size with the consequent increment of overall error [43]. In contrast to the other methods, OBIA classification is based on the integration of spectral and geomorphological factors, which increase the accuracy of the resultant classification map [44]. Nevertheless, its outcomes look really promising if medium-or fine-resolution data are used as input [45].
Thus, none of the listed techniques allows us to generate optimal outcomes in all conditions. Therefore, the approach to apply should be selected considering multiple aspects, such as data type, spatial resolution, accuracy, operator skills, speed, classifier interpretability, and knowledge of ground truths. In [13], it showed that an index-based classification approach is efficient and effective for automatically extracting LC/LU information in multitemporal and multisensory analysis perspectives. The index-based approach involves the combination of two or more spectral bands, in order to classify Earth's features. Each coverage, indeed, showed a specific spectral signature, commonly recognized as their fingerprint, according to their ability of absorbing, transmitting, and reflecting the energy [28]. Thus, properly integrating particular wavelengths, distinctive of a specific element, allows us to detect LC/LU classes. Although several indices have been introduced in literature, we are still lacking an index-based method suitable for classifying the whole study area by using different Landsat satellite images. In fact, each index is based on the integration of different spectral bands in order to address a specific need and to extract a certain LC/LU class [46][47][48][49].
The objective of this paper is to introduce a new classification algorithm to process Landsat images (Landsat Images Classification Algorithm: LICA) in GEE environment to automatically extract LC/LU information. This method was implemented after the analysis of the performance of 82 indices, commonly applied in literature, to detect land cover classes processed in a more efficient way and by increasing the accuracy of final results. LICA is composed by the computation of two new indices, SwirTirRed (STRed index) and SwiRed, introduced in this paper for the first time: The former aimed to detect water, mining areas, and sparse and dense vegetation while, the latter, builtup areas. LICA reliability was tested on the pilot site of Siponto using 12 Landsat images, belonging to missions 5, 7, and 8, as input data. Those images were acquired in different seasons and years, covering a period of about 17 years, in order to demonstrate that it produces a baseline information suitable for performing multitemporal, multiseasons, and multisensory change detection analysis.

Study Area
The method was tested along the coastline of Siponto in the Apulian Region (Southern Italy), studying an area bordered by the Mediterranean relief of Gargano to the north, the marshland to the south, the Candelaro estuary river and the Adriatic Sea to the west and east, respectively ( Figure 1). The area, located about 2 km far from the city center of Manfredonia, was selected as an experimental site both because of its historical relevance and the changes suffered by its landscape over the years. This choice allowed us to test the performance of the proposed algorithm and to assess its accuracy on a zone characterized, over the years, by different features, configurations, and issues, such as the erosion process. Founded in 194 BC, Siponto became a crucial commercial and maritime hub during the Roman period, as proven by the Archaeological Park of Siponto. Its relevance gradually slackened as a result of the depopulation process that followed the swamping of its seaport and two devastating earthquakes in 1223 and 1255. From then on, as highlighted by [49], its territory was essentially earmarked to agricultural purposes, exploiting the dense network of irrigation ditches available in that environment. This trend was only inverted over the last few years as tourism started to develop, being encouraged by the beauty of the local landscape and favorable climate conditions. These elements were not the only triggering factors of the soil erosion process suffered by this area. The construction of the new port in Margherita di Savoia in 1952 was, in fact, currently recognized as its main cause [50]. Although such problems are well known and about 80% of the shoreline conservation activities performed in the Apulia Region have addressed the investigated area, erosion issues are still not solved [50].

Landsat Image Classification Algorithm (LICA)
Classification methods allowed us to generate thematic maps, assigning each pixel to the proper belonging class. As proposed by [46], an index-based approach was efficient to quickly reveal LC/LU classes from satellite images and, therefore, in this case, it was preferred to other classification approaches. By mixing spectral bands' information, spectral indices are able to bring out Earth's features capacity in absorbing, reflecting, and transmitting the energy [51]. For this purpose, 82 consolidated indices, commonly applied in literature, were computed to extract LC/LU information (Table 1). Twenty-six of them were selected to detect bare soil and built-up areas, while the remaining 56, called vegetation indices (VIs), were aimed at identifying vegetation. Conventional indices were tested to bring out the potentiality of the strongest and weakest bands in extracting land cover types by verifying their reliability in the area under investigation. While all the algorithms were easy to implement, just three of them provided accurate results, i.e., Optimized Soil Adjusted Vegetation Index (OSAVI) [52] (Equation (1)), Green Optimized Soil Adjusted Vegetation Index (GOSAVI) [53] (Equation (2)), and Normalized Difference Bareness Index (version 2) (NDBaI2) [54] (Equation (3)).
where NIR is the near-infrared band, R is the red component, G is the green band, SWIR and TIR are the shortwave infrared and the thermal infrared bands. The first two indices (OSAVI and GOSAVI) are included in the VIs group and, consequently, they are suitable for classifying dense and sparse vegetation. Conversely, NDBaI2 can correctly classify a higher number of categories: Built-up areas, mining areas, water, bare soil, and dense and sparse vegetation. Considering the number of LC/LU classes detected by each index and their best overall accuracy (Table 1), NDBaI2 appeared as the most reliable index and was consequently used as the starting point to develop LICA procedures. NDBaI2 is based on the combination of SWIR1 and TIR1 (Equation (3)) and, therefore, this led us to believe that those bands should be the most essential to classify the whole study areas. This consideration was also supported by literature review since LC/LU classes are strongly affected by TIR [48] and SWIR, usually applied to distinguish bare soil and built-up areas [55,56]. Moreover, SWIR also allowed us to distinguish sparse and dense vegetation because of its dependency from the amount of water content in leaves [50,51]. Then, [57][58][59] enhanced the importance of the red band since it is linked to the energy absorbed by chlorophyll. In addition, these data were also integrated with the information retrieved through the spectral signatures' examination of each LC/LU category existing in the study area ( Figure 2). SWIR1 band showed a great difference among mining areas, water, and sparse and dense vegetation. On the contrary, TIR1 displayed different values among water, bare soil, mining, and built-up areas ( Figure  2). In addition, Figure 2 enhances the contribution of red band as well to distinguish water, bare soils, mining, and built-up areas.  Therefore, SWIR, TIR, and R were integrated to classify water, mining areas, and sparse and dense vegetation. Conversely, just SWIR and R were combined to detect built-up areas. The first index, called SwirTirRed index (STRed index), is reported in Equation (4). The second one, named SwiRed index, is described by Equation (5).
The workflow of Landsat Images Classification Algorithm (LICA) is reported in Figure 3. LICA is generated by the sequential computation of the two introduced new indices (STRed and SwiRed index) on the outcome of a cloud masking procedure, performed on atmospherically corrected Landsat images. Once their implementation was completed, thresholds to identify each LC/LU class were set (Table 2) and the resultant maps were merged. Figure 3 describes the suggested workflow to be set.

Database Construction in GEE Platform
GEE (https://earthengine.google.com/) is a cloud computing environment designed and released by Google in the last few years to overcome desktop platforms' limitations related to the storage and the management of a huge amount of geospatial data [25]. Such a platform is characterized by a dedicated high-performance computing (HPC) infrastructure that provides an interactive developing environment directly connected to the available open data, such as Landsat and Sentinel images archive, as well as digital elevation models, vector, socio-economic, topographic, and climate layers sets [20]. Therefore, these data can be directly downloaded both in raw and preprocessed format, minimizing their acquiring and processing time, in GEE platform. To meet the purpose of our research, 12 scenes covering a period of 17 years, from 2002 to 2019, radiometrically and atmospherically corrected, belonging to LANDSAT missions 5, 7, and 8, referring to the experimental area of Siponto, were selected (Table 3). Particularly, four images were collected for each mission, each of them belonging to a different season (winter, spring, summer, and fall). The collected images were provided in the Universal Transverse Mercator (UTM) projection and the World Geodetic System (WGS84) datum. As shown in Table 3, cloud cover information was also considered: Only scenes characterized by a cloud cover value lower than 20% were taken into account in the data selection phase. Where needed, clouds were subsequently masked through the adoption of proper filters, based on the exploitation of the information provided by the quality assessment (QA) band, already implemented in GEE, as suggested by [122] and [123]. In this way, the cloudy pixels were rendered transparent and, therefore, excluded from further algorithm implementation.
On the contrary, selected images were not orthorectified since the geometric accuracy provided by USGS was satisfactory. Therefore, the developed classification algorithm was directly computed on the outcome of the cloud cover masking procedure, as described in the workflow reported in Figure 3. Landsat archive analysis, cloud masking process, and all the further processing phases were performed on the cloud, exploiting GEE interactive environment.

Implementation of Classification Indices and LICA in GEE
Once the images were downloaded and preprocessed, the JavaScript application programming interface (API), implemented in the GEE, was used to integrate the spectral bands and estimate the indices, commonly used in literature to classify satellite images. In Section 2.2 the calculated indices were described in detail. The documentation for combining spectral bands is reported at https://developers.google.com/earth-engine (accessed 2 September 2019). Subsequently, the proposed workflow (Figure 3) for automatically classifying Landsat images was implemented and LICA images were then generated. Class distinctions were obtained using LICA thresholds ( Table 2).

Strategies to Evaluate the Accuracy
A multitemporal reference dataset based on a stratified random sampling point was generated to assess the accuracy of the proposed approach [124,125]. A total of 11,245 pixels as testing samples, proportionally distributed in each class according to their extension, were selected. Therefore, 1328 pixels were used to verify the accuracy of water, 492 pixels for built-up areas, 151 pixels for mining areas, 3165 pixels for mining areas, and 755 and 924 pixels were implemented to verify the accuracy of sparse and dense vegetation categories, respectively. Subsequently, a manual interpretation was performed to label samples according to their allocation. Samples were overlapped on the corresponding original Landsat data, manually interpreted in order to detect land cover information, and assigned to a specific class. This procedure was separately implemented on each resultant classification map.
The metrics of overall accuracy (OA), producer's accuracy (PA), and user's accuracy (UA) were next computed to perform a per-pixel accuracy assessment of classification procedure outcomes [29,[126][127][128][129]. OA, PA, and UA showed a value between 0 and 1: The higher the values, the better the accuracy.
Finally, the performance of the introduced algorithm was compared to the one achieved by each index commonly applied in literature to verify its advantages and disadvantages.

Classification Results
This section is dedicated to the classification procedure outcomes obtained through the application of indices consolidated in literature (Figures 4-6) and the proposed LIC algorithm (Figures 7 and 8). Traditional indices didn't show satisfying results, except for OSAVI, GOSAVI, and NDBaI2, which were presented. Moreover, since their performance was similar for all the Landsat missions considered, for the sake of brevity, just the outcomes generated from the processing of Landsat 8 (17 March 2019) are reported.
OSAVI algorithm distinguishes three classes (water, and dense and sparse vegetation) ( Figure  4). Nevertheless, the classification was not accurate since some misclassified pixels could be pinpointed between dense and sparse vegetation, as highlighted on the right side of Figure 4. This means that it cannot correctly detect different types of vegetation, its density, or health status. This is confirmed by analyzing the accuracy of its performance, reported in the following section (see Section 3.2). GOSAVI algorithm demonstrated a similar trend, as it could only distinguish three classes (water, and dense and sparse vegetation) as well. Like OSAVI, it presented some misclassified pixels, reported on the right side of Figure 5, yet it did not show problems in classifying dense and sparse vegetation. This improvement was due to the introduction of a green band in the OSAVI computation to register the information of leaf pigments. The observed issues were related to water detection. Its classification accuracy is reported in Section 3.2. In contrast to OSAVI and GOSAVI, NDBai2 allowed us to detect more classes: In addition to water, and dense and sparse vegetation, mining areas and built-up areas were also distinguished ( Figure 6). Despite the improved performance, its accuracy was lower and some issues were detected on the resultant map: Built-up areas were generally classified as mining areas, whereas dense vegetation was confused with sparse vegetation and water ( Figure 6). This was confirmed by its confusion matrix (see Section 3.2) As described in Section 2.2, LICA consisted of two different steps: The former intended to classify water, mining areas, and dense and sparse vegetation (Figures 7-9); the latter aimed at identifying built-up areas (Figures 10-12). The first phase was performed by applying the new STRed index, while in the second phase the novel SwiRed index was implemented. Thus, the resultant maps of the proposed algorithm provided information on the same number of classes retrieved by NDBaI2 ( Figure 13). However, LICA showed higher accuracy than NDBaI2, as demonstrated through the confusion matrix described in the following section, since misclassified pixels were drastically reduced.        Tables 4-12 provide OA, UA, and PA of resultant classification maps obtained through the computation of OSAVI, GOSAVI, and NDBIaI2 on the 12 atmospherically corrected Landsat data. On the contrary, just the best OA related to the outcomes generated by the remaining 78 indices are shown in Table 1. Although the best OA value was quiet high for the three indices (88.91, 89.89, and 82.59, respectively), their accuracy matrices bring out the difficulties encountered in classifying the study area, e.g., OSAVI incorrectly identified sparse vegetation pixels; similarly, NDBIaI2 can just detect 30% of pixels included in built-up areas. Therefore, although their results were satisfying, they cannot be used to extract accurate information related to the land cover of the experimental area.       Table 11. OA, PA, and UA obtained through the application of Normalized Difference Bareness Index (version 2) (NDBIaI2) on the data acquired by Landsat 5 mission. UA, user's accuracy; PA, producer's accuracy; OA, overall accuracy.   Tables 13-15 describe the UA, PA, and AO of STRed index computed on the images acquired by Landsat 7, 5, and 8, respectively. STRed index performance was satisfying since the OA was higher than 80.95 for all the selected images. Moreover, UA and PA showed a satisfying value for all the data as well, regardless of the sensors and period under investigation. Indeed, their value was higher than 62.93 with the exception of UA (54.57) for the dense vegetation class extracted from the data acquired on 21 January 2002 (Landsat 7) (Table 13).  UA, PA, and AO of SwiRed index computed on the images acquired by Landsat 7, 5, and 8, respectively, are shown in Tables 16-18, respectively. SwiRed index shows satisfying outcomes as well. Indeed, the OA observed was higher than 85%, while UA and PA were on average equal to 72.56, except for the built-up areas extracted by Landsat5 on 25 August 2011 (58.21). Table 16. OA, PA, and UA obtained computing SwiRed index on the images acquired by Landsat 7 mission. UA, user's accuracy; PA, producer's accuracy; OA, overall accuracy.

Discussion
This paper proposed a new classification algorithm to automatically extract LC/LU information from Landsat satellite open data: Landsat Images Classification Algorithm (LICA). Although no classification method showed optimal performances in all conditions, the index-based method was efficient and robust in detecting LC/LU classes in a short time using satellite images provided by several sources, as highlighted by [32]. Therefore, the index classification method was selected as the benchmark approach to develop the Landsat Images Classification Algorithm, introduced in this paper. LICA integrates two new indices, namely STRed and SwiRed, obtained by combining ad hoc spectral bands in order to classify the whole study area. As shown in Section 2.2, the selected bands were chosen by examining the literature review describing the role of each spectral band [44][45][46], the performance of 82 widely spread indices (listed in Table 1), and the specific spectral signature of each class existing in the study area, revealed in the experimental area under investigation. Therefore, the former index integrated Swir, Tir, and Red bands (Equation (4)), identifying water, mining areas, and sparse and dense vegetation ( Figure 7); the latter, instead, combined Swir and Red bands (Equation (5)), distinguishing built-up areas (Figure 8). LICA was tested on 12 satellite images related to the experimental site of Siponto, an historical municipality in the Apulian Region, Southern Italy ( Figure  1). One image for each season was selected from three Landsat missions (5, 7, and 8) for a total of 12 images. The 82 conventional indices were applied to the study area as well. Among them, just three traditional indices showed quite satisfying outcomes: OSAVI (Figure 4), GOSAVI ( Figure 5), and NDBaI2 ( Figure 6). However, the first two indices (OSAVI and GOSAVI) could just distinguish water, and sparse and dense vegetation, while the third, in addition to those, also identified built-up, bare soil, and mining areas. Their outcomes are confrontable with that one obtained by other research activities. OSAVI and GOSAVI belong to the vegetation indices (VIs) category and, therefore, they are aimed at identifying vegetation class [52,53]. VIs group is composed by many indices, which must be chosen according to the environmental features since each of them is suitable for meeting a specific purpose [130]. Commonly, Vis combining visible and NIR bands show a better sensitivity in detecting green areas [130]. This paper confirmed these assumptions; indeed, both GOSAVI and OSAVI integrated visible and NIR bands. Moreover, GOSAVI showed a higher accuracy than OSAVI, thanks to the introduction of the green band, which is more sensitive to the presence and vitality of vegetation [130]. Conversely, NDBaI was proposed to discriminate different LC/LU categories even if it showed some difficulties in recognizing the bare rock areas and in distinguishing agricultural from urban areas in the zones where the urban heat phenomenon is serious [131]. Therefore, it was partially modified and NDBaI2 was introduced to improve its performance. Here, both of them were able to classify the whole study area, even if the best OA of NDBaI2 (82.59) was higher than NDBaI OA (67.93) (Table 1). Nevertheless, NDBaI2's accuracy was strongly influenced by its difficulties in distinguishing built-up areas and sparse vegetation (Tables 10-12) in all the collected images. The Automated Water Extraction Index (AWEI) was able to discriminate the different kinds of categories as well as NDBaI and NDBaI2, but its accuracy was considerably lower than NDBaI2 OA value. The worst performance was shown by Misra Yellow Vegetation Index (MYVI) [93] and Triangular Greenness Index (TGI) [112] since they were not able to discriminate any LC/LU categories in the experimental site (Table 1). MYVI was based on empirical methods without considering atmospheresoil-vegetation interactions. Therefore, it was particularly affected by soil brightness, encountering some difficulties in extracting land cover information [130]. Although TGI was proposed to assess vegetation zones, it was strongly affected by the scale and by chlorophyll content, showing promising results only applying high-resolution images as input. This resulted in its inability in identifying vegetation using medium-resolution data provided by Landsat missions [132]. Moreover, although Automated Water Extraction Index (shadow version) (AWEIsh) and Ashburn Vegetation Index (AVI) had really high overall accuracy, equal to 91.46 and 99.78, respectively, they could detect a few of LC/LU classes: The former detected water and built-up areas, the latter detected only water.
In view of their performance, NDBAI2 was chosen as a base to develop the new algorithm. Thus, NDBaI2 and the proposed algorithm were the only ones able to extract the maximum number of LC/LU classes with a high overall accuracy. Moreover, LICA went beyond NDBaI2's limitations: Both STRed and SwiRed showed a higher OA than NDBaI2, solving the issues encountered by the last one in classifying built-up areas and sparse vegetation (Tables 13-18). This was due to the introduction of R band, required to improve index performance in detecting vegetated areas since R is sensitive to the energy absorbed by chlorophyll [52]. Moreover, SWIR and TIR1 bands were also combined to distinguish bare soil and built-up areas [53]. This resulted in an optimal OA of STRed and SwiRed, equal to 94.71% and 97.76%, respectively. Besides maximizing the number of categories to be detected and improving classification accuracy, LICA was designed in order to be applied on all Landsat missions, equipped with different sensors, so that multisensors, multitemporal, and multiseason analyses, which are essential in environmental monitoring and planning management, could be performed. Moreover, users can apply the whole algorithm or just one of the two proposed indices, according its needs.
To automatize LC/LU extraction, the algorithm was implemented in GEE environment, a new platform recently designed by Google (https://earthengine.google.com/). Thanks to its parallel processing capacity, already shown in previous research activities [20], LIC algorithms can be run in a few minutes, even if computation times increase with the amount of data to be handled. Therefore, using GEE allowed overcoming desktop system limitations due to excessive processing time needed to process geospatial big data. This paper confirmed the great potentiality of the GEE platform in processing geospatial big data, as already shown in previous research works [18][19][20]113].

Conclusions
In this study, an automated algorithm for extracting land cover information from multitemporal and multisensors open data in the GEE platform was introduced. The procedure did not need any external training datasets, which are time consuming (collection time must be considered) and may be affected by human errors. On the contrary, LICA used the integration of two novel indices (STRed and SwiRed) which allowed us to analyze land covers from Landsat images, maximizing the number of classes to be extracted and increasing classification accuracy, compared to the conventional indices commonly applied in literature. Landsat images were selected to test LICA in order to exploit the huge amount of open data available from 1972 and ensuring its reliability in multitemporal and multisensors analyses in order to provide information that could be used to perform land cover change analyses, which are essential to guide future planning strategies.
All computational steps were implemented in the GEE cloud computing platform, thereby avoiding the necessity of excessive desktop processing power to handle geospatial big data and automating the whole procedure. Therefore, the integration of the LIC algorithm and GEE environment allowed us to quickly extract accurate land cover information. Thus, adopting the proposed method helps to provide more contemporary information while also reducing costs, acquisition, and processing times.