The study area is the Municipal territory of Acquaviva delle Fonti (40°53′N, 16°50′E). It is located in the province of Bari (Northern Puglia, Southeastern Italy) (Figure 2) and covers approximately 117 km².

The terrain is flat and lies in the hilly life zone with an elevation between 189 and 441 m above sea level. The climate is typically Mediterranean with warm, dry summers and mild, moist winters. The average annual temperature is 15.3°C and the average rainfall is 645 mm year⁻¹.

Eight land cover classes were identified after site inspections during spring of 2021. Cereals, mainly winter wheat, and vineyard under plastic films, constitute, respectively, the winter wheat and plastic films LULC classes. Other crops (crops in the rest of the article) was the class with the greatest extension and fragmentation. It mainly hosts legume and vegetable cultivations, olive and orchards. To a lesser degree, strips of oak (Quercus pubescens Willd. and Q. trojana Webb) and coniferous woods (Pinus halepensis Mill.) are found as forest. The municipality of Acquaviva delle Fonti and its industrial site fall within the investigated area as built-up. The mining activity (bare soil as land cover class) is limited to a site close to the urban center. A highway road also crosses the territory, as well as numerous photovoltaic fields.

After identification of the land cover types which dominate the territory (Table 1, first column), a reference data collection (ground truth or training samples and validation samples) was created. The choice of reference dataset is an important step addressed at the best classification accuracy. The training dataset allows the selected classifier algorithm to learn the relationship between the pixel values of the various bands of composite/normalized indices/radar images and the land cover classes, with the aim of attaining optimal accuracy. The validation dataset allows to test the accuracy of the classifier.

For these purposes, the training and validation dataset was designed combining imagery observation and field investigation. The reference dataset was determined using a random sampling approach, selecting the number of samples proportional to land cover classes surface (Table 1, second column), and taking into account the purity of samples (without mixed pixels).

To select the most homogenous pixel representative of every class and to reduce the effect of spatial autocorrelation, the reference samples were uploaded via GEE. Furthermore, data balancing technique, such as the augmentation of minority classes, was used to minimize class imbalance in the training dataset. Accurate visual analysis of Sentinel-2 RGB (2021) and Google Earth high-resolution (2018) images were employed to create reference samples. The assigned labeled classes were cross-checked by visually inspecting Google Earth imagery. The Google Earth images were used as a georeferenced base map to collect reference samples and to address field inspections, in case of inconsistencies between the visual information provided by the different images. In instances of uncertainty, class assignment was confirmed by checking normalized index values derived from optical images or visually inspecting ESA WorlCover maps (GEE image collection: ESA WorldCover v100 and v200; 2020, 2021). Moreover, samples with persistent high uncertainty underwent removal or were subject to field verification carried out in late spring 2021.

A total of 2492 reference samples were iteratively delineated for the period under consideration (Table 1). 80% of the randomly selected reference samples were used as training data, and the remaining 20% were used as validation data. Different data splitting approaches can be used to ensure a separation between the training and test data, including the holdout, the bootstrapping, and the cross-validation procedures. To avoid differences in accuracies, the same samples and split process were used to train and validate each dataset.

Using multispectral/multi-source remote sensing data over large areas requires considerable technical expertise and high computational complexity (Shaharum et al., 2019). Those challenges can be effectively addressed by Google Earth Engine (http://earthengine.google.com), a cloud-based platform providing access to a multi-petabyte catalogue of Earth observation data, analysis algorithms and many other digital products (Padarian et al., 2015) that are ready-to-use with an explorer web app. It was designed to take advantage of Google’s computational infrastructure for planetary-scale access, storage, monitoring, analysis and visualization of geospatial data. An Application Programming Interface (API) enables users to create custom JavaScript and Python algorithms, while an associated web-based Interactive Development Environment (IDE), the Earth Engine Code Editor, provides high-speed parallel processing for prototyping, analysis and visualization of large-scale geospatial data (Gorelick et al., 2017; Tamiminia et al., 2020).

The GEE public data repository includes scientific datasets and several decades of historical images from multiple satellites. For this work, we used remote sensing data for which cloud computing and archived processed data in GEE had advantages for long time-series mapping, such as monitoring land type changes and detecting the extent of human settlement (Johansen et al., 2015). Moreover, machine-learning algorithms integrated into the GEE API allow short-time information extraction from satellite data, avoiding local memory consumption and intensive data transfers (Huang et al., 2017).

The remote sensing dataset consisted of 63 radar images and 16 optical images (Table 2) selected over the study area for the vegetative period 01/01/2021–31/05/2021. We used data from two optical sensors (Sentinel-2 MSI and Landsat 8 OLI) both to increase the number of optical images and to exploit the different bands resolution.

The GEE cloud platform provides online access to a centralized catalogue of georeferenced and atmospherically corrected satellite multispectral images (Bunting et al., 2019), like the USGS and Copernicus archives that include Landsat 8 OLI/TIRS and Sentinel-2 (level 2A) data, respectively. We chose the atmospherically corrected surface reflectance scenes with the aim of minimizing the effects arising from the differences between satellite sensors (Landsat 8 and Sentinel-2) or the different acquisition dates. In addition, from the Copernicus collection, we selected the Sentinel-1 Synthetic Aperture Radar (SAR) (Table 2).

Several normalized spectral indices (NI(s)) have been proposed in the literature to extract the land cover types from remote sensing imagery. Data normalization has the advantage of reducing image noise when multi-sensor and multi-temporal images are used (Yang and Lo, 2002; Zurqani et al., 2018). Furthermore, these spectral indices could improve image classification accuracy, thereby emphasizing the detection of vegetation reflectance signatures or reducing hill shade and building shadows. Supplementary Figure S1 shows a Sentinel-2 derived NDVI time series for the eight land cover classes detected in the study area. Winter wheat had one NDVI peak in 2021 indicating one harvest per year. Deciduous oaks (Quercus trojana Webb. and Q. pubescens Willd.) have a growing stage in spring after the new foliage. Coniferous trees (Pinus halepensis Mill. and P. pinea L.) show a somewhat varying index, such as the urban dwelling. This illustrated the importance of using seasonal information to discriminate vegetation covers from infrastructures.

Some of the standardized indices are sensitive to different factors related to seasonality, study area location, and image resolution. However, given their simplicity and easy formulation, such methods have been widely used for LULC monitoring and mapping.

In this study, the normalized indices calculated from the composite bands of Sentinel-2 and Landsat 8 imagery, and Sentinel-1 imagery are reported in Table 3:

The Normalized Difference Vegetation Index (NDVI) (Rouse et al., 1974) is an index of the vegetation greenness that derives from measurements of the optical reflectance of sunlight in the red and near-infrared wavelengths. It is not a physical property of the biome’s cover, but its very simple implementation makes it widely used for ecosystem monitoring (https://land.copernicus.eu/global/products/ndvi). NDVI is suitable for estimating leaves’ vigor during the early vegetative stages.

The Green Normalized Difference Vegetation Index (GNDVI) (Gitelson et al., 1996) is more sensitive than NDVI to the proportion of chlorophyll absorbed radiation since the Red band is replaced in the calculation by the Green band. It is used to assess photosynthetic activity, water content, and nitrogen concentration in more advanced development stages of dense plant canopies (Vizzari, 2022).

Most of the normalized vegetation indices confuse built-up with bare soil surfaces (Valdiviezo et al., 2018); to reduce such limitations, a built-up index was used in this work.

The Normalized Difference Built-up Index (NDBI) (Zha et al., 2003) uses increased reflectance values from SWIR and NIR bands to highlight manufactured built-up areas. In thermal infrared bands, the higher emissivity and albedo of built-up structures, compared to bare soil, waterbodies, and vegetation, is best detected by this index (Ali and Nayyar, 2021). It is based on separating the built-up area from the background and mitigating the effects of terrain illumination as well as atmospheric effects. According to this index, the greater the value of a pixel in the derived image, the higher the possibility of the pixel being a built-up surface (He et al., 2010).

To analyze the multispectral data combinations on the land cover classification process, we aggregated optical data, radar data and spectral normalized indices into five different input data and combined these data in thirteen different datasets, seven without NI(s) and six with NI(s) (Table 4).

The workflow of the methodological approach adopted in this study for land cover classification is reported in the Graphical abstract. First, the Copernicus (i.e., Sentinel-2 and Sentinel-1) and Landsat 8 image collections comprising the Region of Interest (ROI) were selected via GEE. The winter-spring period of 2021 was chosen to distinguish winter wheat from bare soil land cover. The selected collection data were then processed with GEE to manually detect the main land types in the study area and generate training and testing data. The Random Forest (RF) classifier was then used to generate the LULC map. Lastly, a validation step and an accuracy assessment were performed before the digital restitution of the final land cover map for the study period. The detailed methodology is described in the following paragraphs.

The RF algorithm performance in correctly classifying a random set of features (sample points) was measured by the accuracy of the fourteen classification maps produced. Therefore, once a classification was run, the accuracy of the process was determined. As explained in Section 2.2, each of the land cover map produced was then validated using a dataset obtained from 20% of the total reference locations, via GEE. To ensure the separation between train and test data, we used the cross-validation procedure, in which each sample is included exactly once in the test set and each sample in the test is not used to train the classifier. The optimal tuning parameters percent distribution was determined based on multiple percent distribution of training data to evaluate the effects of different parameter configurations on the performance of the classifier. The definitive classification dataset created for this purpose includes samples randomly selected for each land cover type. The cross-validation comparison produces a table, the confusion matrix, which contains the number of pixels classified correctly or incorrectly in each class. The matrix allows the computation of several accuracy metrics: Overall Accuracy (OA), measuring the proportion of correctly classified pixels; Producer’s accuracy (PA), evaluating how well reference (ground-truth) samples were classified; User’s accuracy (UA), assessessing classification reliability from a user perspective. The goal is to obtain not just the higher OA, but also a good balance of PA and UA, ensuring that classification performance was not biased toward dominant land cover types. To depict the balance of PA and UA, the F-score was calculated as a harmonic mean of PA and UA (Sokolova et al., 2006).

The McNemar’s statistical test (McNemar, 1947) was applied as an additional measure to compare the significance of difference between classification results, with p-values lower than 0.05 regarded as significant.

The kappa coefficient (K) was deliberately excluded because its use was recently discouraged in assessing land cover classification accuracy (Foody, 2020).

Lastly a variable importance analysis was performed for all spectral bands and indices, providing insights into which features contributed most to the classification model. The result paragraph are reported in Section 3.3.

Based on a dataset collection of 1,980 training samples and 512 validation samples, for 2,492 reference samples (Table 1), the accuracy of classifications was characterized by the Overall Accuracy (OA) always above 0.7773, except for S-1 data where the accuracy score was 0.613 (Figure 4; Table 5). The numeric values of Overall Accuracies index reported in Table 5, whose bar chart is displayed in Figure 4, provide a significant insight into the overall contribution of multi-sensor imagery to the LULC classification quality. Also, Figure 4 shows the F-score graphic bars.

OA increases from radar to optical images and with the decreasing of the image pixel size (from Landsat OLI to Sentinel-2). For classifications obtained using spectral indices, the combination of the satellite imagery slightly increases the statistical metrics, with the exception of S-2, which shows an inverse trend, decreasing OA with NI(s) (Table 5; Figure 4). The dataset that combined the complete satellite imagery (S-1/S-2/L8) and NI(s) computation, achieved the highest accuracy with OA = 0.8964 (Table 5a) and F-score = 0.7806, followed by the dataset combining the Copernicus sensor data with NI(s), where OA = 0.8828 (Table 5a) and F-score = 0.7768. In the investigated period, the cumulative number of images used demonstrated that a combination of Landsat OLI, Sentinel-1, and Sentinel-2 performed better than a single sensor. The accuracy increases with the increasing number of bands used. The results illustrated in Figure 4 and Tables 5a, b show that Sentinel-2 always had better accuracy than Landsat 8.

Graphic and numeric reports of accuracy metrics do not hold information on the classification of individual land cover classes. Exhaustive details on LULC classes accuracies can be seen in Table 5, reporting the Producer’s accuracy (PA), a percent index that synthesizes the probability that a certain land cover on the ground has been correctly classified on the map produced, and the User’s accuracy (UA), the ratio between the correctly classified sites and the total number of classified sites.

PA and UA accuracies, with the exception of S-1, range between 0.61–1 and 0.48–1, respectively. The lowest Producer’s accuracy was observed for highway class (0.61), after L8 classification with NI(s), and the lowest User’s accuracy for forest (0.48), after S-1/S-2 classifications without NI(s).

Table 5 shows that the highway class was mapped with the highest PA (1) when the classification was based on S-1/L 8 and S-1/S-2/L 8 scenes without Ni(s). The highest PA and UA (1) were performed for the photovoltaic class using, respectively, S-1/S-2/L8 with NI(s) and S-2. High values of PA (0.98) was calculated for the plastic film class with Copernicus data combinations.

Although the use of normalized indices slightly increases the Overall Accuracy of classification, except for S-2 (Table 5; Figure 4), some deviations can be observed for the individual land cover types. The crops type, which was the predominant class, took advantage from using the complete set of optical and radar data when producer and user accuracies, for that class, were 87% and 90%, respectively. Built-up was better detected by the combination of radar and optical data (with NI(s)), both for PA and UA perspective. Forest and bare soil classes show the highest UA for classification performed combining all sensors dataset and normalized indices. Winter wheat PA best results were performed after the combination of Copernicus (S-1/S-2) datasets and NI(s). Photovoltaic class reaches the highest UA scores (1) for S-2 dataset, both with or without adding NI(s). The highest PA score of bare soil class (0.89) was achieved for S-1/S-2/L8 combination with NI(s), while the highest UA (0.95) was performed several times: combining both optical data (S-2/L 8) or radar and optical data (S-1/S-2/L 8) with NI(s), combining Copernicus data without NI(s).

The statistical significance of the difference in OA between the classifications performed with NI(s) was assessed using the McNemar test, with p < 0.05 (Table 5).

Matrix of McNemar test showing the statistical significance of differences between every classification pairs. The McNemar test revealed that the difference between the classification using the S-1/S-2/L8 and S-1 datasets were always statistically significant (p-value <0.05) in every comparison. Lesser values of significance were assessed for S-2/L8 when compared with S-1/S-2. Moreover, the higher value of significance was detected between S-2 and S-1/L8 classification (p-value = 0.922).

For the most accurate classifications (S-1/S-2/L8 and S-1/S-2), the accuracy of each land cover class was compared, both considering NI(s) (Figure 5). Differences in the producer’s accuracies between the quoted datasets were less than 6.25 for all classes, except for the photovoltaic and highway class, whose PA differences were 18% and −17.95%, respectively. Differences in the User’s accuracies were less than 8.35% for all classes, except for the photovoltaic and forest classes, whose UA differences were −12.33% and 20%, respectively.

In both land cover classifications, crops had the same UA value, so the percentage difference accuracy was equal to zero and does not appear in the UA chart (Figure 5). Built-up, bare soils, and crops, were more accurately classified with the complete combination of radar and optical datasets. Plastic films and highway classes were better classified with the Copernicus dataset combination.

Figure 6 shows the LULC classification maps on three example sites (A, B and C), performed with Sentinel-2 multispectral data (row 1), Sentinel-2 data combined with L 8 (row 2), with S-1 (row 3), and with S-1/L 8 (row 4), all with NI(s). The classified maps are graphically very similar for each inset. Border pixels between adjoining class types often show different boundaries drawn depending on the different spatial resolutions of source images and spectral contamination of pixels near field edges.

A heterogeneous area with agricultural crops, bare soils, plastic film covers, photovoltaic fields, and a highway is displayed in (Figure 6A). All the classified maps detect the photovoltaic field well, highlighted by the white circle, as well as the plastic films, represented in light blue. The bare soil areas are correctly mapped by every satellite combination. With S-2 classification (Figure 6 A1), a bridge over the highway (with a white arrow) and an area inside bare soil land (white square), were wrongly classified as built-up classes. The multispectral combination between S-2/L 8 (Figure 6 A2) and S-2/S-1/L 8 (Figure 6 A4) showed an increase of visual accuracy.

An urban area with its surrounding agricultural fields and bare soil is displayed in (Figure 6B). The confusion of built-up class (Figure 6 B1) was progressively reduced by the addition of a single Landsat 8 (Figure 6 B2) or Sentinel-1 data (Figure 6 B3) to Sentinel-2. Further reduction was performed by the addition of radar (Figure 6 B4) to optical combination. This reduced confusion can be seen directly in the confusion matrices for the whole municipal territory mapped (Table 6). High-density plantations (white arrows) were classified as forests instead of crops in S-2 and S-2/S-1 maps (Figure 6 B1 and B3). S-2/L 8 (Figure 6 B2) and S-2/S-1/L 8 (Figure 6 B4) were instead more accurate.

A less heterogeneous agricultural landscape, with crops, winter wheat, bare soil, and forest covers, is represented in Figure 6C. Land cover stands were well detected from every satellite single sensor or multispectral combination. A misclassifying of built-up instead of photovoltaic class was observed (white circle) in the Sentinel-2 derived map (Figure 6 C1). Here, each multispectral combination considered showed an improvement of accuracy respect to the use of S-2 (Figure 6 C2–C4).

In all the maps produced, an incorrect classification of built-up as bare soil or crops can be attributed to the diffuse presence, in several sites of the investigated territory, of outcropping rock.

Figure 7 displays the LULC classification of the whole territory of Acquaviva delle Fonti performed using S-1/S-2/L 8 band combination, with normalized indices. The detailed insets (Figure 7) and the whole territory classifications (Figure 7) demonstrate a surprisingly smooth result of the classified LULC, given that no prior segmentation technique was applied.

The large amount of multi-source input predictors data required a detailed assessment of single band importance for classification efforts. Figure 8 displays the Gini importance (Breiman et al., 1984), as standardized values, of the input features used for LULC classification. This variable gives quantitative information about the contribution of different bands to the Random Forest classifier used. A specific predictor plays a significant role in the classification when associated with a high Gini importance; instead, a low importance means that a specific feature determines only limited improvements.

For Sentinel-2 image (S-2), short wave infrared 2 band (B12) was the highest performing band, fully in line with other studies (Schuster et al., 2012; Immitzer et al., 2016), while red edge 3 band did not perform as well.

The near infrared band (B5) was the variable with less importance in Landsat 8 classification, both as a single sensor and combined with Sentinel-2 data (S-1/L 8 and S-2/L8). For the S-1/L 8 bands combination, vertical visualization (VV) of C radar data was the band with the higher contribution to detect LULC types of investigated territory. The Sentinel-2 green band (B3) did not perform as well for all optical and radar data combinations (S-1/S-2/L 8).

The importance of the normalized difference indices was quite different in the classifications performed. Green Normalized Difference Vegetation Index (GNDVI) was the higher performing band in optical combination (S-2/L8), even in combination with radar data (S-1/S-2/L8).

Focusing on visible bands, it appears (Figure 8) that those bands were unimportant in the Landsat OLI single sensor classification (L8) and in all the classifications performed in combination with that sensor data (S-1/L 8, S-2/L 8, and S-1/S-2/L 8), where the visible bands of the Sentinel-2 sensor also did not perform as well.

For Sentinel-1 image, VV band performed a variable importance slightly higher than VH band or VH and VV bands combination.