Evaluation of Agricultural Bare Soil Properties Retrieval from Landsat 8, Sentinel ‐ 2 and PRISMA Satellite Data

: The PRISMA satellite is equipped with an advanced hyperspectral Earth observation technology capable of improving the accuracy of quantitative estimation of bio ‐ geophysical variables in various Earth Science Applications and in particular for soil science. The purpose of this research was to evaluate the ability of the PRISMA hyperspectral imager to estimate topsoil properties (i.e., organic carbon, clay, sand, silt), in comparison with current satellite multispectral sensors. To investigate this expectation, a test was carried out using topsoil data collected in Italy following two approaches. Firstly, PRISMA, Sentinel ‐ 2 and Landsat 8 spectral simulated datasets were obtained from the spectral resampling of a laboratory soil library. Subsequently, bare soil reflectance data were obtained from two experimental areas in Italy, using real satellites images, at dates close to each other. The estimation models of soil properties were calibrated employing both Partial Least Square Regression and Cubist Regression algorithms. The results of the study revealed that the best accuracies in retrieving topsoil properties were obtained by PRISMA data, using both laboratory and real datasets. Indeed, the resampled spectra of the hyperspectral imager provided the best Ratio of Performance to Inter ‐ Quartile distance (RPIQ) for clay (4.87), sand (3.80), and organic carbon (2.59) estimation, for the spectral soil library datasets. For the bare soil reflectance obtained from real satellite imagery, a higher level of prediction accuracy was obtained from PRISMA data, with RPIQ ± SE values of 2.32 ± 0.07 for clay, 3.85 ± 0.19 for silt, and 3.51 ± 0.16 for soil organic carbon. The results for the PRISMA hyperspectral satellite imagery with the Cubist Regression provided the best performance in the prediction of silt, sand, clay and SOC. The same variables were better estimated using PLSR models in the case of the resampled hyperspectral data. The statistical accuracy in the retrieval of SOC from real and resampled PRISMA data revealed the potential of the actual hyperspectral satellite. The results supported the expected good ability of the PRISMA imager to estimate topsoil properties.


Introduction
There is an acknowledged requirement for up-to-date, spatially precise and inexpensive soil-mapping methods to monitor the soil properties of both agronomic and environmental interest and their related processes [1,2]. The spatial and temporal monitoring of soil properties is also highly significant, from both environmental and economic perspectives.
The physical link between soil components and the electromagnetic spectrum has allowed the for development of promising soil spectroscopy techniques for the estimation of soil properties in the laboratory, as well as progressively from airborne and satellite platforms. Satellite hyperspectral images offer an important advantage for monitoring soil properties from the field to the regional scale. The relationships between the soil reflectance spectrum in the optical range and soil mineralogy, chemical composition [3], soil organic matter content [4], soil water content [5], and soil particle size distribution [6] have led to the development of promising data-driven or physical-based methods of estimating soil properties [7]. This is due to the interaction of soil components with visible and infrared radiation that show specific absorption features in this spectral range [8][9][10]. For example, clay minerals have typical spectral features in a shortwave infrared region (SWIR) between 2170 and 2360 nm. Soil organic carbon (SOC) considerably influences the shape and nature of soil reflectance spectra [11]. It has specific spectral features [12] in the visible region with a wide absorption characteristic reported at around 664 nm [13], and in the SWIR region, principally due to the two principal organic compounds that influence the reflectance of lignin (between 1600 and 1800 nm and around 2100 nm) and cellulose (around 2100 nm) [13][14][15]. Moreover, the biochemical components of soil organic matter [2] have effects on the reflectance in the visible (400-700 nm) and NIR-SWIR (700-2500 nm) intervals [2,16] of the spectrum.
Spectral chemometrics models are particularly suitable for the precise quantification of soil properties [17]. These approaches make use of the whole spectrum to assess soil parameters, leading to an optimal selection of the bands and spectral range. The most widely used calibration techniques for developing these models are those relying on the Partial Least Squares Regression (PLSR) [2,8,[18][19][20][21], which has been exceptionally valuable for developing models that comprise a significant number of predictor variables while considering the relationship between spectra and soil characteristics [22]. In order to handle big and spatially varied spectral libraries to estimate soil characteristics over a reduced area, numerous approaches, based mainly on the intensification of the estimation capacity of local samples, were considered [23][24][25]. Cubist Regression and random forest have been progressively utilized to model and predict soil variables at both local and regional scales [22,26].
Numerous soil spectroscopy studies have been carried out in the laboratory and proved effective in the estimation of soil properties, with a significant level of accuracy [9,[27][28][29][30]. Spectral chemometrics for the prediction of soil properties was demonstrated to be relevant and accurate in the simultaneous prediction of several soil variables in the laboratory [31][32][33]. The results of [34] provided good SOC content predictions with PLSR models. The authors of [30] concluded that accurate and globally stable PLSR models for predicting soil variables can be developed from spectra acquired in the laboratory. While this approach provides an estimation of soil properties at specific sampling points, geostatistical or interpolation techniques are essential to produce continuous spatial data, i.e., maps; therefore, the prediction accuracy tend to decrease in larger areas.
The strength of the relationships between spectral features and soil properties decreases from laboratory to satellite imaging spectroscopy [35]. This is mainly due to the combination of weather conditions and atmospheric attenuation, soil conditions, the scale, and the spectral resolution of the sensors. As new challenges arise in the quantification of soil properties, there is a continuing effort to meet these challenges and advance model precision [1]. Many recent scientific papers demonstrated the ability of, and the encouraging results obtained for, soil properties' prediction and mapping from Sentinel-2 and Landsat 8 optical data [31][32][33][34].
The extensive Sentinel-2 and Landsat 8 time series also strengthen the possibility of detecting bare-soil areas in croplands and, accordingly, tracking soil properties across big districts. The respective short and moderate revisit times of the two satellite constellations increase the likelihood of obtaining cloud-free images and a high number of bare soil fields. This is particularly important in the soil imaging spectroscopy context, due to the narrow time window in which it is possible to observe bare soil in croplands. The use of multispectral imagery in small study areas, and considering only bare soils, proved to be successful. However, the quantitative prediction of soil parameters through the bare soil images attained from multispectral remote sensors was not optimized by a deficient spectral resolution [32,33], principally due to the lack of narrow bands in the SWIR region (i.e., the spectral interval influenced by the soil chromophores). This is of great importance for the prediction, for example, of clay, which can be estimated by exploiting welldistinguished spectral features between 2200 and 2400 nm that are related to the O-H stretch and metal-OH bends of the clay lattice of the most clay minerals in soil (e.g., kaolinite, montmorillonite and illite) [36]. These specific spectral features can only be detected by a hyperspectral sensor with an adequate spectral resolution in this region.
In this regard, hyperspectral remote sensing is a suitable technology to produce spatially precise topsoil mapping [15,37]. With the current and forthcoming availability of hyperspectral satellite sensors with an adequate signal to noise ratio (SNR), such as with the launch of PRecursore IperSpettrale della Missione Applicativa (PRISMA) [2,8] and the upcoming launch of EnMAP [1], the mapping of soil variables over large extents using imaging spectroscopy can be achieved. These sensors can regularly upgrade soil maps in bare-soil zones, a task that can only be obtained from existing low spectral resolution sensors at present [38]. Furthermore, several bands in the SWIR spectral interval of the forthcoming satellite hyperspectral imagers will presumably permit a more accurate prediction of soil parameters in comparison to current multispectral sensors.
There is an increasing number of papers concerning the estimation of soil variables exploiting resampled and/or simulated satellite hyperspectral data using statistical models [8,11,39]. Gomez et al. [11], while comparing prediction models of SOC using hyperspectral proximal and remote sensing data, revealed that the SOC map predicted using PLSR models showed a similarity with field observations, demonstrating the potential of combining hyperspectral remote sensing datasets and chemometrics models. Novel approaches to soil information retrieval have been proposed in the literature, based on simulated data from satellite hyperspectral sensors, which have already been launched or are to be launched soon. Castaldi et al. [21] showed that using simulated PRISMA datasets with an a priori recognition of the soil humidity classes can lead to a reduction in the error in clay estimation. Steinberg et al. [40] showed that simulated EnMAP hyperspectral satellite data allows for the prediction of iron oxide, clay, and SOC in topsoil with a minor reduction in soil estimate accuracy compared to the airborne scale, while still retaining, with PLSR, a satisfactory prediction accuracy. Castaldi et al. [20] investigated the abilities of seven multispectral and simulated hyperspectral satellite datasets to predict soil variables from PLSR models. Their results showed that, in the absence and with the addition of noise, the performances of hyperspectral data, presented by Ratio of Performance to Inter-Quartile (RPIQ) range, were significantly better compared to those of multispectral imagers for topsoil clay, sand, silt and SOC estimation. The authors of [2] obtained a reasonable SOC prediction accuracy from a local PLSR approach applied to simulated EnMAP satellite imagery.
Even though previous studies explored the abilities of actual PRISMA data in nonphotosynthetic vegetation monitoring [41] and land-use monitoring [42], the ability of the PRISMA hyperspectral sensor in the soil domain has not been reported to date.
The objective of this study is to assess the potential of the PRISMA hyperspectral data for the prediction of topsoil properties over agricultural bare soil areas, in comparison with Sentinel-2 and Landsat 8 multispectral data. We selected two experimental agricultural areas in Italy and compared different processing approaches using both the soil spectral measurements elaborated in the laboratory and extracted from satellite imagers at the field level. For actual remote-sensing data, bare soil fields were identified based on pre-defined spectral indexes thresholds. Subsequently, the usefulness of spectra extracted from single date images was assessed.

Materials and Methods
The flowchart of the methodology is illustrated in Figure 1 and, hereafter, the following steps are described. These include: (A) soil sample collection and analysis (Section 2.2), (B) spectral laboratory measurements (Section 2.3); (C) pre-processing of the hyperspectral (Section 2.4.1) and multispectral (Section 2.4.2) remote sensing satellite data and (C) the prediction of soil properties based on bare soil spectra from both laboratory resampled data (Section 2.5.1) and collected satellite data (Section 2.5.2).

Study Area
Two agricultural test sites were selected for the estimation of topsoil properties in Central and Southern Italy.
For Central Italy, the area under investigation is located within the Maccarese S.p.A. farm (latitude 41°52′18" N, longitude 12°14′05" E, altitude 8 m a.s.l.) near Rome ( Figure  2a). The study area underwent reclamation works in the early 1920s, which led to considerable changes to the landscape of the coastal marshland. The soil of this area is classified, according to a recent regional soil map [43], mainly as Dystric Arenosols (dunes and rear dune and recent alluvial and aeolian deposits) and Eutric Endochromic Luvisols (coastal plain reclaimed on fluvio-marsh sediments and recent fill), using the WRB FAO soil classification [43]. The Maccarese farm area is characterized by groups of beach ridges corresponding to eight homogeneous complexes that can be followed in both the Northern and in Southern part of the delta of the river Tevere and coastline dynamics. In particular, the soil parent materials of the Northeastern part of the studied site are alluvial deposits of the Arrone river. The main soil texture is sandy clay loam, which becomes more clayey towards the Northeast of the site; consequently, no carbonate topsoil is present in the area, while, on dunes, shell remnants can be present. The principal crops grown in the area are maize, followed by durum wheat, winter wheat, fava bean and forage crops.
The investigated area in Southern Italy is Pignola (Figure 2b) in the Basilicata Region (Latitude 40°33′45" N, Longitude 15°45′40" E). The sampled area altitudes vary from 720 to 760 m a.s.l. According to the Regional pedagogical map [44], which follows the WRB FAO soil classification, the Pignola soils correspond to Eutri-Vertic Cambisols (fluvio-lacustrine deposits of a local basin). The site is generally characterized by a moderately differentiated profile due to the partial removal of carbonates and browning of the horizons with free drainage. The main crop classes consist of arable land, pastures, and coppice woodland.

Field and Soil Sampling
Field campaigns were performed in both sites: for Maccarese on 11 November 2019, 28 January 2020 and 17 February 2020; for Pignola on 20 January 2020. The same protocol was followed, collecting soil samples located according to an Elementary Sample Unit (ESU) scheme set up to fit the PRISMA spatial resolution (Table 2 Table 1), i.e., according to 30 by 30 m quadrats. Each ESU contained 10 sampling points, positioned 15 m apart. On each point, three soil samples were collected within a 1 m radius, with an Ejkelkamp auger in the 0 -10 cm depth layer. To optimize the sampling strategy on the PRISMA pixel size, the samples collected in the 10 ESU points were subsequently processed in the lab with two different sets of analyses. Soil samples from points 1 to 5 and samples from 6 to 9 were kept separate at the time of the laboratory analysis, thus obtaining two average soil properties values per ESU. The data of the points from 6 to 10 could be used as a 15 by 15 m sub-ESU, suitable for higher resolution sensors, e.g., compatible with the Sentinel-2 spatial scale, and to evaluate the variability within the ESU. We collected soil samples from a total of 61 ESU (43 from Maccarese and 18 from Pignola). Soil samples were analyzed to determine clay, silt, and sand content, using the pipette method, according to the USDA textural thresholds, and SOC content using the Walkley-Black method. Table 1 shows the descriptive statistics of soil properties in the ESU and sub-ESU of the two experimental locations.

Laboratory Dataset
After air-drying, crushing, mixing, and sieving at 2 mm, a subsample of dry sieved soil was extracted from all soil samples that were originally collected in the field. The subsamples were placed in labeled small aluminum bowls, previously painted in black, on the interior. Reflectance measurements were performed across the 350 -2500 nm range, using an Analytical Spectral Device (ASD, SN:6255) Field Spec Fr Pro spectroradiometer (ASD Inc., Denver, Boulder, USA), with a spectral resolution of 3 nm in the Vis-NIR range (350 -1050 nm) and of 10 nm in the SWIR range (1050 -2500 nm). The ASD was equipped with a contact probe containing a high-intensity quartz-halogen lamp (50 W) ( Figure 3). Radiance across the 350 -2500 nm spectral range was measured and converted to reflectance using a calibrated white Spectralon (Labsphere, Inc., USA) reference panel. The protocol used for the measurements accurately followed the approach detailed in [28] and [1].
To prepare and control the experimental conditions, the ASD was turned on at least 60 minutes before the readings started, to warm up the spectrometers and the lamps. The room was kept dark during the readings to avoid any interference with the instrument. To prepare the soil sample readings, the following operations were performed: (i) mixing and flattening the soil within a bowl with a glass surface, (ii) bringing the samples to the ASD contact probe (CP) using an elevator, while the CP was held firmly in place, (iii) putting the soil in contact with the CP and perform the reading, and (iv) repeating the first three steps a total of 5 times per soil sample. Figure 3 illustrates the main phases of these operations. To correctly process the measurements, a check was performed every five sample readings to assess whether the white Spectralon panel had returned to 100%. If not, the white reference was reacquired to set it back to 100% to correct the detector's drift. The soil standard samples of Lucky Bay and Wylie Bay [28] were also reacquired to standardize the ASD readings according to the standard protocol [28], to minimize the systematic effect of the measurements. The standardization followed the protocol set-up of Ben Dor et al. [28]. The five spectral replicates for each soil sample were averaged, and each spectrum was associated with the sample ID and its physical and chemical properties, constituting the laboratory dataset.

PRISMA Hyperspectral Data
The "PRecursore IperSpettrale della Missione Applicativa" (PRISMA), put in orbit from 2019 by the Italian Space Agency (ASI), is a technology demonstrator satellite mission finalized to develop user-driven [45] pre-operational products, taking advantage of the combination of a hyperspectral sensor and a panchromatic imagery. PRISMA mission is fully funded by ASI, while the payload was developed by the Italian Leonardo company. The main mission characteristics are described in the work of [46,47] and are briefly listed in Table 2. PRISMA payload is a push-broom design based on a prism-based concept, acquiring 234 spectral bands in the 400 -2500 nm spectral range. The first spectrometer provides 66 bands in the Vis-NIR spectral range, while the second spectrometer operates in the SWIR. The sun synchronous orbit obtains a swath of 30 km images with a Ground Sampling Distance (GSD) of 30 m, while PAN provides a co-registered image at 5 m/pixel. The revisit time is 29 days (from nadir), while shorter revisit times, of up to 7 days, could be managed by the available re-look capability. Mission lifetime is 5 years. The use of a prism as diffraction elements determines a non-uniform bandwidth across the Focal Plane Array (FPA). Moreover, as the center wavelength (CW) could change due to FPAs temperature each image is correlated by ancillary information indicating the proper CW and FWHM. Images are provided by ASI (https://prisma.asi.it/) in the Hierarchical Data Format -Earth Observing System (HDF-EOS5) format at different processing levels: L1 TOA radiance units, L2C ground reflectance without geometric correction, and L2D, corresponding to geocoded ground reflectance. Details of the ASI standard processing are available in the PRISMA Products Specification Document (ASI, PRISMA Products Specification Document Issue 2.1, accessed on 22/11/2021). The L2D processing level was selected for all the images utilized for the scope of this work.
The results presented in this study are based on the L2D (geocoded version of the level 2C "At-surface Reflectance Product"; [47,48]) data current version (i.e., 2.0.4) of PRISMA products distributed to the registered users by ASI. The requested L2D images (PRS_L2D_STD) were downloaded from the PRISMA mission catalog website (http://prisma-i.it). The PRISMA images used for this study were collected on the two test sites between November 2019 and August 2021 (see Table 3) with observer (roll ±) angles varying from 4° to 18°. The overall weather conditions were stable, i.e., clear-sky conditions, during the days of acquisition (PRISMA image cloud coverage between 0.01% and 7.14%) with a very low aerosol presence during the field surveys. The original HDF-EOS5 files were converted to Band Interleaved by Line (BIL) file format using an in-house tool developed in the IDL/ENVI 8.7.2. (L3Harris Technologies, USA) by CNR IMAA (Italy) researchers. Before merging the two extracted Vis-NIR and SWIR hyperspectral datasets into a single image file, it was necessary to remove the overlapping reflectance bands in the 930 -998 nm spectral range and remove the bands with a low SNR or those affected by noise.
Moreover, even though L2D PRISMA images are supplied and geocoded (in our case with Datum: WGS-84 and Projection: UTM 33 N), there is a residual, non-constant small shift (maximum 5 pixels) in all the images when compared to the ancillary digital cartography and vectorized shape file of the fields of interest. Therefore, we resampled all the images on the same grid of the field corners vectors, using a Nearest Neighbor algorithm implemented in the THOR Change Detection tool in ENVI 5.5.2. (L3Harris Technologies, USA). This was performed to co-register and align all the images (a) to set-up the data stack layer for the subsequent multitemporal analysis and (b) to allow for the accurate extraction of PRISMA spectra in correspondence to the ground-collected field soil samples.

Sentinel-2 and Landsat 8 Multispectral Data
Sentinel-2 Multi-Spectral Instrument (S2/MSI) and Landsat 8 Operational Land Imager (L8/OLI) data were obtained within the Google Earth Engine (GEE) environment for the two study areas, at dates very close to the PRISMA acquisitions ( Table 3). The images were selected from the S2/MSI level 2A (COPERNICUS/S2_SR in GEE) and L8/OLI level 2, collection 2, Tier 1 (LANDSAT/LC08/C02/T1_L2 in GEE) image collections. Sentinel-2 pixels affected by clouds and clouds' shadows were masked using the 'QA60' bitmask band for cirrus and opaque clouds, the 'MSK_CLDPRB' for clouds, and the 'MSK_SNWPRB' for snow. The cloudy Landsat 8 pixels were removed using the Pixel Quality Assessment Band. We

Soil Properties' Estimation
Soil property estimation models were built, starting from SOC, clay, sand, and silt content (dependent variables or outcomes) measured by chemical and physical analysis of the samples collected in the Maccarese and Pignola and the spectral data acquired in laboratory conditions or by satellite sensors (PRISMA, S2/MSI and L8/OLI), representing the independent variables or predictors. All estimation models were calibrated by testing two algorithms: Partial Least Square Regression (PLSR) [49] and Cubist [50]. The Cubist is a rule-based regression algorithm that allows a specific linear regression model to be defined for a subset of the calibration dataset. Therefore, the Cubist algorithm makes a partition at each node, characterized by a condition (rule) involving the independent variables, i.e., the soil spectral bands.
To reduce the noise of the actual PRISMA data and enhance the spectral features related to soil properties, the following spectral pre-treatments were tested: Savitzy-Golay filter (SG) [51], standard normal variate (SNV) [52] and continuum removal (CR) [36]. The model parameter tuning was conducted set up using the caret package in R [53], testing different numbers of components for PLSR and different combination of committees and neighbors for Cubist models. The best model was determined according to the lowest 10fold cross-validation Root Mean Square Error (RMSE, Eq. 1)

∑
(1) where ym are the soil properties' measured values and yp are the values predicted by the models. Each observation (measured value) appears in only one of the 10 folds selected for the cross-validation process. Thus, no overlap exists between groups, and the observations were split according to a pseudorandom number generator that allows for the output of the random fold-splitting to be reproduced to cross-validate all the predictive models.
To compare the estimation accuracy between soil properties and the other literature results, the coefficient of determination (R 2 , Eq. 2), Ratio of Performance to Deviation (RPD, Eq. 3) and Ratio of Performance to Inter-Quartile distance [54] (RPIQ, Eq. 4) were computed where RSS is the sum of squares of the residuals, TSS is the total sum of squares, and Std and IQ are the standard deviation and the observed values, respectively. For each combination of outcome and predictors (e.g., SOC estimation using S2/MSI data), only the best model in terms of average 10-fold cross-validation Cubist and PLSR is reported in the manuscript. Moreover, in order to better compare the prediction accuracy obtained by bare soil satellite spectra, the standard error of the ten cross-validation folds is also provided for RPIQ and RPD values.

Laboratory and Resampled Spectra
The spectral library acquired by ASD Fieldspec Pro spectroradiometer in our laboratory under controlled conditions according to the Internal Soil Standard (ISS) [28] protocol was included in this work, as it might represent the best conditions for the estimation of soil properties, especially due to the very short distance from the sensor (contact probe) to the target object. Thus, the laboratory spectra were applied to evaluate the spectral resolution ability of the PRISMA, S2/MSI and L8/OLI sensors in ideal conditions, i.e., without disturbing factors affecting the satellite signal [31].
Four different spectral datasets were retrieved from the laboratory spectral library, and as many estimations models as soil properties (SOC, clay, sand, and silt): (i) all the spectral data provided by the ASD Fieldspec Pro spectroradiometer (hereinafter referred to as LAB), (ii) ASD spectra resampled according the PRISMA bands, excluding noisy bands, as described in Section 2.4.1 (hereinafter referred to as LAB PRISMA), (iii) ASD spectra resampled according to the S2/MSI bands (hereinafter referred to as LAB S2), and (iv) ASD spectra resampled according to the L8 bands (hereinafter referred to as LAB L8).
The resampling process was carried out with a Gaussian model, considering only the central wavelength and the bandwidth of each satellite band, using the spectral resampling function of the package hsdr in R [55]. This procedure provided a common method for all the satellites since a spectral response function was unavailable for PRISMA.

Bare Soil Satellite Spectra
For each soil sample collected in the Maccarese and Pignola areas, we selected the PRISMA image, among those listed as calibration in Table 3, that corresponded to the bare-soil conditions for the sampling areas. The selection was carried out by avoiding clouds, green and dry vegetation, thanks to the (a) agricultural management information, mainly provided by direct observations in the field, (b) visual check of the satellite images, (c) computation of the NDVI using a threshold of 0.35 to remove green vegetation pixels, and (d) the burn index 2 (NBR2; Eq. 5) to detect and remove non-photosynthetic vegetation (NPV) for S2/MSI and L8/OLI images 2 where is the reflectance value of B11 for S2/MSI (B11) and B6 for L8/OLI instrument, while corresponds to B12 for S2 and B7 for L8. Instead, we used the normalized cellulose absorption index (nCAI; Eq. 6) (nCAI < 0.03) for the PRIMA hyperspectral data. Compared to NBR2, this narrow band index is better suited to exploiting the of hyperspectral data's ability to identify spectral features related to dry matter and cellulose [2,56,57] and, therefore, to detect NPV in the field. 0.5 * 0.5 * The bare-soil spectra were extracted according to the selected dates at different sampling point locations. For this study, we used the S2/MSI and L8/OLI images with the closest acquisition dates to the PRISMA dataset (Table 3).

Spatial Assessment of Topsoil Maps
An assessment of the correspondence of estimated maps of SOC, clay, sand, and silt, with the pattern of spatial soil variability, was performed for an agricultural field within the Maccarese area (referred to as "Field B071" in Figure 2a), for which ground data on soil's apparent electrical resistivity were available. The best models previously obtained for each soil property and satellite sensor were applied to independent satellite images, i.e., not those used for the calibration of the models (Table 3). For this purpose, we selected satellite images acquired in mid-spring 2021 for S2/MSI (15 h May) and PRISMA (17 th May), while, for L8/OLI, no cloud-free images were acquired in the same time window; consequently, we selected an image acquired at the same time of the year (mid-spring) in 2018, showing similar field conditions. Each soil properties predictions map was spatially compared with the map obtained from a geoelectrical survey carried out on the field of the Maccarese farm area, since some soil properties, such as clay and SOC, are related to electrical resistivity/conductivity [58][59][60][61]. The geological survey acquired an apparent electrical resistivity (Ω m) of the 0 -0.5 m topsoil layer at the meter resolution, using an Automatic Resistivity Profiling (ARP) multi-electrode continuous profiler, pulled by a quad. More specific information on the geoelectrical survey is provided by Casa et al. [62]. Resistivity maps were obtained by spatial interpolation of the ARP measurements using block kriging with block dimensions corresponding to the resolutions of PRISMA and Landsat (30 m) and Sentinel-2 (10 m). The satellite-derived maps of predictions were used as underlying grids for the interpolation. Figure 4 shows the ARP map obtained at Sentinel-2 resolution. The spatial correlation between the resistivity maps and the maps of satellite-predicted soil properties was assessed using the Dutilleul t-test for the spatial processes [63,64] of the SpatialPack package implemented in R [65]. This test allows for an evaluation of the consistency of satellite-predicted maps, with the spatial patterns highlighted by the geophysical survey.

Soil Properties Estimation Using Laboratory and Resampled Spectra
The accuracy of the models resampled using laboratory spectra (LAB) was very high for clay (RPIQ: 4.46), sand (RPIQ: 3.79) and silt (RPIQ: 5.00), and lower for SOC (RPIQ: 2.59; RPD: 1.81; R 2 : 0.75) ( Table 4). The PLSR technique, in the case of the laboratory spectra, provided the best models for the retrieval of the four soil properties.
By resampling the laboratory spectra according to the multispectral Sentinel-2 (LAB S2) and Landsat 8 (LAB L8) bands, the accuracy decreased for all the soil properties and particularly for sand in the case of LAB S2 (RMSE: 7.76; RPIQ: 3.32) and SOC in the case of LAB L8 (RMSE: 0.35; RPIQ: 1.70; RPD: 1.22; R 2 : 0.50). The models resampled to the Sentinel-2 data showed a good estimation accuracy for silt (RPIQ: 4.56; RPD: 3.08; R 2 : 0.92). The same property estimated using the LAB L8 data also provided a good degree of accuracy (RPIQ: 3.54; RPD: 2.39; R 2 : 0.86). Otherwise, LAB S2 data showed the best performances among multispectral datasets for all variables. For the LAB S2 and LAB L8 multispectral data, the best prediction accuracies were obtained using the Cubist method, except for SOC, where the PLSR provided the best results for the LAB S2.

Soil Properties Estimation Using Bare Soil Satellite Spectra
As expected, the estimation accuracy with real satellite data decreased for all the soil properties as compared to the laboratory resampled spectra ( Table 5). The worsening was particularly apparent for clay. However, the SOC model from PRISMA presented even better results than those obtained from LAB PRISMA data (RPIQ ± SE: 3.51 ± 0.16; RPD ± SE: 2.43 ± 0.11 for PRISMA and RPIQ: 2.59; RPD: 1.81 for LAB PRISMA).
The difference in statistical accuracy between hyperspectral and multispectral datasets became more apparent when using real satellite data, as compared with resampled laboratory data. The statistical accuracy obtained using the multispectral imagery was lower compared with those obtained using the hyperspectral image for all soil variables, except for sand, where the Sentinel-2 estimation model provided a slightly higher prediction accuracy (RPIQ ± SE: 2.21 ± 0.10; R 2 : 0.79) in comparison with PRISMA (RPIQ ± SE: 2.15 ± 0.05; R 2 : 0.77). The soil property models obtained using the Landsat 8 image provided a sufficient degree of accuracy (RPIQ ± SE: 2.01 ± 0.08; RPD ± SE: 1.71 ± 0.07 in the case of clay); however, all the L8/OLI statistics were worse than those obtained by S2/MSI data. Table 5. Estimation accuracy of the best-performing prediction models for soil properties derived from Figure 2 and Landsat 8) spectra of bare soil. RMSE: Root Mean Square Error; R 2 : coefficient of determination; RPD: Ratio of Performance to Deviation; RPIQ: Ratio of Performance to Inter-Quartile distance; SE: standard error.

Variable
Satellite  The RPD values obtained from hyperspectral and multispectral resampled data (Figure 6) for clay (RPD: 4.12 for LAB PRISMA and 3.37 for LAB S2) and sand (RPD: 3.58 for LAB PRISMA and 3.13 for LAB S2) estimation were significantly higher than those obtained with real spectra (case of clay: RPD is equal to 1.94 for PRISMA and 1.87 for Sentinel-2; case of sand: RPD is equal to 2.01 for PRISMA and 2.05 for Sentinel-2). Conversely, the estimation accuracy of SOC for the hyperspectral imager (RPD: 2.43) was statistically higher compared with the other laboratory (RPD: 1.81), resampled (LAB PRISMA RPD: 1.81) and real (Sentinel-2 RPD: 1.52) datasets.
In the absence of noise, and considering only the spectral characteristics, the performance of resampled hyperspectral data was significantly better than those of the resampled multispectral datasets. The PRISMA hyperspectral satellite data display a very similar behavior with respect to the laboratory spectral data when estimating all soil texture and SOC content variables (case of sand and SOC: RPD is equal to 3.58 and 1.81 for both LAB and LAB PRISMA, respectively).
For the SOC, clay and silt soil variables, the RPD confirmed the better estimation accuracy of the hyperspectral imager compared to the multispectral imagers. An exception was observed for sand estimation with Sentinel-2 data, which showed slightly higher results (RPD: 2.05) than those achieved using PRISMA (RPD: 2.01). In more detail, we tested the application of the PLSR and Cubist Regression to laboratory resampled and real multispectral and hyperspectral satellite data to quantify the soil texture and SOC content. The results of the use of PRISMA hyperspectral satellite imagery with the Cubist Regression model provided the best-performing models for the prediction of clay, silt, sand and SOC. The variables silt, sand and SOC were better estimated using PLSR models in the case of the multispectral data. In this case, PLSR models performed better than those obtained using the Cubist Regression algorithm, which nevertheless produced the best soil texture and SOC maps when using resampled multispectral data. The PLSR algorithm provided the best predictor for use as an algorithm for soil texture and SOC content mapping over the study area using both the hyperspectral LAB and the resampled (PRISMA) datasets.
The best validated models using the different approaches (Cubist Regression and PLSR) were applied to the independent mapping images ( Table 3) Table 6 presents an analysis of the spatial correlation between the resistivity map and the estimated soil properties. Dutilleul's correlation between the ARP map and the estimated four soil properties from the hyperspectral sensor presented values of 0.58 for sand, -0.51 for silt, and -0.67 for clay. Clay and sand satellite-predicted maps presented a higher spatial correlation with the electrical resistivity when derived from PRISMA (Table 6) compared to those derived from S2/MSI. However, using L8/OLI, the spatial relationship was slightly higher for the case of sand-mapping results with a Dutilleul's correlation of 0.65. The highest predicted sand contents, as well as medium predicted values at the center of the field (Figure 7 d,f), highlighted a similar pattern for Landsat 8 and PRISMA platforms. Considering our results, the SOC map derived using the S2/MSI shows the highest spatial correlation, while a moderate correlation was obtained for SOC mapped from PRISMA data. Both multispectral (S2/MSI) and hyperspectral (PRISMA) sensor-derived SOC maps showed high spatial relationships with the geophysical survey results.
Comparing the spatial distribution maps of SOC derived from both spaceborne sensors, data showed a similar pattern (Figure 7 l,m).

Discussion
The results obtained in this study are a first attempt to assess the level of accuracy of the newly launched hyperspectral imager PRISMA for the estimation and mapping of bare-soil properties. For this purpose, PRISMA, S2/MSI and L8/OLI actual bare-soil spectra acquired in agricultural areas were used to build clay, silt, sand and SOC estimation models, and the three sensors were compared in terms of model prediction accuracy.
Topsoil clay content can generally be monitored by its mineralogical abundance and granulometry, and a high spectral resolution can be essential for detecting its specific wavelengths in the SWIR spectral region [29]. In this regard, the characteristic absorption features linked to the clay lattice are diverse, according to the composition of clay minerals in the soil, e.g., the main absorption features of kaolinite and montmorillonite are centered at 1400, 1900, and 2200 nm [66], while illite has two bands centered near 2300 and 2400 nm [20]. Therefore, as expected, as the S2/MSI and L8/OLI have both only two wide bands in the SWIR, their clay prediction accuracy is lower than that observed with laboratory resampled and real PRISMA imagery data (Tables 4 and 5). In their research, Fongaro et al. [67] confirmed that spectral sensors with reduced bands, such as IKONOS, have low outcomes for clay estimate, as shown by the lower R 2 values. Castaldi [31] confirmed that neither S2/MSI nor L8/OLI could precisely predict the clay content, because of the low spectral resolution and large bandwidth (between 85 and 187 nm) of their SWIR bands, which prevent the exploitation of the narrow spectral features related to clay. The higher spectral resolution of PRISMA imagery could thus offer improvements in the assessment and mapping of clay content as compared to multispectral sensors. However, the improvement shown by PRISMA in this work, in terms of clay estimation, is more remarkable for resampled data (Table 4) than for actual satellite data (Table 5; Figure 6). This could be due to the lower SNR of the hyperspectral sensors in comparison to broad-band multispectral sensors. This SNR lowering is mainly due to energy accumulation in the narrow spectral bands of the hyperspectral sensors [20]. Figure 8 shows the percentage of times where each PRISMA band was used in a condition and/or a linear model of the Cubist algorithm for the estimation of soil properties. The clay prediction model used only visible and SWIR bands (Figure 8a), and wavelengths between 2100 and 2200 nm were largely used for both models and conditions, as expected, for the presence of the clay mineral spectral features in this spectral region. Table 4 showed how LAB PRISMA spectra also provided the best results for sand and silt estimation, very close to those provided by LAB spectra. However, for clay, when moving from laboratory to actual satellite data, the differences between PRISMA and the two multispectral sensors became smaller (Table 5). For actual PRISMA data models, sand and silt showed an increased involvement of the spectral bands as compared to clay, for both Cubist conditions and models (Figures 8b and c). However, the visible region is crucial for all three soil texture properties. This is mainly due to the existing correlation between hue of soil and clay and sand percentage [20]. Generally, under equal moisture conditions and parental material, soil brightness (high reflectance values) is related to high sand percentage, while a dark hue (low reflectance values) can be related to high clay content. Fongaro et al. [67] used the Cubist model to test the bands' importance for mapping soil attributes and highlighted the important contribution of the visible bands to sand estimation, as sand content is strongly related to the high albedo of quartz. In this regard, the study of [30] found a strong increase in reflectance in the NIR band of the Landsat 5 sensor for sandy soil, and this outcome is confirmed by Figure 8b, which shows how the NIR region for the PRISMA sensor is important for both Cubist conditions and models.
The SOC also affects the soil spectra in the visible range in the same direction of clay content, and Figure 8d clearly confirms the importance of the visible region for SOC estimation. In this regard, Reference [68] detected strong correlations between organic matter and reflectance, measured at 450, 590 and 664 nm. Shi et al. [18] explored the possibility of using hyperspectral remote sensing imagery to rapidly produce an SOC map at the regional scale, observing high variance importance projection (VIP) values in their PLSR model in the visible region between 400 and 600 nm. However, in Figure 8d, NIR and SWIR regions are also of great importance for the prediction; this is mainly due to the presence of overtones related to lignin between 1600 and 1800 nm and cellulose at around 2100 nm [68,69]. Moreover, other organic compounds of the organic matter, such as amide and aliphatic groups, can affect the spectral response at around 2300 nm.
The high importance of the SWIR bands in the prediction models, especially for SOC and clay, confirmed the good quality of the PRISMA sensor in this spectral region. This is an important improvement over the previous satellite hyperspectral sensors, particularly Hyperion, which showed strong limitations for clay content estimation due to their very low signal to noise ratio (SNR) in the SWIR, and especially at around 2200 nm, where the clay mineral features are located [20]. In our results, the best SOC estimation accuracy was obtained using PRISMA satellite imagery. PRISMA images showed the highest potential in terms of SOC prediction ability, even better than laboratory tests. This is probably a consequence of the conditions in which the image was acquired, where the pattern of soil moisture in the field may have played a role in the prediction of the SOC. In fact, spectral data collected under laboratory conditions were only computed under dry sample conditions. Casa et al. [39] showed that the presence of the spatial variability of soil moisture, at the time of remote sensing image acquisition, could improve the estimation of soil properties as compared to uniform dry soil conditions. While, for clay and sand, the difference in accuracy between PRISMA and the other sensors is not too evident, the hyperspectral sensor provided a strongly improved estimation accuracy for SOC (RPIQ: 3.51 ± 0.16) as compared to S2/MSI (RPIQ: 1.8± 0.20) and L8/OLI (1.98 ± 0.12). It should be noted that this difference is even more evident for high SOC values (> 1.5%), which Sentinel-2 and the Landsat 8 model underestimated (Figure 5d).
Concerning the spatial distribution of the obtained soil maps, the analysis of spatial correlation between the resistivity data and the estimated soil properties showed a significant correlation with the geophysical survey for all the PRISMA maps, although the Landsat 8 and Sentinel-2 outperformed PRISMA for sand and SOC correlations, respectively (Table 6). Sentinel-2 SOC maps ( Figure 7m) has a good inverse correlation with resistivity, probably due to the higher spatial resolution (10/20 m) in comparison with PRISMA and Landsat 8. However, the SOC patterns shown in Figure 7m are very similar to those obtained by PRISMA at a spatial resolution of 30 m (Figure 7l). Consequently, the spatial resolution of the PRISMA sensor does not seem to be a limiting factor for the estimation of soil properties for agricultural areas, especially where the spatial variability of soil properties does not occur within a very short range.
The high spectral resolution in the SWIR region of the PRISMA sensor allowed for the estimation of topsoil properties from bare soil images within the studied agricultural area; however, an investigation of the accuracy at other sites should be conducted. Future work will include different test areas in different soil regions to test the ability of PRISMA hyperspectral data to retrieve soil texture, SOC and other topsoil properties under different conditions (e.g., climate, soil types and soil moisture). The application of the PRISMA hyperspectral satellite to broad-scale mapping would be of great importance, considering other environmental covariates for a more robust modelling [70].

Conclusions
This research considers the ability of the hyperspectral sensor PRISMA in comparison with Sentinel-2 and Landsat 8 to predict and map soil texture and SOC content in two agricultural areas in Italy.
We obtained high prediction accuracies using simulated hyperspectral PRISMA datasets (RPIQ: 4.87 and 3.80 for clay and sand, respectively). The results were very similar to those of the full spectra (LAB) datasets, while a slight worsening was observed in terms of accuracy using spectra extracted from the real hyperspectral and multispectral spaceborne sensors.
A comparison between a resistivity map on a field of the study area and the soil texture and SOC maps retrieved from the best models applied to the PRISMA imagery revealed that the predicted maps are consistent with the patterns of variation in soil prop-erties that were observed in the field, as confirmed by the high spatial correlation coefficients, especially for the case of clay and silt. Consequently, the medium spatial resolution of the PRISMA images (30 m) is not a limiting factor for mapping soil properties in the investigated area.
Our results indicate that PRISMA satellite hyperspectral data improved, in relative terms, the accuracy of soil variables' estimation from bare-soil imagery, as compared to existing multispectral sensors with few and broad bands in the SWIR region. The higher spectral resolution of the hyperspectral imagery can exploit the well-defined narrow absorption spectral features characteristically linked to the functional sets in the Vis-NIR and SWIR spectral ranges associated with clay minerals and most important compounds, which constitute the organic matter in the soil. We conclude that hyperspectral data from the available PRISMA and the forthcoming satellite missions (e.g., EnMAP and CHIME) can advance the mapping and monitoring of soil texture and SOC as compared to the currently available imagers. Funding: This research was partly funded by the European Space Agency (ESA) within the project "Copernicus Hyperspectral Imaging Mission CHIME -Mission Requirements Consolidation", grant number N°4000125506/18/NL/IA. This study was co-funded by the Italian Space Agency with PRISCAV projects, grant number N°2019-5-HH.0. This work was also supported by the European Union's Horizon H2020 research and innovation European Joint Programme Cofound on Agricultural Soil Management (EJP-SOIL grant number 862695) and was carried out in the framework of the STEROPES of EJP-SOIL.

Data Availability Statement:
The data presented in this study will be available under request.