Evaluation framework for the generation of continental bare surface reflectance composites

Contents lists available at ScienceDirect Geoderma journal homepage: www.elsevier.com/locate/geoderma Evaluation framework for the generation of continental bare surface reflectance composites Paul Karlshoefera,∗, Pablo d’Angelo a, Jonas Eberle b, Uta Heiden a aGerman Aerospace Center, The Remote Sensing Technology Institute, Münchener Str. 20, Weßling, 82234, Germany bGerman Aerospace Center, Earth Observation Center, Münchener Str. 20, Weßling, 82234, Germany A R T I C L E I N F O Handling Editor: Budiman Minasny Keywords: Bare soil Spectral disturbances Temporal composite Regionalization Europe Evaluation Spectral index Threshold definition Cloud cover SCMaP A B S T R A C T Soils play a pivotal role in supporting ecosystems, human health, food security, and climate regulation. Since several years, temporal composites of bare soil reflectances derived from multispectral satellite data are used as input for soil property modeling. Due to the importance of these model inputs, the quality of the surface reflectance composites (SRC) is essential. The quality depends on the precise selection of pixels that are free of green and dry vegetation, cloud contamination and other atmospheric disturbances. However, there is a lack of suitable concepts and tools to evaluate the impact of the diverse processing parameters for the generation of SRC, especially for large areas such as continents. This study presents a novel approach to evaluate the process of computing bare SRCs across large geographical areas. It can estimate the theoretical limit achievable with defined processing parameters (spectral indices, thresholds, specific filtering, etc.) and it is also suitable to compare the performance of different SRC concepts from the literature. The performance is derived from the angular spectral distance between reference spectra derived from the LUCAS survey and the SRC spectra. It is demonstrated that a linear combination of two spectral indices complemented with a regional threshold dataset keep the complexity of threshold data sets low while performing well across Europe. The results also show that regionalization is as crucial as the choice of the index itself. The additional outlier removal focusing on clouds and haze marginally improved the SRC at the continental scale but can be very effective for areas with more frequent clouds. The proposed method offers two main advantages. First, it allows for parameter customization tailored to the region of interest, or, at minimum, to areas well represented by the reference data. Second, it facilitates the systematic evaluation of successive adaptations in the SRC generation process, eliminating the labor-intensive and error-prone task of visually comparing images to assess improvements in the SRC final product. The final bare surface reflectance composite for Europe and adjacent regions provids a robust foundation for future large-scale soil and bare surface monitoring. 1. Introduction The properties of bare soils play a pivotal role in supporting ecosystems, human health, food security, and climate regulation (Montanarella et al., 2015). Earth Observation can be used to directly model bare soil characteristics based on their spectral reflectance characteristics (Chabrillat et al., 2002; Ben-Dor et al., 2009; Gerighausen et al., 2012; Ward et al., 2020) or based on environmental covariates (Mulder et al., 2011; Minasny and McBratney, 2016; Poggio et al., 2021). Recently, temporal composites of bare soil reflectances derived from space-borne satellite data over multiple seasons are used as input for soil property mapping (Dvorakova et al., 2021; Zepp et al., 2021; Broeg et al., 2024; Tziolas et al., 2024). The advantage of this technique is to overcome the temporal coverage of soils with vegetation by selecting only bare soil reflectance pixels from the multispectral ∗Corresponding author. E-mail address: [email protected] (P. Karlshoefer). time stack (Diek et al., 2017; Rogge et al., 2018; Demattê et al., 2018; Roberts et al., 2019; Heiden et al., 2022). This way, the area for which soil properties can be mapped is increased. The quality of bare surface reflectance composites (SRC) depends on the selection of mostly undisturbed pixels representing bare soils and non-vegetated surfaces. Disturbances may arise from partial coverage of the soil by photosynthetically active and non-active vegetation (Rogge et al., 2018; Demattê et al., 2018) and by soil roughness and moisture (Dvorakova et al., 2023; Vaudour et al., 2021). The majority of studies use one or more spectral indices to disentangle such disturbances (Heiden et al., 2022) with the choice of indices and associated thresholds tailored to the local or regional conditions of the area of interest (Castaldi et al., 2023; Dvorakova et al., 2023; Broeg et al., https://doi.org/10.1016/j.geoderma.2025.117340 Received 5 February 2025; Received in revised form 8 May 2025; Accepted 11 May 2025 Geoderma 459 (2025) 117340 Available online 7 June 2025 0016-7061/© 2025 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). P. Karlshoefer et al. 2024). Studies aiming at larger extents, such as continental (Safanelli et al., 2020) to global scales (Rizzo et al., 2023) are rare due to the complexity of capturing the nuanced variations of soil properties across different regions (Gallo et al., 2018) and the computational complexity. Regionalization of data models has been shown to improve results in various global environmental studies, and it is a key component of the framework presented in this work. Notable examples include the mapping of forest canopy height (Ku et al., 2021), biodiversity (Coops et al., 2018), bare ground gain (Ying et al., 2017) and streamflow of river basins (Odusanya et al., 2022). To the best of the authors’ knowledge, this has not yet been tested for the European-wide or continent-scale SRC generation. The reason might be the lack of concepts that rigorously evaluate the performance of various spectral indices paired with a regionalized set of thresholds derived from local characteristics in the reflectance data. Another significant factor affecting SRC quality is cloud cover. Most temporal compositing approaches use the cloud masks produced during conversion from radiance values (Level 1) to reflectance values (Level 2) or additional software such as FMask (Baetens et al., 2019). Cloud masks, while effective in detecting highly affected pixels, must be conservative by design to avoid excessive data exclusion. While thick clouds can be detected with low uncertainties, the flagging of thin and semi-transparent clouds remains difficult (Skakun et al., 2022). Studies also reported spectral confusion between clouds and high reflectance built-up areas (Corbane et al., 2020) as well as in areas with low vegetation and high albedo characteristics, typical for Mediterranean soils (Ben-Dor, 1994). Kempeneers and Soille (2017) reduce the cloud contamination by selecting images based on quicklook RGB information and a minimum index value for the Sentinel-2 blue band. Other methods typically address this by taking the median of reflectance values along the time axis, assuming that cloudy observations are outliers (Simonetti et al., 2021). As an alternative to the median for the temporal averaging, the mean along the temporal axis of the reflectance values can be taken (Gasmi et al., 2021). A fundamental advantage is the greater asymptotic efficiency of the mean over the median in statistically consistent and normally distributed data. For instance, given 10 samples (samples being bare surface pixels), the relative efficiency of the median is only approximately 80% that of the mean (Serfling, 2011). The observed phenomenon is assumed to be relatively stable in time; however, the dynamics in the atmosphere (e.g. haze) and surface (e.g. surface moisture) lead to a high variance in the observed reflectance data and therefore a tail-heavy distribution. The mean is inherently more biased towards outliers and, therefore, it has to be combined with robust outlier detection to specifically remove pixels affected by residual clouds. In addition to its better convergence to the true mean of the distribution, the mean enables the derivation of statistical products centered around it, such as the standard deviation or the confidence interval, which allows to characterizes the spectral variability uncertainty of the SRC pixels. All of the issues described above influence the quality of SRC and thus affect the subsequent derivation of soil properties. To develop and test new concepts to improve the SRC at large scales such as continents, a robust and efficient evaluation concept is necessary. In many studies, the quality of different SRC versions are indirectly evaluated using the prediction accuracy of the derived soil properties such as described in Vaudour et al. (2021), Žížala et al. (2019) and Tziolas et al. (2020a). However, this strategy introduces additional uncertainty related to the machine learning models. For local studies, the resulting SRC spectra have been directly compared with laboratory spectra using the soil line and spectral dispersion using principle component analysis (Demattê et al., 2018). Heiden et al. (2022) use spectral distance measures between the SRC and laboratory spectra of the Copernicus Land Use-Land Cover Area Frame Survey (LUCAS) for Bavaria in Germany . The evaluation of SRC at large scale such as the European continent is still unsolved due to missing reference data that represent averaged spectral information across different seasons and years. However, Nocita et al. (2015) proved the statistical relation between spectral reflectance measurements of soils and their chemical properties. As a consequence, many of the soil parameter mapping approaches have used the spectral and chemical analysis collected in local laboratory and field surveys (Castaldi et al., 2018; Tziolas et al., 2020b; Mzid et al., 2022) or in the LUCAS survey (Ward et al., 2020; Castaldi et al., 2019; Safanelli et al., 2020) as a reference. Thus, for an Europeanwide evaluation of SRCs, the LUCAS survey seems to be a suitable reference that contains chemical and spectral measurements of the soil following the same protocol (Orgiazzi et al., 2017). The challenge is the different measurement method between the SRC and LUCAS spectra. In particular, the spectral measurements of the LUCAS survey are derived from dried and sieved soil samples at high radiometric resolution, making direct comparisons with space-borne soil reflectance under natural field conditions challenging. To allow comparison with Sentinel-2-based surface reflectance, the high-resolution LUCAS soil spectra are often resampled to match the multi-spectral resolution of Sentinel-2 (Castaldi et al., 2019). This resampling process smooths out the distinctive narrow spectral absorption features, effectively filtering out high-frequency information. However, the remaining broad spectral shape — like convexity and overall curvature can be used to compare SRC and LUCAS spectra. This study aims to introduce a systematic evaluation framework designed to quantify improvements in SRC generation methods, moving beyond traditional visual assessments. By optimizing key processing parameters (e.g. spectral indices and cloud-filtering) through an iterative process that maximizes concordance between SRC spectra and resampled LUCAS spectra at corresponding pixel locations, we aim to identify good processing settings with reasonable computational effort. We hypothesize that this framework can determine the choice of the spectral index and its corresponding thresholds (regional or global) and quantify the benefit of further subroutines in the SRC generation. This framework is applied to (1) evaluate SRC generation enhancements introduced in this work, while building upon developments from the Soil Composite Mapping Processor (SCMaP) as described in previous studies (Rogge et al., 2018; Heiden et al., 2022) and (2) it is used to compare other SRC generation concepts that have been published. The developed framework can provide quantitative and objective means to define the best SRC processing parameter for a specific area. 2. Material and methods This section outlines the methodological foundation of the evaluation framework designed to quantify the influence of processing parameters on the quality of bare surface reflectance composites (SRCs). A schematic overview of the evaluation framework is provided by Fig. 1. The main idea is to search for the optimal set of processing parameters (blue box in Fig. 1) identified by minimizing the angular distance between the derived SRC spectra and the reference spectra through an iterative refinement of those parameters. The processing parameters are generic and can be defined as proposed in this study or taken from the literature such as provided by Diek et al. (2017), Demattê et al. (2018) and Broeg et al. (2024). In this study, the proposed parameters include the selection of spectral indices and additional bare surface detection routines suitable to minimize the influence of green and dry vegetation and to filter out cloudy and hazy pixels. A common temporal stack of Sentinel-2 L2A data (Section 2.1.1) is used ensuring consistency across all optimization runs and enabling a fair comparison of the performance of the spectral indices. The reference spectra are obtained from the European LUCAS 2015 survey (Section 2.1.3). The SRC generation itself remains intentionally generic within this framework, with its key components determined during the optimization process (Section 2.2). Fig. 2 illustrates the generation of a SRC image, where the parameters found in the optimization loop are used Geoderma 459 (2025) 117340 2 P. Karlshoefer et al. Fig. 1. The evaluation framework: The spectral reference data is taken from the LUCAS survey. The temporal stack of the Sentinel-2 L2A data is identical for any SRC, that is being generated in this study. Processing parameters are iteratively optimized via gradient descent until the spectral angular distance between reference and SRC spectra no longer decreases, providing an estimate of the performance of the chosen parameters. Fig. 2. SRC generation enhancements: Pixel-based surface thresholds are derived using HISET and a refined landcover data set, that is based on ESA Worldcover. The resulting surface thresholds are then used in combination with the entire S2 temporal stack and the processing parameters found prior (blue), to generate the SRC. (blue). Instead of limiting the evaluation to LUCAS points, this approach enables the generation of a spatially continuous SRC image. A key component in this process is the HISET algorithm (Section 2.2.1), which is employed to construct surface thresholds essential for an areawide SRC derivation. This procedure relies on cropland and grassland masks, initially based on the ESA WorldCover dataset (Section 2.1.2) and subsequently refined to align with the specific requirements of this study (Section 2.2.2). This approach contrasts with existing studies, where the indices and thresholds are predefined. The proposed method offers two main advantages. First, it allows for parameter customization tailored to the region of interest, or, at minimum, to areas well represented by the reference data. Second, it facilitates the systematic evaluation of successive adaptations in the SRC generation process, eliminating the labor-intensive and error-prone task of visually comparing images to assess improvements in the final product (Section 2.3). In this work, the term ‘‘bare surface’’ encompasses non-vegetated terrestrial surfaces including soils, barren rock, and sandy surfaces. Artificial and urban surfaces, snow and water bodies are excluded. 2.1. Data preparation 2.1.1. Sentinel-2 The foundation of the SRC is a temporal stack of Sentinel-2 multispectral imagery. This archive is well-suited for the analysis due to the Sentinel-2 constellation’s ability to consistently capture data with a five-day revisit period (European Space Agency (ESA), 2015), providing a dense temporal resolution critical for capturing the rare observation opportunities of non-vegetated, bare surfaces. The selection process for the input data used to generate all tested SRCs in this paper is filtered based on several criteria: •Metadata Filtering: Images with a maximum cloud coverage of 80 percent are considered. This value is strategically chosen to maximize the density of the temporal stack, acknowledging that high cloud cover could still yield valuable data. Furthermore, the sun has to be elevated more than 20 degrees between the senor and the horizontal plane. •Temporal Range: The study period spans from January 2018 to December 2022. Geoderma 459 (2025) 117340 3 P. Karlshoefer et al. Table 1 Sentinel 2 bands used in this work. (Copernicus, 2020). Band Central wavelength (S2A) Central wavelength (S2B) B02 (blue) 492 nm 492 nm B03 (green) 560 nm 559 nm B04 (red) 664 nm 665 nm B05 (red edge 1) 704 nm 704 nm B06 (red edge 2) 741 nm 739 nm B07 (red edge 3) 783 nm 780 nm B08 (NIR) 833 nm 833 nm B8A (red edge 4) 865 nm 864 nm B11 (SWIR 1) 1614 nm 1610 nm B12 (SWIR 2) 2202 nm 2186 nm •Spectral Band Selection: Only bands with a spatial resolution at ground level of 20 m or better are included (see Table 1). The spatial organization of the input images adheres to the Military Grid Reference System (MGRS) (Meyers, 2015), which facilitates systematic processing and analysis. The data were pre-processing to Level-2A (L2A) reflectances using the MAJA Atmospheric Correction processor (Hagolle et al., 2010). This step adjusts for atmospheric effects and cloud identification and ensures radiometric consistency across the image stack. The used scenes are corrected for adjacency and topographic effects (denoted by the suffix FRE) (Hagolle et al., 2018, 2021). For this study, the analysis is confined to pixels classified as ‘0’ by the geophysical mask (*_MG2_R2), effectively excluding areas represented as open water, clouds, snow, and those obscured by topographic shadows or unsuitable sun elevation angles. This selection criterion ensures that the focus remains on the surfaces that are most relevant to the SRC. In total, a 1145 MGRS tiles and a total of 445806 individual Sentinel-2 scenes (each comprised of 10 bands) contributed to the computation. 2.1.2. Landcover mask The integration of land cover data, specifically the ESA WorldCover map (Zanaga et al., 2022), serves as an input for the estimation of surface thresholds and is therefore only part of HISET (Section 2.2.1). Therefore, it functions as an indirect input for the final composite rather than as a direct masking layer. Designed for the year 2020, ESA WorldCover is centrally positioned within the temporal stack of Sentinel-2 data. With a spatial resolution of 10 m and precise classification of critical land-cover types, such as grasslands and croplands, it is - with refinement - a central input to the threshold derivation, which is detailed in Section 2.2.1. 2.1.3. LUCAS spectral library The LUCAS topsoil spectral database (2015) is used as reference for the general shape of soil spectra at their respective sample sites. Important for this work are the more than 22,000 topsoil spectra from European soil samples. To increase the likelihood of bare soil visibility for optical remote sensing, our analysis focuses specifically on the 8823 data points located within agricultural fields (labeled as ‘‘Cropland’’). Another 149 points have been removed from the set due to low spectral variability in time (according to Eq. (6)) in their local vicinity, in an effort to exclude sites which have turned into inactive cropland since 2025. On such fields, their might be permanent vegetation with no option to observe the bare surface. To allow for comparison between the LUCAS spectra, 𝑠𝑙𝑢𝑐𝑎𝑠 and the Sentinel-2 spectra, the LUCAS spectra are projected into Sentinel-2 data space using the spectral response function 𝑠𝑟𝐴 and 𝑠𝑟𝐵 for Sentinel-2 A and B, respectively: 𝑞 =𝑠𝑟𝐴+𝑠𝑟𝐵 2𝑠𝑙𝑢𝑐𝑎𝑠 (1) Using Eq. (1), the high-resolution LUCAS spectra 𝑠𝑙𝑢𝑐𝑎𝑠 are resampled to match the spectral resolution of Sentinel-2 by the average of the spectral response functions of Sentinel-2 A and B (Castaldi et al., 2019; Okujeni et al., 2024). Due to the radiometric resampling to ten bands, only the general shape of the LUCAS spectra is preserved, while detailed absorption features are eliminated. This, combined with an angular distance metric that is invariant to the scaling of the spectral vector, ensures that the primary characteristic compared between composite spectra and 𝑞 is their overall shape of the spectra. 2.2. Procedure to generate bare surface reflectance composites Starting from the foundational element, a single observation, its multi-spectral cube is 𝑂∈R𝑢×𝑣×𝑏 for an image of spatial dimension 𝑢×𝑣 and 𝑏 spectral bands. Given that observations are geometrically co-registered and mapped onto the Sentinel-2 tile grid, thus, each image sharing uniform dimensions, the tensor can be directly extended by a temporal component 𝑆= (𝑂0, 𝑂1,…, 𝑂𝑛) over 𝑛 observations within a tile. The synthesis of a SRC 𝐶 from 𝑆 is executed through two pivotal functions. Firstly, the classifier 𝑓𝑏𝑎𝑟𝑒(𝑆)→{0|1}𝑢×𝑣×𝑏×𝑛, discriminates each element within 𝑆 whether it is representative of undisturbed bare surface (1) or not (0). Using the mask will result in the set 𝑆𝑏𝑎𝑟𝑒 ∶= {𝑠∈𝑆|𝑓𝑏𝑎𝑟𝑒(𝑠)>0}.(2) In our case, 𝑓𝑏𝑎𝑟𝑒 is primarily based on a spectral index, 𝑧, thus: 𝑓𝑏𝑎𝑟𝑒 ∶= {1,if 𝑡0,(𝑖,𝑗)< 𝑧(𝑠𝑖,𝑗,∶,𝑚)< 𝑡1,(𝑖,𝑗), 𝑖 ∈ {1..𝑢}, 𝑗 ∈ {1..𝑣}, 𝑚 ∈ {1..𝑛} 0,otherwise (3) Eq. (3) masks all pixels indexed by the coordinates 𝑖 and 𝑗 that are not in between their local set of thresholds 𝑡0,(𝑖,𝑗) and 𝑡1,(𝑖,𝑗) for index 𝑧. It is extended by an outlier detection and removal which is detailed in the next chapter. Secondly, the aggregation 𝐶(𝑆)→R𝑢×𝑣×𝑏 condenses the temporal dimension to a single multi-spectral image using the definition of 𝑆𝑏𝑎𝑟𝑒 (Eq. (2)) 𝐶(𝑆𝑏𝑎𝑟𝑒) = 1 |𝑆bare|∑ 𝑟∈𝑆bare 𝑟. (4) In this work, we deliberately chose the mean over the median for its greater asymptotic relative efficiency; thus, its tendency to converge faster to the true mean of the distribution (Serfling, 2011). The delineation of bare surfaces from vegetated areas in our analysis hinges on the precise determination of thresholds applied to a carefully chosen spectral index, crucial for identifying natural bare surfaces. These thresholds, 𝑡0 and 𝑡1, are defined within the bounds of a spectral vegetation index 𝑧 calculated only from a reflectance vector 𝑠, which is composed from the spectral band values 𝑠= (𝑠𝐵2, 𝑠𝐵3, 𝑠𝐵4, 𝑠𝐵5, 𝑠𝐵6, 𝑠𝐵7, 𝑠𝐵8, 𝑠𝐵8𝐴, 𝑠𝐵11, 𝑠𝐵12)𝑇: 𝑡0< 𝑧(𝑠)< 𝑡1(5) The spectral index 𝑧∶R10 →R serves as the basis for this classification, with higher index values indicating an increased presence of vegetation. The selection of 𝑡1 is particularly critical as it delimits the transition from bare to vegetated surfaces. We hypothesize that these thresholds require careful calibration with respect to the local characteristics of soil and crops. Within the scope of this work, spectral indices commonly found in the literature are investigated (Table 2). Single and multi-staged indices are tested, like NDVI/NBR2, NDVI/NBR, BCC/NDVI/NBR2, VNSIR/NDVI/NBR2, where the indices are applied to a spectrum consecutively, each with its own threshold. VNSIR and BCC are only tested as proposed by the authors staged with filtering by NDVI and NBR2, first. The BSI has been inverted so lower values correspond to bare Geoderma 459 (2025) 117340 4 P. Karlshoefer et al. Table 2 Spectral indices tested in this work for their suitability to generate bare surface composites. Abbreviations of the indices and their original source (column 1), their respective formulas in terms of Sentinel-2 bands (column 2) and their full names and usage (column3) are shown. Index Formula Usage NDVI (Rouse et al., 1974)(𝐵8−𝐵4 𝐵8+𝐵4)Normalized Differences Vegetation Index used by Urbina-Salazar et al. (2021) NBR2 (Deventer et al., 1997)(𝐵11−𝐵12 𝐵11+𝐵12 )Normalized Burn Ratio 2 used by Vaudour et al. (2021) NBR (García and Caselles, 1991)(𝐵8−𝐵12 𝐵8+𝐵12 )Normalized Burn Ratio BSI (Rikimaru et al., 2002) −(𝐵12+𝐵4−𝐵8𝐴−𝐵2 𝐵12+𝐵4+𝐵8𝐴+𝐵2)Bare soil index (inverted) used by Diek et al. (2017) NDVI+NBR (𝐵8−𝐵4 𝐵8+𝐵4+𝐵8−𝐵12 𝐵8+𝐵12 )Linear combination of NDVI and NBR used by Heiden et al. (2022) NDVI/NBR2 ((𝐵8−𝐵4 𝐵8+𝐵4),(𝐵11−𝐵12 𝐵11+𝐵12 )) Staged NDVI and NBR2 with 2 separate sets of thresholds used in Dvorakova et al. (2021) NDVI/NBR ((𝐵8−𝐵4 𝐵8+𝐵4),(𝐵8−𝐵12 𝐵8+𝐵12 )) Staged NDVI and NBR with 2 separate sets of thresholds VNSIR (Demattê et al., 2020) 1−((2 ∗ 𝐵4 − 𝐵3 − 𝐵2) +3 (𝐵12 − 𝐵8))∕10 000 Visible to Shortwave Infrared Tendency Index staged with NDVI/NBR2 3 sets of thresholds, used by Demattê et al. (2020) BCC (Gillespie et al., 1987)(𝐵2 𝐵4+𝐵3+𝐵2)Blue Chromatic Coordinate staged with NDVI/NBR2 3 sets of thresholds, used by Broeg et al. (2024) surfaces (similar to NDVI, NBR2 etc.). Linear combinations are also possible, where two indices are combined and hence can be treated with a single threshold (like NDVI+NBR). 2.2.1. Surface thresholds (HISET) This study employs two distinct strategies for deriving thresholds. The first strategy involves an optimization process, as detailed in Section 2.3, which computes a set of thresholds at each reference point (LUCAS sample site). The second strategy derives thresholds for SRC generation using an extended version of the Histogram Separation Threshold (HISET) algorithm, introduced by Heiden et al. (2022). These thresholds, termed surface thresholds, are applicable to entire areas rather than singular points. Ultimately, the surface thresholds are utilized for the SRC generation, as the thresholds computed at the LUCAS points are specific to the sparse grid of reference points and are only used to determine the best spectral index. The general principles of HISET, along with the modifications required to ensure its applicability across diverse conditions, are elaborated in this section. The surface thresholds adhere to the Sentinel-2 tile-grid. It is chosen to be the spatial skeleton, where thresholds are determined for each tile in the grid. The granularity of roughly 100 km2 ensures a balance in scale that avoids both the over-generalization of large regions and the over-fitting risks of excessively small ones. To refine spatial continuity and minimize edge effects within the resulting data, these thresholds were smoothed employing a radial basis function centered on the centroid of each tile, see result in Fig. 6. This allows for a seamless integration of tiles into the final composite. Heiden et al. (2022) details the HISET algorithm in the general case. It is basically designed to account for the spectral similarity for the non-photosynthetically active vegetation (NPV) and bare soil (Daughtry and Hunt, 2008; Dennison et al., 2019, 2023). Since Sentinel-2 only has two spectral bands in the SWIR, the usage of spectral absorption features to detect NPV is not an option. . Therefore, HISET measures the spectral overlap between bare soils and NPV and defines a threshold that minimizes the spectral overlap. HISET encompasses four primary steps: 1. Spectral Index Calculation: Compute the spectral vegetation index for each image and pixel in the temporal Sentinel-2 series 2. Minimum Index Selection: Selecting the minimum spectral index value across the temporal sequence per pixel ensures representation in its least vegetated state, typically indicating bare surfaces and non-photosynthetic vegetation (NPV) for croplands and grasslands, respectively. 3. Histogram Aggregation: All pixels that belong to cropland or grassland according to the land-cover map (Zanaga et al., 2022) are compiled into two separate histograms. Each histogram is normalized to an area of 1, effectively transforming it to a discrete probably density function. 4. Threshold Identification: The threshold is chosen to be the center of the histogram bin that minimizes the larger of both areas, that is intersected, at either side (see Fig. 4). In addition to the value for the threshold, the area that has been minimized in step 4, gives an indication to the degree of separation of both histograms. Due to normalization, this value is between 0 for two disjointed histograms and 0.5 when one is entirely contained within the other. Multiplied by 100, this value is the quality score of the threshold. It serves as an indicator for the effectiveness of the surface threshold. 2.2.2. Enhancements to the surface threshold derivation The initial methodology (HISET) for threshold derivation, as outlined in the previous Section 2.2.1, was further refined to address the challenges posed by the complexity of bare surfaces in Europe. A notable concern was the presence of mixed pixels, that embody areas where agricultural land intersects with non-agricultural elements, such as roads or cart tracks, within the same pixel. These elements interfere with the spectrum and consequently the spectral index, introducing variance to the histograms. This effect leads to blurred histograms and consequently elevated scores of the threshold, that indicate a degraded or inaccurate thresholds. To reduce the effect of mixed pixels, we introduce a measure of spectral variability: 𝑣∶= 1 𝑛− 1 𝑛−1 ∑ 𝑖=0 ||||| 𝑧(𝐻𝑖+1) − 𝑧(𝐻𝑖) 𝑑𝑖+1 −𝑑𝑖||||| (6) Geoderma 459 (2025) 117340 5 P. Karlshoefer et al. Fig. 3. The high-resolution google earth image of an agricultural area in southern Germany (a) shows the overestimation of cropland pixels by the ESA Worldcover landcover map (b). By excluding pixels with low spectral variability in the vegetation index (darker areas in the (c)), the land-cover map is refined (d). Some fields without active management are excluded by this process as well as many small roads. Here, 𝐻= (𝑠0, 𝑠1,…, 𝑠𝑛) is a sequence of averaged spectra over two months in time with 𝑑 indicating the respective central dates. The initial temporal averaging of two months of the reflectance data eliminates signals at a higher frequency, that are not related to the crop lifetime cycle. The spectral variability 𝑣 estimates the frequency and intensity of changes to agricultural surfaces per pixel and therefore the likelihood that a bare state is exposed at some point in time. It is comparable only at local scale, as climatic and seasonal effect or land management policies do impact it. Greater spectral variability translates to actively manged croplands, that are less likely to intersect a stable surface. Thus, the land-cover map for cropland is refined by this additional criterion. Depending on the overall abundance of cropland pixels, a percentile (up to the 20th-percentile) of less active pixels is removed or the pixel’s contribution to the histogram is weighted according to its variability. Beyond considering spectral variability, pixels that contribute to expansive, contiguous areas are prioritized. In contrast, isolated pixels and small patches — regardless of classification — are excluded from contributing to the histograms. They are found through a parallel flood-fill algorithm (van der Walt et al., 2014). The resulting land-cover map can be seen in Fig. 3. It gives an example of the identification of both mixed pixels at field edges and fields with comparatively low variability in the vegetation index (red circle). Both are excluded, so the surface threshold estimation is based on more pixels, which actually represent the bare surface. While the HISET methodology shows robust performance across the majority of Europe, challenges arise in tiles with limited cropland or grassland coverage, leading to imbalanced or erratic density functions and ultimately, poor surface thresholds. Such cases are notably prevalent in coastal regions and islands where a large portion of the tile is water, as well as northern regions with inherently sparse cropland. Various tests indicated, that both land-cover classes should be represented by at least two percent of all pixels in the tile, (at least 2𝑒5), thereby a tile is classified as fit. In instances where this criterion is not met, alternative approaches are employed to estimate thresholds: •Proximity Interpolation: For tiles adjacent to one or more ‘fit’ tiles, thresholds are linearly interpolated, utilizing the centroid of each tile as a basis for distance calculation. This method is particularly relevant for coastal regions. •Histogram Aggregation: In the absence of nearby ‘fit’ tiles, as often found with islands, pixel data from multiple adjacent tiles are combined into a singular histogram. The resulting threshold is used for all participating tiles. •Bio-geographic Aggregation: In the case of northern Scandinavia, the system of multiple climatic zones led to an increased variance in the combined histograms. Here, pixels are aggregated by bio-geographic zone (EEA, 2016). Thresholds are then calculated for each zone and applied proportionately across the tiles based on the area of each bio-geographic zone. These tailored approaches ensure the methodology’s effectiveness, even in regions where traditional threshold derivation faced significant challenges. Geoderma 459 (2025) 117340 6 P. Karlshoefer et al. 2.2.3. Detection of residual clouds and haze Eq. (3) can be extended by a set of cloud filters, integrated to specifically address the simpler challenge of separating cloudy and hazy pixels from those representing a bare surface. This excludes pixels that are affected by haze or clouds. Given the high brightness of clouds in contrast to most soils, this is important. The first filter exploits the distinct spectral characteristics between soil and cloud pixels particularly in the SWIR. While soil spectra have generally higher values in the short-waved infrared (SWIR) compared to the visible to near-infrared (VNIR) (true for more than 99% of LUCAS soil spectra), the opposite is true for clouds. Therefore the condition can be formulated as 𝑠𝐵12 −𝑠𝐵8 ! >0, 𝑠 ∈𝑆𝑏𝑎𝑟𝑒,(7) where the subscript denotes the selected Sentinel-2 band. The second filter within 𝑓𝑏𝑎𝑟𝑒 focuses on the removal of pixels affected by haze or thin clouds. The influence is very prominent in the blue band (Band B02), which is consequently used for detection (B02 is also extensively used in ESA (2021)). Given the challenge that the reflectance in the B02 band of the underlying soil is unknown, this filter detects outliers along the temporal axis in relation to the median of the dataset in the blue band. The normalized mean absolute deviation (NMAD) for in the blue band 𝑆𝑠𝑜𝑖𝑙,𝐵2= {𝑠𝐵2,∀𝑠∈𝑆𝑏𝑎𝑟𝑒} is given by:  𝑆𝑠𝑜𝑖𝑙,𝐵2=𝑚𝑒𝑑𝑖𝑎𝑛𝑡(𝑆𝑠𝑜𝑖𝑙,𝐵2) Absolute Deviations = {|𝑥− 𝑆𝑠𝑜𝑖𝑙,𝐵2|∶𝑥∈𝑆𝑠𝑜𝑖𝑙,𝐵2} NMAD =𝑘⋅𝑚𝑒𝑑𝑖𝑎𝑛(Absolute Deviations)(8) where 𝑘= 1.4826, a constant scaling factor for data adhering to a normal distribution and 𝑚𝑒𝑑𝑖𝑎𝑛𝑡 denotes the median along the temporal axis. Eventually, reflectance values are considered only if they fulfill the condition: 𝑠𝐵2− 𝑆𝑠𝑜𝑖𝑙,𝐵2 ! < 𝜎 ⋅NMAD, 𝑠 ∈𝑆𝑏𝑎𝑟𝑒 (9) Typically, 𝜎 is set to 3. The conditions (7) and (9) serve to adapt the scope of cloud exclusion of the MAJA processor to a stricter set of rules. Consequently, only pixels that satisfy these criteria are selected for inclusion in the SRC, contributing to a more accurate reflection of the spectral signature of the actual surface. 2.3. SRC evaluation Determining the optimal spectral vegetation index, 𝑧(𝑠), is critical for the 𝑓𝑏𝑎𝑟𝑒 classifier, ensuring it accurately distinguishes bare soil and surface pixels from other pixels using only the Sentinel-2 spectral bands’ reflectance vector, 𝑠. The index 𝑧 is parameterized by a set of thresholds, so it is able to adapt to the diverse climatic and surface conditions. By finding the optimal set of thresholds for various choices of 𝑧 at each evaluation pixel, one can draw conclusions about the suitability of the chosen spectral index. Ultimately, this performance estimation can serve as a baseline for testing the whole SRC generation and therefore can quantify improvements made to the processor. To evaluate 𝑧, we measure the congruence between the SRCs generated using a particular 𝑧 and its associated thresholds, and data measured on the surface. These reference spectra 𝑞 are derived from the LUCAS topsoil database (2015) (Orgiazzi et al., 2018), which was resampled according to Section 2.1.3. To quantify the difference between a single composite spectrum 𝑐 and 𝑞, literature proposes the spectral angle, which is derived from the cosine similarity: 𝛼(𝑐, 𝑠𝑠2) = 𝑐𝑜𝑠−1 (𝑐 ⋅𝑞 ‖𝑐‖‖𝑞‖).(10) The measure is invariant to global changes in brightness, e.g. the length of the vector in its vector-space (van der Meer, 2006), which is important since the reference spectra are scaled differently to the remote sensing spectra. The establishment of a reference dataset, alongside a mechanism for its comparison against Sentinel-2 data sets the stage for an optimization problem. Our objective is to find, for a choice of 𝑧 at each LUCAS point an ideal set of thresholds, that ensure maximum concordance with the shape of the LUCAS reference. For the 𝑗th LUCAS point, a set of thresholds that are used to compute the composite 𝐶(𝑧, 𝑔) are required. The dimension of this set depend on the form of 𝑧 (for example 𝑧 can be a set of multiple spectral indices). The set of (𝑡0, 𝑡1) pairs is called 𝜏 to avoid confusing subscripts. The task then becomes identifying the optimal 𝜏𝑗= (𝜏(𝑗,(1)), 𝜏(𝑗,(2))..)𝑇 for each point that minimizes the angular distance between the SRC and the ground truth: 𝜙𝑗(𝑧) = 𝑚𝑖𝑛 𝑗(𝛼(𝐶(𝑧, 𝜏𝑗), 𝑞𝑗)).(11) 𝐶(𝑧, 𝜏) and in result Eq. (11) are not smooth nor differentiable. However, for the optimization, this is necessary. Therefore, the set notation of Eq. (2) is replaced by a soft exclusion using the sigmoid 𝜎(𝑥) = 1 1+𝑒−𝛽𝑥 . The soft exclusion is based on the distance between the current guess of the thresholds of 𝜏 to the actual value of the spectral index for all 𝑠∈𝑆.  𝐶(𝑧, 𝜏) = ∑𝑖𝜎(𝑧(𝑋𝑖) − 𝜏)𝑠 ∑𝑖𝜎(𝑧(𝑋𝑖) − 𝜏) ,∀𝑠∈𝑆(12) The parameter 𝛽 is the sharpness of the soft exclusion. It is set to 15 in this case. The performance 𝑝 of 𝑧 using its optimal set of thresholds 𝜏𝑜𝑝𝑡 is the average angular distance over all 𝑚 LUCAS points 𝑝(𝑧) = 1 𝑚 𝑚 ∑ 𝑗=1 𝜙𝑗(𝑧).(13) This is the central measure to evaluate the SRC. The 𝜏𝑜𝑝𝑡 𝑗 that minimizes Eq. (11) is found iteratively at each LUCAS point, using gradient descent: 𝜏𝑛+1 𝑗=𝜏𝑛 𝑗−𝛾∇𝛼(𝐶(𝑧, 𝜏𝑛 𝑗), 𝑞𝑗),(14) with ∇ = (𝜕 𝜕𝜏𝑗,(1) ,𝜕 𝜕𝜏𝑗,(2) ,…) and a learning rate 𝛾 > 0. The partial derivatives are approximated numerically using finite differences. In addition to Eq. (11), there is an additional condition that the number of spectra composing the final composite must be at least five for each pixel. This requirement reduces the likelihood of overfitting the thresholds 𝜏𝑜𝑝𝑡. For normalized spectral indices that are used in this work, 𝛾 is set to 𝛾=𝑡𝑚𝑎𝑥 −𝑡𝑚𝑖𝑛 200 , 𝑧 ←←→ [𝑡𝑚𝑖𝑛, 𝑡𝑚𝑎𝑥].(15) Assuming that Eq. (11) is locally convex, the search described in Eq. (14) is terminated upon reaching a local minimum. By restarting the algorithm with different initial guesses that converge to the similar 𝜏𝑜𝑝𝑡, we assume that Eq. (11) is indeed convex in the vicinity of its minimum. 2.4. Software and tools The computations related to the optimization were performed on the high-performance data analytics platform terrabyte hosted at the Leibniz Supercomputing Center. The source code is mostly written in Python and leverages the power of vectorized computations implemented in numpy and scipy for calculations, as well as gdal to handle IO. A spectral library of all Sentinel-2 observations in time at all LUCAS points allows to generate single pixel SRCs quickly, crucial for quick convergence during the gradient descent method. The final SRC at continental scale was computed using the Soil Composite Mapping Processor (SCMaP) from DLR. Geoderma 459 (2025) 117340 7 P. Karlshoefer et al. Table 3 Performance 𝑝 according to Eq. (13), computed from a bare surface discriminator 𝑓𝑏𝑎𝑟𝑒 that is based on the spectral index 𝑧. Lower values of 𝑝 indicate greater spectral similarity. For readability, 𝑝 was normalized to the performance of NDVI+NBR. Spectral index 𝑧[𝑡𝑚𝑖𝑛, 𝑡𝑚𝑎𝑥]Initial guess 𝜏0𝑝(𝑧) normal. NDVI (𝐵8−𝐵4 𝐵8+𝐵4)[−1, +1] (−0.25,0.25)𝑇1.55 NBR2 (𝐵11−𝐵12 𝐵11+𝐵12 )[−1, +1] (−0.3,0.1)𝑇1.26 NDVI/NBR2 [−1, +1] ((−0.25,0.25),(−0.3,0.1))𝑇1.11 BSI −(𝐵12+𝐵4−𝐵8𝐴−𝐵2 𝐵12+𝐵4+𝐵8𝐴+𝐵2)[−1, +1] (−0.1,0.1)𝑇1.09 NBR (𝐵8−𝐵12 𝐵8+𝐵12 )[−1, +1] (−0.3,0.1)𝑇1.05 NDVI/NBR2/VNSIR [−1,+1]2, R((−0.25,0.25),(−0.3,0.1),(0,1))𝑇1.04 NDVI+NBR (𝐵8−𝐵4 𝐵8+𝐵4+𝐵8−𝐵12 𝐵8+𝐵12 )[−2, +2] (−0.6,0.25)𝑇1 NDVI/NBR (𝐵8−𝐵4 𝐵8+𝐵4,𝐵8−𝐵12 𝐵8+𝐵12 )[−1,+1]2((−0.6,0.25),(−0.5,0.1))𝑇0.98 NDVI/NBR2/BCC [−1,+1]2, R((−0.25,0.25),(−0.3,0.1),(0,0.3))𝑇0.96 3. Results 3.1. Selection of the spectral index First, common choices for a spectral vegetation index 𝑧 from publications are investigated (Heiden et al., 2022; Broeg et al., 2024; Demattê et al., 2020; Diek et al., 2017). Using the method outlined in Section 2.3 to obtain the best set of thresholds for each index at each LUCAS point, the Table 3 shows the performance 𝑝(𝑧) of these indices in calculating an accurate SRC spectrum at these points. The performance is measured according to Eq. (13), thus, lower values signify greater spectral similarity. For readability, the last column is normalized to the value of NDVI+NBR. Table 3 illustrates that NDVI (1.55 times worse than NDVI+NBR) and NBR2 (1.26 times worse than NDVI+NBR) on their own are poor choices for the classification of bare surfaces. Combining them into a two-stage spectral index with separate, independent thresholds improves their performance to 1.11 times performance of NDVI+NBR. Demattê et al. (2020) proposed to add a third stage to the classifier called VNSIR: 1 − ((2 ∗ 𝐵4 − 𝐵3 − 𝐵2)+ 3 (𝐵12 − 𝐵8)) ∕10000 (see Table 2). This improves the combined index substantially to just 1.04 time less accurate compared to the NDVI+NBR index. The BSI (Rikimaru et al., 2002; Diek et al., 2017) manages to produce SRC spectra about 1.09 times less accurate than the NDVI+NBR index and the NBR index performs better, at just 1.05 times worse compared to NDVI+NBR. Eventually, by adding the NBR and the NDVI (NDVI+NBR), such as proposed by Heiden et al. (2022) , we obtain the most effective, single-stage index combination of this list in selecting bare surface pixels at the LUCAS sites. The two-staged NDVI/NBR index with distinct thresholds for both components or the three staged NDVI/NBR2/BCC (Broeg et al., 2024) improve performance (factors 0.98 and 0.96, respectively). However, further analysis will be based on the NDVI+NBR index. It has the highest performance of all indices that have only a one-dimensional set of thresholds (see column 2), which substantially reduces complexity and robustness in the actual SRC generation implementation. 3.2. Validating enhancements to the soil detection and threshold derivation The thresholds 𝜏𝑜𝑝𝑡 computed in the previous section reflect for the NDVI+NBR index the optimal threshold set at each LUCAS site. These thresholds are tuned to individual pixels and potentially biased to the selection of the LUCAS points; therefore, they cannot be applied at broader scales, such as that of Sentinel-2 tiles. Nevertheless, we can use the 𝜏𝑜𝑝𝑡 to compare them to the performance 𝑝 of the SRC generation. Incorporating successive refinements into 𝑓𝑠𝑜𝑖𝑙, including additional cloud filter and land cover adjustments for threshold derivation, has quantifiably enhanced the accuracy of SRCs.Table 4 shows the Table 4 Gradual improvement of the performance of the SRCs built by surface thresholds. Values are normalized to the limit of the index. Lower values correspond to greater spectral similarity (via lower angular distance) and in result, a better SRC. Optimizations for SRC creation Performance 𝑝 of surface thresholds Basic threshold database (HISET) Section 2.2.1 1.16 + additional cloud filters for SRC creation Section 2.2.3 1.13 + land-cover optimization Section 2.2.2 1.07 Theoretical limit of the index NDVI+NBR 1 reduction in the angular distance from 116% to 107% of the theoretical limit achievable with the NDVI+NBR index at pixel level. Fig. 4 illustrates the effect of the optimization of the land cover map to estimate the thresholds. Using the vegetation variability defined in Eq. (6), areas with little change are excluded from HISET. The variability mask reveals that most often the removed areas are mixed pixels covering unmanaged fields, which consequently have less chance in exposing the bare surface. Removing these pixels (magenta region) from the threshold estimation allows for a cleaner separation of the two histograms (smaller score) and in result, stricter thresholds (solid line compared to dashed line). Table 4 confirms that throughout Europe the addition of cloud filters and a sharper determination of thresholds result in SRCs closer to the reference data. 3.3. Performance of regional thresholds The difference for the computed SRC at the LUCAS points between a set of regional thresholds versus a common, global threshold is illustrated in Fig. 5. This figure plots the performance of all existing composite spectra using the specified threshold versus the percentage of valid pixels in the bare surface composites. An invalid pixel means that less than five pixels in the temporal stack are labeled as bare surface, which is deemed insufficient for the composition (Dvorakova et al., 2023). Ideally, a high performance (low angular distance between the LUCAS and the SRC spectra), close to the LUCAS reference data, and a high percentage of valid pixels, are desired (upper left corner). As shown in Fig. 5, the HISET surface thresholds (+) manage to compute composite spectra for over >82% of LUCAS points, while maintaining a high performance (an average angular distance of 0.058). In contrast, all choices for a global threshold (x) either result in sparse SRCs (thresholds are excessively strict and exclude too many Geoderma 459 (2025) 117340 8 P. Karlshoefer et al. Fig. 4. Histograms of the spectral vegetation index NDVI+NBR in the Sentinel-2 tiles 30TYQ and 34UED for grassland (green) and croplands (orange) pixels (both normalized to an area of 1). The magenta region marks the pixels that have been removed by a refined landcover mask using the vegetation variability 𝑣, resulting in a sharper separation. Fig. 5. Performance 𝑝 of composite spectra using the regionalized surface HISET (marked with a plus (+)) versus a static global threshold (marked with an (x)) for the NDVI+NBR index. The star (*) indicates the theoretical upper limit, constrained by the condition to maintain at least 5 bare surface observations per SRC pixel. The fraction on the y-axis cannot exceed 1 (indicated by the horizontal dashed line). The dashed line connecting the (x) marker indicates that any choice for a global 𝑡1 is further away from (*) than (+). observations) or exhibit a low performance (larger values) to the reference, thus being less accurate. The SRC using the pixel-based (optimal) thresholds, that are computed at each LUCAS point, are located at the star, representing the theoretical performance limit of the NDVI+NBR index. The star is closer to the HISET surface thresholds than any choice for a global threshold, indicating that HISET surface thresholds is always a better compromise than any globally fixed threshold. In practice, comparing the performance of the regional HISET surface thresholds for the NDVI+NBR index with global thresholds found in the literature reveals a more significant improvement in quality than initially suggested by Fig. 5. Table 5 illustrates, that soil spectra produced by our surface thresholds are substantially closer to the reference data, compared to the method proposed by Demattê et al. (2020) in row three. This demonstrates that even with a well-chosen combination of indices, such as the case of the approach used in row three (the index performs only slightly worse than the NDVI+NBR, see Table 3), a global threshold can significantly impair results. Broeg et al. (2024) (fourth line of Table 5) use a fixed, global threshold initially that is refined on a pixel-level in a second stage. While this work is designed for Germany, it was applied in this study at the European scale. Here, the performance is by a factor of 1.14 times lower for the European continent. 3.4. Threshold distribution for NDVI+NBR Fig. 7 presents the resulting SRC from Sentinel-2 data between 2018 to 2022 displayed in true color (blue: B02, green: B03, red: B04). Pixels shown in white indicate areas with no data, where soil was not sufficiently often exposed during the observation period. The interpolation process to display this large dataset may give the impression of a dense map at larger scales in some places like central Europe. However, upon zooming in, as highlighted in the window on the right, it becomes evident that the map is fragmented. This fragmentation is expected, as many surfaces, such as permanently vegetated areas or urban developments, do not expose bare, natural surfaces. Fig. 6 shows the distribution of HISET surface thresholds for the NDVI+NBR index in Europe. The map reveals a diverse and structured distribution of thresholds, spanning from −0.12 to 0.6. 4. Discussion 4.1. SRC evaluation method The framework enables an automatic estimation of the performance of various processing parameters involved in the SRC generation Geoderma 459 (2025) 117340 9