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

Abstract

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.

1. 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 European Commission’s Joint Research Center (JRC) in 2000, the Globcover product realized by the European Space Agency, and the Moderate-resolution Imaging Spectroradiometer (MODIS) Collection 5 Land Cover database show a resolution of about 1 km at the equator (larger at higher latitudes) [4,5], 300 m [6,7], and 500 m [8], respectively. The 300-m global Climate Change Initiative Land Cover (CCI-LC) maps covering the period from 1992 to 2015 were utilized to assess the quality of the unchanged training sample pixels in five time periods. The Global Land Cover Characterization database (GLCC), produced by the effort between the U.S. Geological Survey (USGS), University of Nebraska Lincoln (UNL), and the JRC [9], shows a resolution of 1 km as well. Recent research activities have shown that, due to their coarse resolution, most of the listed datasets are not reliable over urban and agricultural areas since they show substantial disagreement with each other and with national statistics [1,10,11]. A much smaller number of high-resolution LC/LU maps, based on the available Landsat data, were generated at large scale and at various timescales as well. Nevertheless, such maps were produced for forestry purposes and, consequently, they do not report LC/LU information [12]. Similarly, Landsat data were also applied to provide contemporary data on human population distributions in Africa, Asia, and the Americas (WorldPop project) [13]. Only three Landsat-based global land cover maps are currently available: Finer Resolution Observation and Monitoring of Global Land Cover (FROM-GLC) by [14], GlobeLand30 by [15], and Normalized Urban Areas Composite Index (NUACI) derived maps by [16]. FROM-GLC and GlobeLand30 provide LC/LU information for the years of 2000 and 2010. Conversely, [16] provided LC/LU maps for the period 1990-2010 at five-year intervals. The situation changes at a continental, national, and regional scale, where Landsat and Sentinel images were widely used in many applications [17,18,19,20,21,22,23,24].

Landsat archive has provided a limitless well of information since 1972, freely available and accessible and, consequently, suitable for describing Earth surface features. However, as underlined by [18], Landsat satellite implies three main challenges:

(1)

Dealing with a low number of useful images: Just satellite images with minimal cloud cover are acceptable. Thus, the amount of adequate data depends on the weather conditions of the experimental site. Consequently, areas characterized by a high rainfall, such as tropical and subtropical regions, show a lower number of available adequate images;

(2)

Identifying an efficient platform suitable for large image data processing; and

(3)

Developing adequate image classification methods with satisfactory performance.

The introduction of Google Earth Engine (GEE) (https://earthengine.google.org), a cloud processing platform designed and developed over the last years by Google, offered a large number of available tools to face the first two issues [25]). As emphasized by [26], GEE integrates a data catalogue, continuously updated and composed of publicly available geospatial datasets, which may be consulted by users through the application programming interface (API). Therefore, operators can handle hundreds of data sources simultaneously and detect their quality and usefulness to meet their purposes. Moreover, nearly 6000 scenes belonging to active missions are integrated into GEE catalogue daily. As well as private data, such scenes can be processed by applying a set of complex and advanced algorithms implemented in GEE environment, exploiting its excellent computational power. Unlike desktop software, it involves many processors in running custom algorithms, speeding the process up considerably, and deleting the problems linked to the storage, the processing, and the analysis of a large volume of geospatial data [25]. For example, [12] tracked forest cover changed over a period of 12 years (2000–2012) at global scale by analyzing 654,178 Landsat 7 scenes (707 terabytes) on GEE platform. The milestone was achieved in 100 h, while a standard desktop computer would have taken nearly 1,000,000 h to meet the same needs. In addition, GEE is more flexible than the software commonly applied to process geospatial data, such as Environment for Visualizing Images (ENVI) and Earth Resources Data Analysis System (ERDAS) Imagine, since users can implement their own custom codes. Although its potentialities are enormous, GEE is still in development and, consequently, many existing algorithms have not been programmed and integrated into the platform yet [25]. Moreover, its great versatility sets it as the standard to deal with the third challenge implied by Landsat satellite and to carry out several classification approaches introduced over the years.

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 index-based 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, built-up 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.

2. Materials and Methods

2.1. 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].

2.2. 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. Main commonly used classification indices listed in alphabetical order. Indices in bold show the best performance. LC/LU column describes the land cover/land use classes detected from each index. OA column reports the best overall accuracy of each index. LC/LU, land cover/land use; OA: overall accuracy; W, water; DV, dense vegetation; SV, sparse vegetation; MA, mining areas; BS, bare soil; BUA: built-up area; *, water mask is required; -:, no classes were detected.

Spectral Index	Citation	LC/LU	OA (%)
Aerosol Free Vegetation Index version 1.6 (AFRI1.6)	[60]	DV, SV	72.24
Aerosol Free Vegetation Index version 2.1 (AFRI2.1)	[60]	DV, SV	86.02
Atmospherically resistant vegetation index (ARVI)	[61]	W, DV, SV, BUA, BS	59.97
Adjusted Soil Brightness Index (ASBI) *	[62]	DV, SV	66.70
Ashburn Vegetation Index (AVI)	[63]	W	99.78
Automated Water Extraction Index (AWEI)	[64]	W, DV, SV, BUA, MA, BS	68.04
Automated Water Extraction Index (shadow version) (AWEIsh)	[64]	W, BUA	91.46
Build-area extraction index (BAEI) *	[65]	DV, SV, BUA	63.60
Biophysical Composition Index (BCI)	[66]	W, DV, SV	68.23
Built-up Land Features Extraction Index (BLFEI)	[67]	W, DV, SV, BUA, BS	72.03
Bare Soil Index (BSI) *	[68]	DV, SV	73.62
Built-up land (BUI)	[69]	W, DV, SV	69.81
Combinational Biophysical Composition Index (CBCI)	[70]	DV, SV	67.22
Green Chlorophyll Index (CI)	[71]	W, DV, SV	68.40
Davies-Bouldin index (DBI)	[72]	W, DV, SV, BUA, BS	70.59
Dry Bare-Soil Index (DBSI) *	[73]	DV, SV	68.47
Simple Difference Indices (DVI)	[74]	W, DV, SV	69.85
Enhanced Built-up and Bareness Index (EBBI)	[75]	W, DV, SV	64.93
Enhanced Normalized Difference Impervious Surfaces Index (ENDISI)	[76]	DV, SV BUA, MA	67.55
Enhanced Vegetation Index (EVI)	[77]	W, DV, SV, BUA, BS	58.59
Green Atmospherically Resistant Vegetation Index (GARI)	[78]	W, DV, SV, BUA, BS	69.78
“Ghost cities” Index (GCI)	[79]	W, DV, SV, BUA, BS	71.26
Green Difference Vegetation Index (GDVI)	[80]	W, DV, SV	70.59
Global Environment Monitoring Index (GEMI)	[81]	W, DV, SV	67.74
Green leaf index (GLI)	[82]	DV, SV	66.70
Green Normalized Difference Vegetation Index (GNDVI)	[78]	W, DV, SV, BUA, BS	72.48
Green Optimized Soil Adjusted Vegetation Index (GOSAVI)	[53]	W, DV, SV	89.89
Green-Red Vegetation Index (GRVI)	[83]	W, DV, SV, BUA, BS	71.26
Green Soil Adjusted Vegetation Index (GSAVI)	[53]	W, DV, SV, BUA, BS	73.91
Green Vegetation Index (GVI) *	[84]	DV, SV, BUA	57.30
Built-up Index (IBI)	[85]	DV, SV	74.75
Infrared Percentage Vegetation Index (IPVI)	[86]	W, DV, SV, BUA, BS	69.10
Modified Bare Soil Index (MBSI)	[70]	W, DV, SV	73.22
Modified Chlorophyll Absorption Ratio Index1 (MCARI1)	[87]	DV, SV	64.28
Modified Chlorophyll Absorption Ratio Index (MCARI2)	[87]	W, DV, SV, BUA, BS	82.24
MERIS Global Vegetation Index (MGVI)	[88]	W, DV, SV	76.88
Modification of Normalized Difference Snow Index (MNDSI)	[89]	W, MA	76.55
Modification of normalized difference water index (MNDWI)	[90]	W, BUA	74.62
Modified Nonlinear Vegetation Index (MNLI)	[91]	W, DV, SV	77.40
Modified Soil Adjusted Vegetation Index 2 (MSAVI2)	[92]	W, DV, SV, BUA, BS	83.30
Misra Soil Brightness Index (MSBI)	[93]	W, DV, SV, BUA, BS	78.56
Modified Simple Ratio (MSR)	[94]	W, DV, SV, BUA, BS	67.03
Misra Yellow Vegetation Index (MYVI)	[93]	-	-
New Built-up Index (NBI) *	[95]	DV, SV, BUA, MA, BS	71.46
Normalized Difference Bare Land Index (NBLI) *	[96]	DV, SV, BUA, MA, BS	75.51
New Built-up Index (NBUI)	[97]	W, DV, SV	76.39
Normalized Canopy Index (NCI)	[98]	W, BUA	78.34
Normalized Difference Bareness Index (NDBaI)	[54]	W, DV, SV, BUA, MA, BS	67.93
Normalized Difference Bareness Index (version 2) (NDBaI2)	[54]	W, DV, SV, BUA, MA, BS	82.59
Normalized Difference Built-up Index (NDBI)	[99]	DV, SV	71.14
Normalized Difference Impervious Surface Index (NDISI)	[100]	W, MA	97.60
Normalized Difference Moisture Index (NDMI) *	[101]	DV, SV	73.47
Normalized Difference Tillage Index (NDTI) *	[102]	DV, SV	71.57
Normalized Difference Vegetation Index (NDVI)	[56]	W, DV, SV, BUA, BS	73.24
Normalized Difference Water Index (NDWI)	[103]	W, DV, SV, BUA, BS	73.54
Non-Linear Index (NLI)	[104]	W, DV, SV	76.63
Optimized Soil Adjusted Vegetation Index (OSAVI)	[52]	W, DV, SV	88.84
Renormalized Difference Vegetation Index (RDVI)	[105]	W, DV, SV	77.34
Ratio Vegetation Index (RVI)	[106]	W, DV, SV, BUA, BS	72.30
Soil-Adjusted Vegetation Index (SAVI)	[107]	W, DV, SV	72.04
Soil Brightness Index (SBI)	[108]	W, BUA, MA	80.27
Specific Leaf Area Vegetation Index (SLAVI)	[109]	W, DV, SV	83.56
Simple Ratio (SR)	[110]	W, DV, SV	68.93
Transformed difference vegetation index (TDVI)	[111]	W	99.81
Triangular Greenness Index (TGI)	[112]	-	-
Triangular Vegetation Index (TVI)	[113]	W, DV, SV	74.15
Urban Index (UI)	[114]	BUA	76.66
Visible Atmospherically Resistant Index (VARI)	[115]	W, DV, SV	68.34
Visible-Band Difference Vegetation Index (VDVI)	[116]	DV, SV	66.70
Vegetation Index of Biotic Integrity (VIBI)	[117]	DV, SV	66.57
Wide Dynamic Range Vegetation Index (WDRVI)	[118]	W, DV, SV	78.87
Water index 2015 (WI2015)	[119]	W	99.81
Worldview Improved Vegetative Index (WV-VI)	[120]	W, DV, SV, BUA, BS	75.47
Yellow Stuff Index (YVI) *	[121]	DV, SV	66.70

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

$$OSAVI=\frac{1.16\times \left(NIR-R\right)}{NIR+R+0.16}$$

(1)

$$GOSAVI=\frac{NIR-G}{NIR+G+0.16}$$

(2)

$$NDBaI2=\frac{SWIR1-TIR1}{SWIR1+TIR1}$$

(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).

$$STRed index=\frac{SWIR1+R-TIR1}{WIR1+R+TIR1}$$

(4)

$$SwiRed index=\frac{SWIR1-R}{SWIR1+R}$$

(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. Range value of LICA to extract the different land cover classes.

Land Cover Class	Range value (SwiRed)
Built- up areas	0 < value < 0.22
Land Cover Class	Range value (STRed)
Water	value < −0.5
Dense vegetation	−0.05 < value < −0.07
Sparse vegetation	0.07 < value < 0.00
Mining areas	Value > 0.45

) and the resultant maps were merged. Figure 3 describes the suggested workflow to be set.

2.3. 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. Selected Landsat data description. ETM+, enhanced thematic mapper; TM, thematic mapper; OLI-TIRS, operational land imager - thermal infrared.

ID	Landsat Satellite Mission	Sensor	Landsat Images	Acquisition Date	Average Cloud Cover (%)
1	Landsat 7	ETM+	LE07_L1TP_188031_20020121_20170213	21 January 2002	4
2			LE07_L1TP_188031_20020801_20170213	01 August 2002	6
3			LE07_L1TP_189031_2002127_20170128	27 October 2002	1
4			LE07_L1TP_188031_20030414_20170126	14 April 2003	4
1	Landsat 5	TM	LT05_L1TP_188031_20110207_20161010	07 February 2011	1
2			LT05_L1TP_188031_20110327_20161209	27 March 2011	16
3			LT05_L1TP_189031_20110825_20161008	25 August 2011	0
4			LT05_L1TP_188031_20111005_20161005	05 October 2011	1
1	Landsat 8	OLI-TIRS	LC08_L1TP_188031_20171208_20171223	08 December 2017	1.69
2			LC08_L1TP_189031_20180812_20180815	12 August 2018	8.1
3			LC08_L1TP_188031_20180922_20180928	22 September 2018	2.41
4			LC08_L1TP_188031_20190925_201911017	17 March 2019	19.46

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

2.4. 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).

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

3. Results

3.1. Classification Results

This section is dedicated to the classification procedure outcomes obtained through the application of indices consolidated in literature (Figure 4, Figure 5 and Figure 6) and the proposed LIC algorithm (Figure 7 and Figure 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 (Figure 7, Figure 8 and Figure 9); the latter aimed at identifying built-up areas (Figure 10, Figure 11 and Figure 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.

3.2. Accuracy Assessment

Table 4. OA, PA, and UA obtained through the application of Optimized Soil Adjusted Vegetation Index (OSAVI) on the data acquired by Landsat 7 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L7—21 January 2002		L7—14 April 2003		L7—01 August 2002		L7—27 October 2002
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	95.79	100.00	95.04	100.00	90.11	98.99	95.60	99.24
Dense Vegetation	79.57	87.40	72.38	86.21	70.97	98.05	95.56	45.84
Sparse Vegetation	73.89	89.78	90.55	81.91	25.14	33.96	49.30	67.03
Not classified	100.00	78.38	100.00	92.28	91.13	62.43	79.47	98.17
Mining Areas	/	/	/	/	/	/	/	/
Bare Soil	/	/	/	/	/	/	/	/
Built-up areas	/	/	/	/	/	/	/	/
OA (%)	87.83		88.84		74.02		78.23

,

Table 5. OA, PA, and UA obtained through the application of Optimized Soil Adjusted Vegetation Index (OSAVI) on the data acquired by Landsat 5 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L5—07 February 2011		L5—27 March 2011		L5—25 August 2011		L5—05 October 2011
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	95.05	93.35	89.73	100.00	90.30	99.05	66.30	96.02
Dense Vegetation	64.16	53.27	66.49	78.44	48.54	93.55	80.54	69.00
Sparse Vegetation	35.36	58.67	76.25	81.06	43.05	51.90	67.92	80.36
Not classified	92.21	70.75	100.00	71.88	81.90	40.89	93.83	72.77
Mining Areas	/	/	/	/	/	/	/	/
Bare Soil	/	/	/	/	/	/	/	/
Built-up areas	/	/	/	/	/	/	/	/
OA (%)	68.71		84.84		72.44		77.98

,

Table 6. OA, PA, and UA obtained through the application of Optimized Soil Adjusted Vegetation Index (OSAVI) on the data acquired by Landsat 8 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L8—08 December 2008		L8—12 August 2018		L8—22 September 2018		L8—17 March 2019
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	99.13	100.00	98.19	99.85	89.01	99.53	97.30	100.00
Dense Vegetation	76.93	75.39	59.09	91.95	86.49	96.60	80.85	73.47
Sparse Vegetation	22.91	52.00	27.34	51.27	81.74	82.82	55.04	71.57
Not classified	99.46	70.60	99.81	65.33	93.98	83.79	100.00	79.54
Mining Areas	/	/	/	/	/	/	/	/
Bare Soil	/	/	/	/	/	/	/	/
Built-up areas	/	/	/	/	/	/	/	/
OA (%)	81.41		81.00		88.91		83.56

,

Table 7. OA, PA, and UA obtained through the application of Green Optimized Soil Adjusted Vegetation Index (GOSAVI) on the data acquired by Landsat 7 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L7—21 January 2002		L7—14 April 2003		L7—01 August 2002		L7—27 October 2002
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	96.34	100.00	96.70	100.00	86.26	100.00	96.15	97.22
Dense Vegetation	60.57	88.48	70.22	91.18	70.97	93.12	86.69	44.79
Sparse Vegetation	65.78	83.83	94.67	76.05	18.23	39.05	52.09	71.39
Not classified	100.00	70.57	89.16	95.52	97.21	59.50	83.09	97.50
Mining Areas	/	/	/	/	/	/	/	/
Bare Soil	/	/	/	/	/	/	/	/
Built-up areas	/	/	/	/	/	/	/	/
OA (%)	82.86		87.74		73.57		79.12

,

Table 8. OA, PA, and UA obtained through the application of Green Optimized Soil Adjusted Vegetation Index (GOSAVI) on the data acquired by Landsat 5 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L5—07 February 2011		L5—March 27 2011		L5—25 August 2011		L5—05 October 2011
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	94.32	100.00	95.73	100.00	93.56	58.46	74.91	97.61
Dense Vegetation	60.04	43.51	63.98	88.15	52.30	89.29	81.71	54.40
Sparse Vegetation	25.03	35.39	89.31	81.74	43.05	40.25	59.62	85.87
Not classified	59.51	52.17	99.58	100.00	90.26	43.90	95.41	77.60
Mining Areas	/	/	/	/	/	/	/	/
Bare Soil	/	/	/	/	/	/	/	/
Built-up areas	/	/	/	/	/	/	/	/
OA (%)	55.22		89.89		76.37		78.82

,

Table 9. OA, PA, and UA obtained through the application of Green Optimized Soil Adjusted Vegetation Index (GOSAVI) on the data acquired by Landsat 8 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L8—08 December 2008		L8—12 August 2018		L8—22 September 2018		L8—17 March 2019
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	99.49	100.00	98.19	100.00	90.84	98.41	99.81	100.00
Dense Vegetation	75.62	75.29	59.72	89.02	84.80	89.64	92.88	63.74
Sparse Vegetation	40.50	24.23	25.20	43.56	75.00	87.34	17.94	61.86
Not classified	99.46	66.95	100.00	66.31	97.95	84.62	100.00	99.74
Mining Areas	/	/	/	/	/	/	/	/
Bare Soil	/	/	/	/	/	/	/	/
Built-up areas	/	/	/	/	/	/	/	/
OA (%)	79.32		80.89		89.15		80.35

,

Table 10. OA, PA, and UA obtained through the application of Normalized Difference Bareness Index (version 2) (NDBIaI2) on the data acquired by Landsat 7 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L7—21 January 2002		L7—14 April 2003		L7—01 August 2002		L7—27 October 2002
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	95.24	100.00	97.07	100.00	87.55	100.00	95.24	100.00
Dense Vegetation	69.89	54.32	81.48	95.14	41.94	50.61	89.11	79.78
Sparse Vegetation	19.62	51.75	52.36	75.66	69.34	60.19	26.84	27.27
Not classified	/	/	/	/	/	/	/	/
Mining Areas	30.96	92.55	77.58	94.78	90.20	95.04	64.77	95.25
Bare Soil	94.35	63.88	87.17	62.54	91.33	79.71	69.31	64.17
Built-up areas	27.18	31.98	41.04	37.04	31.79	48.44	46.15	46.88
OA (%)	66.89		76.10		76.23		66.61

,

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.

	L5—07 February 2011		L5—27 March 2011		L5—25 August 2011		L5—05 October 2011
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	96.89	100.00	90.59	100.00	82.49	100.00	84.25	100.00
Dense Vegetation	86.56	53.91	74.73	73.42	62.97	45.47	53.70	54.55
Sparse Vegetation	54.76	84.42	38.24	65.05	78.63	78.63	48.11	57.41
Not classified	/	/	/	/	/	/	/	/
Mining Areas	56.86	99.32	86.27	98.65	84.71	99.08	82.56	97.89
Bare Soil	93.81	89.17	92.13	59.38	92.50	83.88	79.33	58.56
Built-up areas	28.21	26.19	16.92	31.43	42.82	38.66	34.87	45.26
OA (%)	78.95		75.07		80.28		61.55

and

Table 12. OA, PA, and UA obtained through the application of Normalized Difference Bareness Index (version 2) (NDBIaI2) on the data acquired by Landsat 8 mission. UA: UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L8—08 December 2008		L8—12 August 2018		L8—22 September 2018		L8—17 March 2019
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	98.55	100.00	100.00	99.93	93.77	99.03	96.72	100.00
Dense Vegetation	91.53	57.10	32.45	70.41	39.86	77.12	98.89	80.93
Sparse Vegetation	27.48	50.30	69.59	62.02	36.96	48.60	37.59	71.83
Not classified	/	/	/	/	/	/	/	/
Mining Areas	41.96	100.00	81.18	99.04	74.02	99.52	79.92	100.00
Bare Soil	93.72	87.34	93.36	80.67	88.88	32.16	30.14	39.25
Built-up areas	29.74	39.27	26.41	32.54	20.00	38.31	70.19	23.10
OA (%)	80.20		82.59		66.31		72.04

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 13. OA, PA, and UA obtained computing STRed on the images acquired by Landsat 7 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L7—21 January 2002		L7—14 April 2003		L7—01 August 2002		L7—27 October 2002
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	95.79	97.39	98.80	100.00	98.35	100.00	99.15	96.88
Sparse Vegetation	85.75	69.37	72.56	92.81	62.98	74.51	92.95	93.78
Dense Vegetation	76.76	54.57	98.66	89.84	75.65	69.52	84.08	76.41
Mining areas	81.57	100.00	75.11	99.41	97.65	88.25	46.44	97.18
Not classified	85.67	95.78	99.73	80.62	89.66	72.90	61.93	83.63
OA (%)	86.49		94.30		80.95		87.73

,

Table 14. OA, PA, and UA obtained computing STRed on the images acquired by Landsat 5 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L5—07 February 2011		L5—27 March 2011		L5—25 August 2011		L5—05 October 2011
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	99.27	98.91	98.72	100.00	94.86	99.62	97.38	100.00
Sparse Vegetation	81.40	90.78	91.45	95.65	89.69	94.63	68.64	55.10
Dense Vegetation	80.11	80.54	98.31	91.26	84.94	78.99	96.04	82.20
Mining areas	69.80	99.44	98.67	99.11	87.45	100.00	94.22	99.53
Not classified	99.32	80.87	99.31	97.17	99.47	89.61	86.83	78.42
OA (%)	87.88		97.76		93.20		85.83

and

Table 15. OA, PA, and UA obtained accuracy computing STRed on the images acquired by Landsat 8 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L8—08 December 2008		L8—12 August 2018		L8—22 September 2018		L8—17 March 2019
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Water	98.63	97.92	97.68	100.00	98.65	100.00	99.43	100.00
Sparse Vegetation	94.21	97.93	89.59	87.13	69.29	72.73	93.53	98.64
Dense Vegetation	86.28	92.92	73.01	74.72	73.74	86.13	99.17	93.93
Mining areas	84.51	96.53	82.75	98.60	61.56	100.00	98.22	97.79
Not classified	99.55	56.08	99.53	79.52	100.00	53.04	98.89	99.19
OA (%)	93.33		87.93		85.08		98.71

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

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.

	L7—21 January 2002		L7—14 April 2003		L7—01 August 2002		L7—27 October 2002
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Built-up areas	72.82	82.75	90.45	73.77	74.10	74.10	80.00	83.88
Not classified	98.27	96.05	93.88	89.07	97.88	97.32	97.55	97.20
OA (%)	91.31		89.28		92.77		94.30

,

Table 17. OA, PA, and UA obtained computing SwiRed index on the images acquired by Landsat 5 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L5—07 February 2011		L5—27 March 2011		L5—25 August 2011		L5—05 October 2011
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Built-up areas	77.44	71.79	76.11	72.56	58.21	87.04	86.11	79.81
Not classified	93.90	93.97	97.43	90.50	98.07	93.09	97.37	97.43
OA (%)	91.03		89.80		94.00		95.56

and

Table 18. OA, PA, and UA obtained computing SwiRed index on the images acquired by Landsat 8 mission. UA, user’s accuracy; PA, producer’s accuracy; OA, overall accuracy.

	L8—08 December 2008		L8—12 August 2018		L8—22 September 2018		L8—17 March 2019
Land Cover Class	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)	PA (%)	UA (%)
Built-up areas	76.92	84.73	72.56	77.70	85.00	72.26	74.78	80.51
Not classified	97.69	98.07	78.81	96.01	97.49	93.73	97.73	92.74
OA (%)	94.71		85.31		91.40		91.20

, 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).

4. 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 (Table 10, Table 11 and Table 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 atmosphere-soil-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 (Table 13, Table 14, Table 15, Table 16, Table 17 and Table 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].

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