Sentinel-2 Data and Unmanned Aerial System Products to Support Crop and Bare Soil Monitoring: Methodology Based on a Statistical Comparison between Remote Sensing Data with Identical Spectral Bands

Abstract

The growing need for sustainable management approaches of crops and bare soils requires measurements at a multiple scale (space and time) field system level, which have become increasingly accurate. In this context, proximal and satellite remote sensing data cooperation seems good practice for the present and future. The primary purpose of this work is the development of a sound protocol based on a statistical comparison between Copernicus Sentinel-2 MIS satellite data and a multispectral sensor mounted on an Unmanned Aerial Vehicle (UAV), featuring spectral deployment identical to Sentinel-2. The experimental dataset, based on simultaneously acquired proximal and Sentinel-2 data, concerns an agricultural field in Pisa (Tuscany), cultivated with corn. To understand how the two systems, comparable but quite different in terms of spatial resolution and atmosphere impacts, can effectively cooperate to create a value-added product, statistical tests were applied on bands and the derived Vegetation and Soil index. Overall, as expected, due to the mentioned impacts, the outcomes show a heterogeneous behavior with a difference between the coincident bands as well for the derived indices, modulated in the same manner by the phenological status (e.g., during the canopy developments) or by vegetation absence. Instead, similar behavior between two sensors occurred during the maturity phase of crop plants.

1. Introduction

For several years, agriculture has been called upon to play its role as a supplier of food and other fundamental raw materials in a context of environmental sustainability and climate change [1,2,3]. This emergency was announced in the reports of the biggest international organizations, such as the Food and Agriculture Organization (FAO) and the United Nations (UN) [4,5,6]. As highlighted by Schiavon et al. (2021) [7], user needs of national and international institutional communities (i.e., policy makers, environmental agencies, regional government and local authorities) as well as private entities (i.e., insurance companies, businesses) require value-added information related to cropland mapping and phenological study, to cope with Common Agricultural Policy (CAP) 2023–2027 objectives, as sustainable agriculture production and natural resources management. For more than 40 years, remote sensing (RS) has contributed to the evolution of agricultural practice to support agricultural policy and economics. Exploiting various types of sensor in optical, thermal, and microwave spectral domains that are mounted on aircraft, satellites, drones or tractors, RS provides repetitive information on crop and soil status throughout the several seasonal steps, at different scales and for different actors [8,9]. Furthermore, the European Commission (EC), beginning in 2018 [10], has encouraged the use of new technologies such as Copernicus Sentinel data, integrated with EGNOS/Galileo, UAS data, geotagged photographs of EO products, possibly taking advantage of Copernicus In situ Component, for the control and granting of CAP payments by local authorities. This promotes open data with a common data-sharing approach. Nowadays, the technological innovation required in agricultural management, so-called smart farming [11], concerns the use of Information and Communication Technology (ICT) [12], big-data [13], Artificial Intelligence (AI) tools, specific sensors and integrated procedures of satellites and UAS data [14], to monitor crop growth and development, plant health, productivity and to manage nutrient and water use optimization programs. This is a great opportunity to boost the free and open market of Copernicus, the EO program of the European Union, as demonstrated by the analysis of the users’ technical and operational requirements for the expansion of the European space component and its services in different application domains [15].

Regarding the use of satellite and UAS data, many applications have been tested in precision agriculture (PA) [16,17,18]. The first satellite products emerged in the 1970s [19], although the spatial resolution was very low, such as Landsat data [20,21]. However, more recent advances in high-resolution satellites with a larger number of wavebands, as well as the advent of new satellite constellations such as Geo-Eye [22] or WorldView [23], have greatly facilitated the application of the technological achievements of remote sensing studies [8,9]. Nevertheless, unfavorable weather conditions, the timeliness of observation, high prices or not having suitable spatial resolution, still represent the most critical aspects in their practical applications for smart farming [9]. On the other hand, the UAS, also defined as proximal sensing, are applied inside the PA [24,25], through the redefinition of several techniques, from a soil plant analysis development (SPAD) method [26] to fluorescence sensing for nitrogen deficiencies [27,28]. UASs offer a series of advantages such as better spatial resolution, the ease of reaching narrow and difficult places accessible to humans and the ability to program the acquisition date [29]. On the contrary, the limitations of UAS are its not very extended tiles, the need for ongoing maintenance of its mechanical parts, problems with the battery or fuel, the short flying time and the small payload.

Considering that the two systems are potentially complementary, several studies have taken into account the hypothesis of integrating satellite and UAS data within the same time series [30]. The comparison of the results of certain elaborations carried out on the two remote sensing products suggest the possibility to integrate both data sources [31,32,33,34] or they specify particular use conditions [35,36,37,38,39]. In some cases, however, the spatial and spectral characteristics of the two considered systems are very different and/or the comparison method applied, such as linear regression only, might not consider all the variables involved.

The study presented focuses on a statistical comparison between data acquired by two multispectral instruments characterized by the same spectral bands, i.e., Sentinel-2 satellite sensor (S2) and MAIA S-2 (MS2) camera (SAL Engineering S.r.l.), mounted on an unmanned aerial vehicle. The purpose was to evaluate the feasibility of robust synergistic use of the two sensors, verifying the correspondence between the S2 and MS2_10 products. The ability to generate more information through sensor cooperation was analyzed and adapted to the crops and bare soil recognition. As already experimented by other authors in different contexts [40,41], the favorable and unfavorable conditions of UAS and satellite data synergistic use were studied. Therefore, the procedure provided for a statistical quantification of eventual differences found and to evaluate under which conditions they can be considered acceptable. The derived advantage is a greater awareness in the application of mixed time series to monitor and to discriminate crops and bare soil. As well as S2 [42,43,44,45], MS2 also found several applications due to its centimetric spatial resolution as PA, crop classification and assessment [46,47,48]. Using a corn crop field as a case study, the accurate statistical comparison between S2 and MS2 imagery was based on the Normalized Difference Vegetation Index (NDVI), [49] and the Soil Adjusted Vegetation Index (SAVI), [50]. Several studies are based on these vegetation indexes for different purposes, such as estimating living biomass in semi-arid areas [51] in pasture fields [52], identifying tree infections [53], assessing crop hail damage [54], comparing two indexes’ reliability with different land use [55] and to cross-compare between two different sensor products [56,57]. Considering different crop vegetal cover density and their time evolutions, the indices were compared using the Welch‘s test [58]. This test can perform when sample variances are not equal and allows easy inferences to be made using a null hypothesis based on the mean. Then, to better understand the results, the same test was also applied on corresponded Near Infrared (NIR) and Red bands.

In the end, as suggested by other studies [31,37,59,60], the linear regression between S2 and MS2 index values was calculated to verify the reliability of its use in the various scenarios analyzed.

2. Study Area and Materials

2.1. Study Area

The study area is located in Madonna dell’Acqua in the municipality of San Giuliano Terme (Pisa), a few kilometers north of Pisa (Figure 1). The field selected for the joint survey (S2-MS2) was cultivated in 2019 with corn (Zea mays) in an area of about 10 ha.

The soils that characterize the area, and in particular the studied field, are mainly Eutric Cambisols, according to the World Reference Base (WRB) taxonomy [61]. These are relatively young soils, featuring a more or less evolved subsurface horizon that presents pedogenetic alterations of a chemical/physical nature. These occurred on site (cambic horizon), and are not due to illuviation processes. The soils are very deep, with an Ap-Bw profile (horiz. Cambico) -Bg. They are non-gravelly, with a sandy to silty texture, are non-calcareous to slightly calcareous, with a neutral to moderately alkaline reaction, with saturation in bases very high and are moderately well drained. The field has conventional transversal drainage channels for the entire extension, with an average depth compared to the ground level of about 50 cm.

2.2. Data

Sentinel-2 comprises two polar-orbiting satellites, in the same sun-synchronous orbit: Sentinel-2A, launched in 2015, and Sentinel-2B, launched in 2017. They carry an optical instrument, called MultiSpectral Instrument (MSI), that collects energy in 13 different spectral bands with several spatial resolutions (

Table 1. Characteristics of the two sensors used, in terms of band name, central wavelength and spatial resolution. The corresponding central wavelengths between the two instruments are highlighted in bold.

Sentinel-2			MAIA S2
Band Name	Central Wavelength (nm)	Spatial Resolution (m)	Band Name	Central Wavelength (nm)	GSD (m)
B1—Coastal aerosol	443	60	S1—Violet	443	0.047
B2—Blue	490	10	S2—Blue	490	0.047
B3—Green	560	10	S3—Green	560	0.047
B4—Red	665	10	S4—Red	665	0.047
B5—Vegetation Red Edge	705	20	S5—Red Edge 1	705	0.047
B6—Vegetation Red Edge	740	20	S6—Red Edge 2	740	0.047
B7—Vegetation Red Edge	783	20	S7—NIR 1	783	0.047
B8—Narrow NIR	842	10	S8—NIR 2	842	0.047
B8A—NIR	865	20	S9—NIR 3	865	0.047
B9—Water vapour	945	60	/	/	/
B10—SWIR - Cirrus	1.375	60	/	/	/
B11—SWIR	1.610	20	/	/	/
B12—SWIR	2.190	20	/	/	/

). The data used were downloaded from the Theia web page (https://www.theia-land.fr/en/satellite-data/ (accessed on 30 September 2020). Theia is a Scientific Expertise Center conducting research on remote sensing, such as mobilization satellites and airborne and in situ data on ‘continental surface’ issues.

The S2 data used are level 2-A. These data were corrected for atmospheric effects using MAJA software [62]. This was developed by the Centre d’Etudes Spatiales de la Biosphère (CESBIO), which developed the methods and a prototype, and by the Centre National d’Etudes Spatiales (CNES) which funded the operational version of the processor. In the case of the S2 satellite, the first step consists in estimating and subsequently correcting for gaseous absorption. The next step is the creation of a composite image with the unclouded pixels from the processed data used as a reference for cloud detection [63]. After a multitemporal test, which is the sign of cloud presence, a final test measures the correlation of the pixel neighborhood in relation to the previous images. If a large correlation is observed, the pixel is, ultimately, not declared as a cloud. After the aerosol optical thickness (AOT) is estimated [64], it is possible to retrieve the surface reflectance, using look-up tables (LUT). Before editing the output image, the software corrects the adjacency effects and the effects of terrain slopes on illumination [65]. The S2 radiometric resolution is 12 bit, meaning the sensor can record 4096 energy levels. This range is converted into floating point type (from 0 to 1). S2 pixel values make ‘original’ reflectance × 10,000 (factor scale), but, at times, this value may be higher than 10,000 due to specific angular reflectivity effects [66].

MAIA S-2 is a multispectral camera that permits the simultaneous acquisition of high-resolution images at various wavelength intervals in Visible (VIS) and NIR electromagnetic spectrum regions. In particular, it is equipped with the same first nine wavelength intervals of the S2 satellite (Table 1). The MS2 system comprises of an array of nine monochromatic sensors that collect energy in unison thanks to global shutter technology, permitting it to be acquired ‘one-off’ for synchronized multiband measurements [67]. It was mounted on an UAV: a DJI S-900. For each performed survey, a single multispectral orthophoto was generated using the 300 acquired images. The orthorectification process was carried out using photogrammetric processing with Metashape-Agisoft software. The software uses Structure from Motion algorithms for the homologous point definition and for the Bundle Adjustment application based on the Ground Control Point (GCP) use, [68]. All flights were programmed and performed in the same way: an altitude of 100 m Above Ground Level (AGL) and strips able to ensure longitudinal overlap of over 80% and a lateral overlap of about 35% between the frames. For each frame, the Ground Sample Distance (GSD) equal to 4.7 cm and the footprint of 63 m × 48 m were due to the combination of study area size (9.1 ha), its almost rectangular shape and multispectral camera sensors of dimensions 4.8 mm × 3.6 mm with pixel size of 3.75 microns. For correct Bundle Adjustment process application, 8 GCPs were placed on the ground (Figure 1), uniformly distributed over the entire study area. All the GCPs, clearly visible on the images, were measured with dual frequency GNSS (Global Navigation Satellite System) geodetic instrumentation (Topcon GB500) equipped dual-frequency geodetic antennas (Topcon PG-A1 L1/L2). To guarantee centimeter accuracy (1–2 cm in the three dimensions), the single base static GNSS method was used. The reference instrument (Master Station) was positioned near the study area (maximum distance from the GCPs of less than 1 km).

To compare MS2 and S2 data, the S2 radiometric resolution was reduced from 12 bit to 8 bit. To achieve this result, the value of each DN was divided by 10,000 and then multiplied by 255 using a Raster Calculator in QGIS. To guarantee an accurate comparison between the two instruments, three acquisition dates were chosen based on the corn phenological cycle. For June and August there was a perfect match between MS2 and S2 with about half an hour’s difference between the two flights. Unfortunately, for July, the contemporaneity of the two acquisitions was hampered as the satellite data was affected by the presence of clouds. The available date closest to the proximity acquisition of the satellite data without cloud cover was therefore chosen (

Table 2. Sentinel-2 and MAIA S-2 acquisition time.

	Acquisition Date	Time (UTC+1)	Sun Azimuth (°)	Sun Elevation (°)
Sentinel-2	3 June 2019	10:18:45	146.13	24.62
	16 July 2019	10:28:42	148.1	25.26
	5 August 2019	10:28:41	151.76	29.32
MAIA-S-2	3 June 2019	12:00:00	193.48	55.96
	11 July 2019	12:00:00	190.43	56.1
	5 August 2019	12:00:00	189.29	51.07

).

Before the statistical test application, to also make S2 and MS2 data geometrically comparable, the MS2 spatial resolution was resampled from 0.047 m to 10 m, using a grid framework that derived from S2. The MS2 resampled data will be later called MS2_10.

3. Methods

3.1. Georeferencing

The S2 orthoimages are in UTM-WGS84 projection [69]. Concerning MS2, to ensure final centimeter accuracy (1–2 cm), the Master Station point geodetic coordinates were determined by a long-term static GNSS survey based on the connection to the three EUREF network permanent stations: GENO, IGMI, ELBA. The chosen geodetic reference system was IGb14—epoch 2019.5, projected in the UTM-WGS84 cartogaphic system. For the processing of GNSS data of all three surveys, the precise ephemeris distributed by IGS (International GNSS Service) was always used. The final orthoimage accuracy, confermed by some GCPs not used for the Bundle Adjustment calculation, was always less than 5 cm, i.e., less than the size of the GSD. This result, also in relation to the GCP number choice, is in agreement with that reported by [70].

A quality control of the relative georeferencing between the MS2_10 and S2 images was required. In fact, if the georeferencing was not correct, the value measured in the pixel $\left(\mathrm{i},\mathrm{j}\right)$ of the first data would not correspond to that contained in the pixel $\left(\mathrm{i},\mathrm{j}\right)$ of the second data as they would cover different portions of land. Therefore, to achieve the correct overlap, it would be necessary to shift one pixel in relation to the other by a certain lag $\left(\mathrm{u},\mathrm{v}\right)$. The square of the radiometric intensity difference (Equation (1)) is a measure of the fit between an image f $\left(\mathrm{x},\mathrm{y}\right)$ and a feature t $\left(\mathrm{x},\mathrm{y}\right)$ shifted by $\left(\mathrm{u},\mathrm{v}\right)$:

$${\mathrm{d}}_{\mathrm{f},\mathrm{t}}^2\left(\mathrm{u},\mathrm{v}\right)={\displaystyle \sum }\mathrm{x},\mathrm{y}{\left[\mathrm{f}\left(\mathrm{x},\mathrm{y}\right)-\mathrm{t}\left(\mathrm{x}-\mathrm{u},\mathrm{y}-\mathrm{v}\right)\right]}^2$$

(1)

where the ∑ is extended on the window containing the feature t translated by $\left(\mathrm{u},\mathrm{v}\right)$. To identify the minimum of the distance ${\mathrm{d}}_{\mathrm{f},\mathrm{t}}^2\left(\mathrm{u},\mathrm{v}\right)$, Lewis (1995) [71] proposed a robust algorithm that search the maximum of the Normalized Cross-Correlation matrix γ$\left(\mathrm{u},\mathrm{v}\right)$, (Equation (2)):

$$\mathsf{\gamma }\left(\mathrm{u},\mathrm{v}\right)=\frac{{{\displaystyle \sum }}_{\mathrm{x},\mathrm{y}}\left[\mathrm{f}\left(\mathrm{x},\mathrm{y}\right)-{\overline{\mathrm{f} }}_{\mathrm{u},\mathrm{v}}\right]-\left[\mathrm{t}\left(\mathrm{x}-\mathrm{u},\mathrm{y}-\mathrm{v}\right)- \overline{\mathrm{t}}\right]}{\sqrt{{{\displaystyle \sum }}_{\mathrm{x},\mathrm{y}}{\left[\mathrm{f}\left(\mathrm{x},\mathrm{y}\right)-{\overline{\mathrm{f} }}_{\mathrm{u},\mathrm{v}}\right]}^2{{\displaystyle \sum }}_{\mathrm{x},\mathrm{y}}{\left[\mathrm{t}\left(\mathrm{x}-\mathrm{u},\mathrm{y}-\mathrm{v}\right)- \overline{\mathrm{t}}\right]}^2}}$$

(2)

where:

f $\left(\mathrm{x},\mathrm{y}\right)$ is an image;

t $\left(\mathrm{x},\mathrm{y}\right)$ is a feature;

$\overline{\mathrm{t}}$ is the mean of the feature;

${\overline{\mathrm{f} }}_{\mathrm{u},\mathrm{v}}$ is the mean of f $\left(\mathrm{x},\mathrm{y}\right)$ in the region under the feature.

The sum is over $\mathrm{x}, \mathrm{y}$ under the window containing the feature t positioned at $\left(\mathrm{u},\mathrm{v}\right)$. The maximum of γ$\left(\mathrm{u},\mathrm{v}\right)$ corresponds to the feature’s most right index [71].

The calculation was performed by comparing the MS2_10 and S2 bands over the study area in the three different phenological stage [72]. The maximum of the cross-correlation matrix normalized was obtained for indexes (8,8) i.e., for a lag (0, 0), considering that the image size is an 8 × 8 pixel. This means that no shifts improve the correlation achieved in georeferencing.

Figure 2 shows an example of the values obtained for the correlation coefficients matrix, for the observable NIR (June date).

3.2. Data Accuracy

For each band, the slightest noticeable radiometric difference between S2 and MS2_10 data depends on the acquisition accuracy level of the two sensors.

The MS2_10 translates the received radiation intensity (photons) into a discrete DN representable on 256 radiometric levels. This range was converted into floating point type, from 0 to 1. Numerous experimental tests, implemented by the authors with technicians from SAL Engineering of Modena (Italy), showed that the MS2 (8 bit) measurement have a repeatability of one bit. The tests were carried out on few hundred images defining the black level of all the CMOS sensors that make up the MAIA-S2 camera. The optical black pixel method [73] was used to define the dark voltage value. The test results have confirmed the stability of DN in correspondence of all the pixels of images, with a black level constant value of DN = 9. Therefore, we can assume an accuracy of τ = 1/256~0.004. This value is considered the resolution of the measurement, i.e., the minimum discriminable quantity.

Concerning S2, the original signal discrimination occurs at 12 bits but a resampling at 8 bits was carried out. This range was converted into floating point type, from 0 to 1. Similar considerations to the previous ones can be made about S2 radiometric accuracy.

Based on the assumed hypothesis, the differences between S2 and MS2 measures were significant if they exceeded τ’= 0.005.

3.3. Test Area Choice and Vegetation Index Applications

Significant within-field spatial variability exists in factors that influence crop yield; this variability can be measured [74]. Consequently, it is important that the resulting information is as correct as possible in order to be used to improve crop management practices, in relation to increasing profit and decreasing environmental impact. In order to analyze several vegetation spatial distribution scenario, three test areas with a different percentage of vegetation cover were identified inside the corn crop field. The areas had to be characterized by a high, medium and low vegetation percentage, respectively. For this purpose, a procedure based on the chlorophyll content map was applied. The chlorophyll content was used as a parameter of maturity and quality of food crops [75,76]. Several authors have showed that reflectance in the green and/or red-edge spectral regions is an optimal parameter for non-destructive estimation of leaf chlorophyll content in a wide range of its variation [77,78,79,80]. They used different vegetation indexes, such as the MERIS Terrestrial Chlorophyll Index (MTCI), [81], Enhanced Vegetation Index (EVI2), [82], and Red-edge Chlorophyll Index (CI_Red-Edge), [80,83]. In this paper, the considered vegetation index was the green Chlorophyll Index (CI_green, Equation (3)), [75]:

$${\mathrm{CI}}_{\mathrm{green}}=\frac{{\mathsf{\rho }}_{\mathrm{NIR}}}{{\mathsf{\rho }}_{\mathrm{Green}}}-1$$

(3)

where ${\mathsf{\rho }}_{\mathrm{NIR}}$ and ${\mathsf{\rho }}_{\mathrm{Green}}$ are NIR and green reflectance, respectively.

As showed by [84,85,86], this index has a strong linear correlation with chlorophyll content at the leaf scale.

The procedure applied a moving window over the CI_Green map and calculated the statistical parameters on the window elements. In this way, the characteristics of homogeneity and intensity were calculated in correspondence with the neighborhood of the current position (i, j). The calculated statistical parameters are the average for the intensity and standard deviation for dispersion. The association of the vegetal percentage was based on the histogram of the average and standard deviation values referred at the floating window centers.

Several suitable positions have been identified. These showed the desired characteristics, i.e., high value and low dispersion of CI_Green (area A), low intensity and high dispersion of CI_Green (area C), and intermediate conditions for area B. Then, an expert operator chose the three definitive zones from the suitable positions. He defined the test area boundaries based on expeditious field surveys and a photo-interpretation of the MS2 data. Square areas measuring 80 × 80 m (equivalent to an 8 × 8 pixel) were chosen to facilitate subsequent comparisons.

Then, NDVI [49] and SAVI [50] indices (Equations (4) and (5)) were applied both on S2 and MAIA_S2 data. For both indices, the bands used were B4 (Red) and B8 (NIR) bands for Sentinel-2 and S4 (Red) and S8 (NIR) bands for MAIA S2_10.

$$\mathrm{NDVI}=\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}}$$

(4)

Concerning SAVI index, the Huete equation was used:

$$\mathrm{SAVI}=\left[\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}+\mathrm{L}}\right]\left(1+\mathrm{L}\right)$$

(5)

where:

L = 0.75 for June data and A, B, C areas;

L = 0.50 for July data and B, C areas;

L = 0.25 for July data and A area.

Considering that the proximity data were resampled to the S2 data resolution, the choice of L values was based on their predominant use in the satellite data processing [50,87]. Considering the total vegetation cover showed in August data, the SAVI was not calculated for this epoch.

3.4. Statistical Test

When the index differences exceeded τ’, it was necessary to evaluate their statistical significance relating their value to the associated uncertainty level.

After the three test areas’ identification with different homogeneous vegetation cover, the obtained radiometric parameters from the S2 and MS2_10 vegetation index values were compared. The comparison is validated by statistical tests.

The radiometric homogeneity of each test area allowed its 64 pixels (n) belonging to the same statistical distribution to be considered. The greater the test area homogeneity, the greater the robustness of this hypothesis. Therefore, both for the 64 S2 pixels and for MS2_10, sample mean and variance $\left({ \overline{\mathrm{x}}}_{\mathrm{k}},{\mathrm{s}}_{\mathrm{k}}^2\right)$ were calculated (Equation (6)):

$${ \overline{\mathrm{x}}}_{\mathrm{k}}={\displaystyle \sum }_{\mathrm{i}}^{\mathrm{n}}\frac{{\mathrm{x}}_{{\mathrm{k}}_{\mathrm{i}}}}{\mathrm{n}};{ \mathrm{s}}_{\mathrm{k}}^2={\displaystyle \sum }_{\mathrm{i}}^{\mathrm{n}}\frac{\left({\mathrm{x}}_{{\mathrm{k}}_{\mathrm{i}}}-{ \overline{\mathrm{x}}}_{\mathrm{k}}\right)}{\left(\mathrm{n}-1\right)}$$

(6)

And ${\mathrm{s}}_{\overline{{\mathrm{x}}_{\mathrm{k}}}}$ was the standard deviation associated with the sample mean (Equation (7)):

$${\mathrm{s}}_{\overline{{\mathrm{s}}_{\mathrm{k}}}}=\frac{{\mathrm{s}}_{\mathrm{k}}}{\sqrt{\mathrm{n}}}$$

(7)

Between the two means with variance estimated from the samples they belong to, the equality test was based on the null hypothesis of the expected values equality (Equation (8)). Namely, the expected value of the difference between corresponding values was considered equal to zero. In this context, the corresponding values were related to pixels of equal position (Equation (8)):

$${\mathrm{H}}_0:\mathrm{E}\left\{{ \overline{\mathrm{x}}}_1\right\}=\mathrm{E}\left\{{ \overline{\mathrm{x}}}_2\right\},{ \mathrm{H}}_0:\mathrm{E}\left\{{ \overline{\mathrm{x}}}_1-{ \overline{\mathrm{x}}}_2\right\}=0$$

(8)

where:

$\mathrm{k}=1$ for S2

$\mathrm{k}=2$ for MS2_10.

The hypothesis was tested through Welch’s t-test. The null hypothesis (Equation (8)) based on the mean offered an intuitive interpretation. The nonparametric tests based on the comparison of ranks, and which generally have less power, have not been chosen. Tests based on trimmed means have the disadvantage that the amount of trimming must be arbitrarily chosen. On the other hand, the Welch test requires the normality of the distribution of the samples, which does not always occur [88].

The considered statistic was the ratio (Equation (9)) between the difference of the index means and its standard deviation estimated by the two samples (Equation (10)):

$${\mathrm{t}}_{\mathsf{\upsilon }}=\frac{\mathrm{d} \overline{\mathrm{x}}}{{\mathrm{s}}_{\mathrm{d} \overline{\mathrm{x}}}}$$

(9)

$$\mathrm{d} \overline{\mathrm{x}}={ \overline{\mathrm{x}}}_1-{ \overline{\mathrm{x}}}_2;{ \mathrm{s}}_{\mathrm{d} \overline{\mathrm{x}}}=\sqrt{\frac{{\mathrm{s}}_1^2}{{\mathrm{n}}_1}+\frac{{\mathrm{s}}_2^2}{{\mathrm{n}}_2}}$$

(10)

where, with i = 1,2:

${\overline{\mathrm{x}}}_{\mathrm{i}}$ = sample mean;

${\mathrm{s}}_{\mathrm{i}}$ = sample standard deviation;

${\mathrm{n}}_{\mathrm{i}}$ = sample size.

In the Welch’s t-test, the degrees of freedom associated is approximated given by Equation (11):

$$\mathsf{\upsilon }\approx \frac{{\left(\frac{{\mathrm{s}}_1^2}{{\mathrm{n}}_1}+\frac{{\mathrm{s}}_2^2}{{\mathrm{n}}_2}\right)}^2}{\frac{{\mathrm{s}}_1^4}{\left({\mathrm{n}}_1^2{ \mathsf{\upsilon }}_1\right)}+\frac{{\mathrm{s}}_2^4}{\left({\mathrm{n}}_2^2{ \mathsf{\upsilon }}_2\right)}};{ \mathsf{\upsilon }}_{\mathrm{i}}={\mathrm{n}}_{\mathrm{i}}-1$$

(11)

Considering that the two compared samples had the same element number ${\mathrm{n}}_1={\mathrm{n}}_2=\mathrm{n}$, the Equations (9) and (11) have become:

$${\mathrm{t}}_{\mathsf{\upsilon }}=\frac{\sqrt{\mathrm{n}}\left({ \overline{\mathrm{x}}}_1-{ \overline{\mathrm{x}}}_2\right)}{\sqrt{{\mathrm{s}}_1^2+{\mathrm{s}}_2^2}};$$

(12)

$$\mathsf{\upsilon }=\frac{\left(\mathrm{n}-1\right)\left({\mathrm{s}}_1^2+{\mathrm{s}}_2^2\right)}{\left({\mathrm{s}}_1^4+{\mathrm{s}}_2^4\right)}$$

(13)

The Welch test t variable follows to Student t distribution at ν degree of freedom (Equations (12) and (13)). The ${\mathrm{t}}_{\mathsf{\nu }}$ value was compared with the assumed critical value ${\mathrm{t}}_{\mathrm{Wc}}$ corresponding to a confidence chosen level.

It is underlined that the previous test concerned the difference between the average values and did not consider the oscillation of the differences around the average observed in the various pixels. Therefore, the differences between the averages reflected effects obviously affecting the entire area.

The used ${\mathrm{s}}_{\mathrm{d}\overline{\mathrm{x}}}$ can help to evaluate the fluctuation extent of the values in the various pixels around the average of the considered area.

3.5. Linear Regression

The evaluation of the correlation between the values obtained by the two sensors in the corresponding pixels was performed by linear regression (according to the least squares procedure). For each regression, $\mathrm{Y}=\mathrm{p}1 \mathrm{X}+\mathrm{p}0$ the following quantities were calculated: regression parameter $\mathrm{p}1$ and $\mathrm{p}0;$ their standard deviations ${\mathsf{\sigma }}_{\mathrm{p}1}, {\mathsf{\sigma }}_{\mathrm{p}0}$; the standard error of the regression ${\mathsf{\sigma }}_0$; correlation coefficients; the determination coefficient R² and linear correlation coefficients $\mathsf{\rho }$.

4. Results

4.1. Test Areas

To identify several potential test areas, a procedure was implemented. It was based on the statistical parameter calculation of the pixels contained in a moving window on the image. The moving window was applied on the MS2_10 CI_Green map that originated from the MS2 data, characterized by higher geometric resolution (GSD 4.7 cm). Only the original July map was used to recognize the three scenarios due to it including a nonhomogeneous spatial vegetation distribution overall. On the other hand, the June and August MS2 CI_Green maps were not suitable: the first showed mainly bare soil due to the majority of corn plants not yet to having sprouted in that month; the second map showed a homogenous plant distribution over the whole crop field and the differentiation among vegetation densities was not possible.

In order to quantify and to discriminate the different vegetation distributions, the considered criteria were: high average (>4.4) and low dispersion (<1.4) identify area A, low mean (<4.0) and low dispersion (<1.1) identify area C; mean values included from 4.05 and 4.15 and high dispersion values (>1.75) describe area B. These values were chosen by analyzing the histograms of mean and dispersion of the CI_Green in the moving windows.

Figure 3 shows the pixels that satisfied both previous conditions (in yellow) in the corn field. To simplify the moving window application, the field surface was included in a larger rectangular area (in purple).

The agronomic expert operator then defined the final size of the test areas around the identified pixels, i.e., 80 × 80 m (8 × 8 S2 pixels), (Figure 4). Concerning test area B, to clearly separate it from test areas A and C, the expert preferred to consider the northernmost zone among those available.

4.2. Vegetation Index and Statistical Test Results

For each analyzed stage and for each test area, the NDVI and SAVI averages were calculated on the 64 vegetation index values included in each test area. For each average value, the corresponding root-mean-square error (r.m.s.e.) values were reported (

Table 3. For each test area and for each stages: NDVI and SAVI average values and corresponded r.m.s.e. for both sensors, S2—MS2_10 difference of vegetation index mean values, Welch’s test output. In bold, the values show the cases in which the equality hypothesis was accepted. In italics, the values show the cases in which the difference were considered not significant.

			NDVI		SAVI
Epoch	Test Area		Sentinel-2	MAIA S-2	Sentinel-2	MAIA S-2
June	A	mean	0.247	0.155	0.160	0.103
		r.m.s.e.	0.043	0.047	0.030	0.032
		difference	0.093		0.057
		W test	11.6		10.4
	B	mean	0.269	0.177	0.160	0.110
		r.m.s.e.	0.047	0.053	0.028	0.033
		difference	0.092		0.051
		W test	10.4		9.3
	C	mean	0.205	0.117	0.124	0.073
		r.m.s.e.	0.038	0.041	0.021	0.024
		difference	0.090		0.103
		W test	12.7		14.1
July	A	mean	0.871	0.776	0.685	0.602
		r.m.s.e.	0.013	0.048	0.020	0.038
		difference	0.095		0.083
		W test	15.3		15.4
	B	mean	0.856	0.728	0.575	0.488
		r.m.s.e.	0.020	0.066	0.033	0.051
		difference	0.128		0.087
		W test	14.8		11.4
	C	mean	0.819	0.617	0.477	0.371
		r.m.s.e.	0.035	0.081	0.045	0.035
		difference	0.202		0.106
		W test	18.2		14.8
August	A	mean	0.947	0.946	Not applicable
		r.m.s.e.	0.005	0.010
		difference	0.001		Not applicable
		W test	0.9
	B	mean	0.956	0.948	Not applicable
		r.m.s.e.	0.007	0.009
		difference	0.007		Not applicable
		W test	4.9
	C	mean	0.955	0.950	Not applicable
		r.m.s.e.	0.013	0.014
		difference	0.006		Not applicable
		W test	2.4

).

The S2—MS2_10 vegetation index differences were calculated based on these average values (Table 3). When these differences exceeded the minimum appreciable radiometric resolution threshold τ’, the equality hypothesis of their expected values was tested. To confirm the reliability of the Welch’s test, it was applied on all test areas and all stages.

The difference values were compared with the assumed critical value ${\mathrm{t}}_{\mathrm{Wc}}$ and the hypothesis of equality of the means is accepted for t_v <${ \mathrm{t}}_{\mathrm{Wc}}$ = 2.6 (confidence level 99%), (Table 3).

To better understand the results obtained from the indices, the Welch’s test was also applied to the bands used, NIR and Red.

For each analyzed epoch and for each test area, the NIR and Red averages were calculated on the 64 vegetation index values included in the test area. For each average value, the corresponding root-mean-square error (r.m.s.e.) values were reported (

Table 4. For each test area and for each stages: NIR and Red average values and corresponded r.m.s.e. for both sensors, skewness values, S2—MS2_10 difference of band mean values, Welch’s test output. In bold, the values show the cases in which the equality hypothesis was accepted. In italics, the values show the cases in which the differences were considered not significant.

			NIR		Red
Epoch	Test Area		Sentinel-2	MAIA S-2	Sentinel-2	MAIA S-2
June	A	mean	0.274	0.268	0.160	0.192
		r.m.s.e.	0.024	0.026	0.028	0.036
		skewness	0.74	0.49	−0.03	−0.06
		difference	0.005		−0.032
		W test	1.1		5.5
	B	mean	0.246	0.244	0.142	0.171
		r.m.s.e.	0.018	0.022	0.013	0.019
		skewness	0.98	0.22	0.02	0.41
		difference	0.002		−0.029
		W test	0.5		−10.0
	C	mean	0.242	0.241	0.165	0.198
		r.m.s.e.	0.032	0.035	0.029	0.038
		skewness	0.80	0.63	0.65	0.59
		difference	0.001		−0.033
		W test	0.1		−5.5
July	A	mean	0.396	0.363	0.027	0.046
		r.m.s.e.	0.019	0.017	0.002	0.011
		skewness	−0.19	0.25	0.08	2.35
		difference	0.033		−0.019
		W test	10.2		−13.6
	B	mean	0.377	0.349	0.029	0.055
		r.m.s.e.	0.028	0.024	0.003	0.012
		skewness	−0.05	−0.16	−0.16	0.34
		difference	0.028		−0.026
		W test	6.0		−16.5
	C	mean	0.340	0.318	0.033	0.076
		r.m.s.e.	0.024	0.016	0.005	0.020
		skewness	−0.08	0.20	1.47	1.29
		difference	0.020		−0.042
		W test	6.1		−16.6
August	A	mean	0.436	0.436	0.012	0.012
		r.m.s.e.	0.015	0.022	0.002	0.001
		skewness	0.41	0.09	0.61	1.37
		difference	0.000		0.000
		W test	0.1		0.8
	B	mean	0.401	0.405	0.011	0.009
		r.m.s.e.	0.019	0.026	0.002	0.001
		skewness	−0.11	−0.02	−0.06	1.18
		difference	−0.004		0.002
		W test	−0.9		−6.2
	C	mean	0.419	0.424	0.011	0.010
		r.m.s.e.	0.014	0.025	0.003	0.002
		skewness	−1.01	−0.50	1.51	1.78
		difference	−0.005		−0.001
		W test	−1.5		3.4

). Skewnness values were also shown in the same table. The skewness values are generally low and comparable for both samples, except for some MS2_10 Red samples.

The S2—MS2_10 band differences were calculated based on these average values (Table 4). When these differences exceeded the minimum appreciable radiometric resolution threshold τ’, the equality hypothesis of their expected values was tested. To confirm the reliability of the Welch’s test, it was applied on all test areas and all epochs.

The difference values were compared with the assumed critical value ${\mathrm{t}}_{\mathrm{Wc}}$ and the hypothesis of equality of the means is accepted for ${\mathrm{t}}_{\mathsf{\nu }}$ < ${\mathrm{t}}_{\mathrm{Wc}}$ = 2.6 (confidence level 99%), (Table 4).

The results of Table 4 show two anomaolous cases where the equality hypothesis was not accepted, despite the corresponded differences being lower than the significance level τ’: B and C test area, for the Red band and August epoch. In these cases, the two data were assumed to be equivalent.

4.3. Linear Regression: MS2_10 NDVI vs. S2 NDVI

The regressions of MS2_10 NDVI values towards those of S2 NDVI present in the corresponding pixels were performed.

Linear regression estimates the linear law $\mathrm{Y}=\mathrm{p}1 \mathrm{X}+\mathrm{p}0$ which transforms $\mathrm{X}$ measures (S2) into $\mathrm{Y}$ measures (MS2_10). The functional relationship is all the more reliable the less the experimental points are dispersed around the regression, i.e., the lower the value ${\mathsf{\sigma }}_0$ and the parameter’s standard deviation ${\mathsf{\sigma }}_{\mathrm{p}1},{ \mathsf{\sigma }}_{\mathrm{p}0}$.

Figure 5 shows some interesting results compared to the Welch’s test finding, as commented in the Discussion section. The graphs report the NDVI value points (S2, MS2_10), the estimated regression line and the lines representing the expected range for a new observation at the confidence level of 0.99. For the same samples,

Table 5. Estimate of the parameters ($\mathrm{p}1, \mathrm{p}0$), their uncertainties (${\mathsf{\sigma }}_{\mathrm{p}1},{ \mathsf{\sigma }}_{\mathrm{p}0}$ ), the standard error of the regression (${\mathsf{\sigma }}_0$), the determination and linear correlation coefficients (R², $\mathsf{\rho }$) relating the five graphs of Figure 5.

Epoch	Test Area	p1	σp1	p0	σp0	σ0	R2	ρ
June	B	0.83	0.04	0.12	0.01	0.018	0.87	0.93
July	B	0.27	0.02	0.66	0.01	0.009	0.79	0.89
July	C	0.40	0.02	0.57	0.01	0.015	0.84	0.92
August	A	0.26	0.05	0.7	0.05	0.004	0.28	0.53
August	C	0.60	0.08	0.39	0.08	0.009	0.47	0.68

shows the estimate of the parameters ($\mathrm{p}1, \mathrm{p}0)$, their uncertainties (${\mathsf{\sigma }}_{\mathrm{p}1},{ \mathsf{\sigma }}_{\mathrm{p}0})$, the standard error of the regression (${\mathsf{\sigma }}_0$) and the determination and linear correlation coefficients (R², $\mathsf{\rho }$).

5. Discussion

Beginning from a short temporal series of proximal measurements derived from a UAV platform sensor compliant to the S2 (MAIA-S2, par. 2.2), a comparison analysis was developed between these data and the mostly contemporary S2 dataset. The study aimed to highlight their temporal behaviors through a robust statistical analysis between the two data and their derived VIs in a corn field (Figure 1 and par. 2.1), from the flowering to the ripening phase (June-August 2019). The main purpose was to verify the relationship between the dynamics of the canopy and the bare soil component, useful for precision agriculture approaches and to optimize the cropping systems in a more sustainable way.

Considering the most resolute S2 pixel size (e.g., 10 × 10 m), and the field features (e.g., phenological phase, microtopography, moisture and agriculture practices) it is realistic to expect that an ever-different mixture of canopy, bare soil and other elements affected the radiance signal provided by individual pixels during the research time span. In this perspective, the resampling of the proximity MS2 data (almost 5 cm) in the S2 satellite data resolution (10 m) was carried out. Precisely, the desire to study a proximity data role in a regional-scale analysis, compared to the mid-resolution, prompted the use of this uncommon approach [37,89]. This choice was made as many crop field activities are carried out on a regional scale, and will be even more so in the future, through free spatial mid-resolution satellite data and with temporal resolution of a few days (e.g., S2 and Landsat). Instead, the proximity data, even if more suitable for local scale studies, are more expensive to acquire and to process, both in time and economic terms.

After the MS2 data resampling (generation of the MAIA-S210 product), the quality of the relative georeferencing of the two data types (par. 2.2.1) was evaluated to correctly proceed with the next experimental design operations. The spatial optimal combination, in this case probably derived from the Madonna dell’Acqua site conditions, i.e., horizontal surface without rugged topography, require that the UAV data acquisition be carried out with the highest reference standards.

The strategy to select three A, B, C test sub-areas (par. 3.1) allowed the adequate exploitation of the homogeneity of pixel radiometric characteristics within the areas and to represent, by their sum, the heterogeneity of the patterns of crop and bare soil in the field under observation. This condition corresponded to the parametric test base hypothesis, i.e., that the measurements were realizations of a single variable. Therefore, the measurements constituted a random sample of realizations of the variable. Mean and variance of the variable were unknown, such as samples made with MAIA S2 and Sentinel-S2. There are three main advantages deriving from the analysis of subareas with different vegetation/bare soil cover:

−

it allows the identification of different behaviors of the relationship between bare soil and crop, and their mixture;

−

it allows the application of the most suitable and diversified analysis strategy within each subarea until the completion of the entire plot. The mobile window used was also essential for characterizing the entire plot and,

−

the transfer of the method to different types of crop systems characterized by areas with different percentages of vegetation cover.

In effect, for a specific multispectral analysis, it is necessary to take into account the realistic possibility of having, within cultivated plots with extensive crops, areas of even limited extensions with a different surface coverage ratio between canopy and soil. This ratio can be due to several factors, such as different nutritional supply, soil composition, irrigation methods and fertilizer spread, or the presence of anthropogenic elements such as ditches, drainage channels, service roads, and in any case, to elements of different nature; all factors that must be understood and managed.

Welch’s test was applied on NDVI average values before, and then on NIR and Red average values, to evaluate the correspondence between the two remote sensing data that nominally work with the same bands, even if at completely different scales of observation. Table 3 and Table 4 show that the two sensors provided statistically similar data, i.e., the mean equality hypothesis ($H_0$) was accepted, for some specific scenarios. However, when the numerical difference of the means was less than τ’ = 0.005, the test application was not necessary because the mean values are indistinguishable. In these cases, the $H_0$ was always accepted. Based on the Welch’s test calculation, the $H_0$ was accepted for June, the NIR data in A, B, C test areas; for August, the NDVI data in A and C test areas, Red data in A test area and NIR data in all test areas. Based on the epochs of acquisition, the analysis of previous results provided much of information [89] for the potential integrated use of the MS2 and S2 sensors.

In June, in dominant bare-ground conditions, the S2 and MS2_10 data could not be considered equivalent. Indeed, despite NIR measurements in A, B and C areas passing the Welch’s test, for Red band and NDVI data the $H_0$ was never accepted.

In July, the tested areas showed nonhomogeneous plant cover conditions due to the corn plants that not yet flowered or that were blooming. Consequently, the bare soil influence is again considerable in the B area and especially in the C area. According to the agronomic expert who supported the entire study, the 5-day phase shift in the two acquisitions in the July stage did not affect the results because the change related to vegetation variation in the field was minimal. Based on this data, it appears that there is a significant difference between the two sensors (p-value < 0.001). Without a total plant cover, the two sensors, looking at the same space, measured it in a different way. Even in the presence of fairly large bare terrain, the S2 data (at 10 m) tend to level off the results, meaning that the most pronounced potential DN differences between pixels are flattened out. This is fixed, on the other hand, by considering proximity data (in this case with original GSD equal to 5 cm) in which, even after aggregation at 10 m, a greater values dynamic is found.

On the other hand, in August, canopy development was almost complete and fairly homogeneous due to the plants reaching the maximum phenological maturity level and to complete the soil coverage. Relating to the NDVI index, the areas A and C passed the test but not the B area. However, in B, the NDVI average value difference of 0.007 is very close to the minimum threshold τ’, i.e., an uncertainty of two bits out of 256, which can be considered negligible. Furthermore, the NIR measurements were similar in the three test areas and for the Red band, the two areas that are statistically different showed differences of less than τ’, not satisfying the applicability conditions of the test itself. Therefore, the results revealed a vegetation in full activity with very strong emission in the NIR and a total absorption in the Red, as expected [47,90].

The knowledge of correspondences and differences between two sensors is important to recognize what Sentinel-2 data can measure with certainty on the one hand and with uncertainty on the other. Consequently, it became useful to understand when to use MS2_10 data as a calibration tool of S2 data and when to use MS2_10 data as an integration of S2 data. In the first case, as in August for corn, the MS2 data become ground truths suitable to calibrate the information extracted from the satellite data [57]; in the second case, as in June for corn, the MS2 data completed and expanded the S2 measurements [89].

Different plant phenological stages gradually define different plant canopy/bare soil ratios. In terms of the relationship between irradiance and reflection, for both sensors, the Red and NIR bands behaved very differently as time varies because they were modulated by the evolution of the phenological phase. The transition from the flowering phase to the ripening phase involves an increase in the leaf surface [65] and a decrease in the bare soil effect, as evidenced by the NDVI and SAVI values. The differences between indices are greater than those of bands. The differences in the Red are almost always greater than in the NIR, especially for periods with a lot of bare soil patches (such as June in our example). Further, Huang et al. [91] also demonstrated that by applying a paired t-test, the differences in the Red channel between different satellite sensors are always highly significant (p < 0.0001). Consequently, reporting the values of Red in the index formulas influences their result. While the reflectance in Red decreased with increasing growth (Table 4), the NIR had the opposite behavior. The vegetation reflects much of the NIR signal.

The use of SAVI [45], in addition to NDVI, derived from the curiosity to test a vegetation index able to minimize the brightness of the soil. Considering the scenario under observation, to make even a temporal comparison between VI products from S2 and MS2_10, the SAVI was calculated specifically for the June and July stages. It did not provide additional information but only a general decrease in values is noted with respect to the NDVI values. The trend of average SAVI do not change, especially for June (A, B, C) and July (B, C), when bare soil is present in a dominant or significant way, even though SAVI was studied to consider and reduce the contribution of the soil beneath the plants. In summary, based on the SAVI, the two sensors see different things, as already noted previously for Red and NIR bands. MAIA-S2, due to its higher spatial resolution, is more suitable to detect soil contribution, and the differences in this index are greater than the threshold τ‘; they are always significant. For these reasons, the authors think that the definition of L coefficient for proximity data should be revised through experimental analyses and described by a specific function.

Concerning linear regression as a tool (cf. 4.3, Figure 5 and Table 5) analyzing the correspondence between S2 and MS2_10 data, its effectiveness was not exhaustive. An initial comment concerns the significant difference between the NDVI ranges found in some cases between S2 values and the MS2_10 values that derive from resampled data. For example, in June, data with a lot of bare soil, in the test area B, most S2 NDVI values are in the 0.18 to 0.33 range (a few up to 0.42). These correspond to lower MS2_10 NDVI values, including between 0.07-0.25 (with a few up to 0.33). Even for area B, in July and in the presence of a variable bare-vegetated soil distribution, the S2 NDVI values were compacted in a much higher and limited range (0.81-0.89) while the corresponding MS2_10 values were distributed in a much wider range (0.58-0.85). This different behavior suggested that MS2_10, despite the resampling, is more able to discriminate different responses linked to mixture patches, while S2 merges the values into a small range [89]. In July, the other test areas also showed a similar trend while, in the August data, the value trend was very different. The data were distributed on NDVI intervals of 2-3 cents for both measurement types. They did not manifest a linear relationship although the estimation process obviously provides the regression line. Even if the statistical tests and the different evaluation between the means indicated a substantial correspondence of the S2 and MS2_10 data, the correlation was very low (for example, in the A and C test areas, R² was 0.28 and 0.47 respectively). Consequently, the value trend seemed to depend more on the measurement variability of the two sensors than on a functional trend.

6. Conclusions

This paper evaluates the integrated use of multispectral data in order to accurately monitor the dynamics of the relationship between bare soil areas and the crop. The protocol developed to compare the two types of data is based on a statistical comparison between Sentinel-2 satellite data and UAV multispectral data, with the same spectral bands.

The activity was carried out in a corn field, considering three stages (June, July and August) and three test areas: A (high vegetation density), B (medium vegetation density), C (low vegetation density). The similarity between the two types of data was accepted:

−

for June, for NIR data in area A, B and C;

−

for August, for NDVI values in test areas A and C, Red data in test area A and NIR data in all test areas.

Based on the acquisition stages, these results provided much information for the potential integrated use of UAV and Sentinel-2 sensors. In particular, the analysis was useful in understanding when to use UAV data as a calibration tool for Sentinel-2 data and when to use UAV data as a supplement to Sentinel-2 data:

−

in August, MS2 data can become ground truths suitable for calibrating information extracted from satellite data;

−

in June, MS2 data can complement and extend S2 measurements.

Scaling up surveys to a local level is very useful when precise or small-area assessments are required. This approach would also allow for parameters that can specialize satellite and proximal data in one or more processing steps, such as spatial adaptation of atmospheric correction algorithms for satellite data, or the possibility of locally characterizing and estimating the SAVI L parameter to account for an actual average percentage of bare soil. All the above considerations will be useful for those who will have to elaborate classifications or detailed analyses aimed at improving the discrimination procedures of crops and bare soils. Finally, the analysis could also be completed in relation to moisture content with the addition of short-wave and thermal sensors.