Remote Sensing Data for Digital Soil Mapping in French Research—A Review

Abstract

Soils are at the crossroads of many existential issues that humanity is currently facing. Soils are a finite resource that is under threat, mainly due to human pressure. There is an urgent need to map and monitor them at field, regional, and global scales in order to improve their management and prevent their degradation. This remains a challenge due to the high and often complex spatial variability inherent to soils. Over the last four decades, major research efforts in the field of pedometrics have led to the development of methods allowing to capture the complex nature of soils. As a result, digital soil mapping (DSM) approaches have been developed for quantifying soils in space and time. DSM and monitoring have become operational thanks to the harmonization of soil databases, advances in spatial modeling and machine learning, and the increasing availability of spatiotemporal covariates, including the exponential increase in freely available remote sensing (RS) data. The latter boosted research in DSM, allowing the mapping of soils at high resolution and assessing the changes through time. We present a review of the main contributions and developments of French (inter)national research, which has a long history in both RS and DSM. Thanks to the French SPOT satellite constellation that started in the early 1980s, the French RS and soil research communities have pioneered DSM using remote sensing. This review describes the data, tools, and methods using RS imagery to support the spatial predictions of a wide range of soil properties and discusses their pros and cons. The review demonstrates that RS data are frequently used in soil mapping (i) by considering them as a substitute for analytical measurements, or (ii) by considering them as covariates related to the controlling factors of soil formation and evolution. It further highlights the great potential of RS imagery to improve DSM, and provides an overview of the main challenges and prospects related to digital soil mapping and future sensors. This opens up broad prospects for the use of RS for DSM and natural resource monitoring.

1. Introduction

Soils are a key component of the Critical Zone and a limited resource. They are essential to addressing the climate crisis we are facing. Soils are at the crossroads of existential issues such as food and water security and safety, climate change mitigation, and adaptation, human health, sustainable energy, and biodiversity protection [1,2]. Yet, around the world, they are under threat [3,4], even though they are essential for achieving most of the sustainable development goals (e.g., [5,6,7,8]). These threats continue to grow under human pressure [3,4]. To address this problem, policy initiatives are emerging in Europe, such as the Soil Thematic Strategy [9] or the E.U. Mission “A Soil Deal for Europe” [10]. There is an urgent need to map and monitor soils at field, regional, and global scales in order to improve their management and prevent their degradation. This remains a challenge due to the high/complex spatial variability inherent in soils. Providing such information by mapping and monitoring soil properties in a reproducible way is now possible through Digital Soil Mapping (DSM, [11]). The concept of DSM has rapidly developed and become operational thanks to the harmonization of soil databases, advances in spatial modeling and machine learning, and the increasing availability of spatial covariates [12,13], including the considerable increase in remote sensing data.

DSM emerged in the early 2000s [11,14] as a new concept for mapping soil properties and/or soil classes using soil information and spatially comprehensive covariates, which builds on the work of Jenny [15]. In their seminal article, McBratney et al. [11] proposed a framework called the scorpan-SSPFe. This framework included a soil spatial prediction function for soil mapping with spatially auto-correlated errors, often referred to as SCORPAN. As defined by McBratney et al. [11], SCORPAN is based on seven predictive factors of soil types or soil attributes, namely: (1) s: soil, other or previously measured attributes of the soil at a point; (2) c: climate, climatic properties of the environment at a point; (3) o: organisms, including land cover and natural vegetation; (4) r: topography, including terrain attributes and classes; (5) p: parent material, including lithology; (6) a: age, the time factor; (7) n: space, spatial or geographic position.

The prediction functions can use a wide variety of statistical tools, including regression, classification, geostatistics, machine learning, etc. It is not the purpose of this review to detail these methods or their advantages and limitations. We refer to McBratney et al. [11], Grunwald [16], Minasny and McBratney [12,17], Zhang et al. [18], Wadoux et al. [19], Arrouays et al. [20], Wadoux et al. [19] and Chen et al. [21].

A recent trend is that the availability of RS products, and in particular those derived from satellite [22], airborne, or unmanned aerial vehicle (UAV) imagery [23], has increased significantly. Thus, almost all DSM approaches rely on RS data to obtain spatially exhaustive environmental covariates [21], as well as legacy soil, geology, geomorphology, and landforms maps. Ustin and Middleton [24] identified and described 48 instruments and 13 multi-instrument platforms that collected data in the 2000s, were recently launched, or are expected to be launched in this decade. Indeed, some of these are now being considered as a direct measurement of some targeted soil properties [25]. Commercial operations are also entering the market, providing competitive products with high spatiotemporal (e.g., Planet, https://www.planet.com/products/planet-imagery/ accessed on 29 May 2023 or EarthDaily https://earthdaily.com/earthdaily/ accessed on 29 May 2023) and even spectral resolution (e.g., Kuva Space, https://kuvaspace.com/ accessed on 29 May 2023). In this review, we mainly focus on satellite and airborne imagery as we aim to review DSM studies covering rather large areas. Nevertheless, we do not totally exclude examples on UAV when they are relevant for field- or farm-scale soil mapping or monitoring.

French research has a long history in the field of soil RS. As early as the early 1960s, they developed innovative approaches for the photo-interpretation of aerial photographs for soil mapping [26,27,28], visual interpretation, and/or digital processing of RS images. It was stimulated by the launch of the first Earth observation satellite, Landsat, in 1972 and, in 1986, the French Satellite Probatoire d’Observation de la Terre (SPOT). Some pioneering French works include general books on RS applications to the characterization of the biosphere or terrestrial environment (e.g., [29,30,31,32,33,34]). The development of soil RS was notably based on French research focused on the spectral characterization of minerals [35,36,37] and on the soil surface status [38,39,40,41,42,43,44,45,46,47,48]. In the latter case, the aim was to define soil spectral indices such as the widely used brightness index [49,50,51], the Transformed Soil-Adjusted Vegetation Index (TSAVI) [52], and the theorization of the soil line concept [53]. This led to the physical modeling of soil bidirectional reflectance [54,55] thanks to the SOILSPECT model of S. Jacquemoud [56], which was followed by the MARMIT model of A. Bablet [57,58]. The same thing has happened in the microwave and thermal infrared domains [59,60].

Observation and hence, understanding the spectral behavior of soil as a function of its composition and structure has been the cornerstone for establishing relationships between the spatial distribution of soil properties at different scales, soil types, or soilscapes, and RS imagery products in the solar (e.g., [61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77]), as well as in the microwave (e.g., [65,70,71,73]) and the thermal infrared (e.g., [59,78,79]) domains. Several studies have also addressed the mapping of soil color (e.g., [46,51,80]) or soil moisture (e.g., [78,79,81,82,83,84,85,86]). Others aimed to consider the bare soil surface reflectance, roughness, rockiness, or soil cover by vegetation or crop residues for several applications: modeling/mapping of soil crusting, surface water runoff, erodibility or erosion risk (e.g., [40,47,87,88,89,90,91,92,93,94,95,96,97,98]), and salinity [41,48,99,100,101].

Vegetation is a problem for direct RS of some primary soil properties that are easier to capture by images acquired on bare soils [102]. Therefore, several works have attempted to develop vegetation indexes that are nearly insensitive to the soil spectral response, unlike the NDVI (e.g., [103,104]). However, (i) vegetation is one of the main controlling factors of soil evolution (e.g., inputs in soil organic carbon and influence on pH); (ii) vegetation information may also be related to some soil properties that are not easily detected using only soil surface response (e.g., soil depth and available water capacity (AWC)) and thus need to be captured for DSM.

Comprehensive reviews of the use of remote sensing for soil and terrain mapping were made from the French research as early as 1974, albeit in French [31,105], and major reviews were later published in English by Ben-Dor et al. [106,107], Ben-Dor [108], and Mulder et al. [109]. The latter introduced the use of RS as primary or secondary soil data to improve DSM. As RS products have become increasingly available and diverse and also offer higher resolutions in space, time, and spectral range, updated reviews have been published by Demattê et al. [110] and Chabrillat et al. [111], focusing on the main soil characteristics monitored by RS and the perspectives offered by new sensors [112].

The idea of an article gathering the specific contributions of RS to DSM from French research originates from a scientific working group founded in 2015 and named Theia “CNS for Cartographie Numérique des Sols” (DSM for Digital Soil Mapping). This French working group aims to federate the efforts of French research laboratories developing digital mapping approaches for perennial soil properties [113,114] under the supervision of Centre National d’Etudes Spatiales (French Space Agency). We believe that such an approach can be useful not only on a national scale, but that it can also be considered as the first step to implement a similar approach on a global scale.

The objectives of this review are as follows:

• Summarize the main soil properties and threats that have been studied on a large scale over the last decade using RS by the French research community;
• Synthesize the main recent methodological advances of DSM related to the use of RS products in France or elsewhere from French research;
• Highlight the complementarity of the new RS products and the other covariates currently used in DSM.

2. General Considerations on the Relative Permanence of Soil Properties

DSM predicts soil properties that may change over different times, from brief sudden events to millions of years. Time intervals vary with soil properties and processes and are highly dependent on the controlling factors for soil formation and anthropogenic influences. These factors can be rather slow, such as weathering and pedogenesis, or short-lived when related to destructive events, such as deforestation, erosion and flooding, or local contamination. RS technology has evolved rapidly over the last 50 years or so, providing greater diversity in spectral domains, but also a drastic increase in their spatial resolution and revisit time. This improves our understanding of the relative permanence of soil properties, yet it is limited to one week to half a century. In this section, we will refer mainly to the international literature, as many concepts have emerged outside the strict French soil science community.

Soils are continually evolving over time, and the rate of change in soil properties is accelerating due to global changes and continued unsustainable land management [2,4]. The “natural” formation time of soils, from the initial weathering of parent materials to pedogenic evolution, may range from about 10,000 to 2,500,000 years. As a result, the evolution and properties of natural soils mainly depend on age (time), climate (including past climates), parent material, and natural vegetation. It is generally accepted that the “age” factor is least used in DSM, except for some long periods of extreme climate (e.g., glaciation and interglacial periods) [16]. Therefore, the influence of past climate may be underestimated by DSM covariates and is often barely captured by RS products, as they do not cover the duration of most pedogenetic processes. Many soil properties have been affected by drastic changes since the impact of humans, especially by gradual or abrupt land-use changes. This began when humans moved from hunting and gathering to agriculture (Neolithic revolution). This then led to (ongoing) deforestation and the conversion of land to arable land, and more recently to global changes, such as anthropogenic climate change, industrialization, urbanization and increases in impervious surfaces, land-take, and changes in land use/land cover due to food demand. To differentiate the effects of long- and short-term periods on soil variability, researchers have proposed a distinction between genosoils and phenosoils [115]. This concept makes it possible to distinguish certain semi-permanent properties inherited from “natural” pedogenesis from properties that have evolved under human pressure, especially for the history of land use. It has led to further development and discussions of which soil properties can be considered more or less “permanent”, e.g., certain characteristics of deep diagnostic horizons [116], while others may change with land use, management (e.g., [117,118,119,120]), or under other circumstances. This concept is quite similar to the concept of soil capability vs. soil condition, as defined by McBratney et al. [1].

In this study, we will therefore make the assumption that certain soil attributes can be considered permanent (e.g., soil texture, soil mineralogy, and parent material) on the time-scale of RS observations (i.e., 50 years). Unless there is evidence that abrupt or destructive events have occurred, such as severe erosion or salinization. The typical duration of the semi-permanent changes can range from decades to centuries (e.g., soil organic carbon stocks and diffuse contaminant contents) to seasonal ones (e.g., deep horizon structure and permeability, pH, rooting depth, and soil temperature). Much more rapid and/or even nearly instantaneous changes are related to human actions or climate events, e.g., local contamination, nutrients content, topsoil structure and roughness, pH in limed cultivated soils, and, of course, soil moisture.

It is therefore important to keep these different times in mind, especially when dealing with soil data and/or RS acquisition dates, or the time-scales relevant to monitoring. Note that some information gathered by RS, and even by other DSM covariates (i.e., relief) might not reflect the variation of SCORPAN covariates over large time-scales.

3. Developments Related to the Extraction of Soil Properties from Remote Sensing Data

3.1. Use of RS Data as a Substitute for Soil Properties’ Measurement

In many cases, RS data are directly used as a measure of soil properties, which Mulder et al. [109] called the “primary soil data source”. Bare soil reflectance, for example, has been correlated with many topsoil properties measured in the lab. A model or a spectral index is built to predict certain properties, a given behavior, or a threat to the soil from its spectral reflectance at one or more wavelengths. Thus, many attempts have been made to determine soil moisture [121], soil carbon [122], soil clay content [123,124], or soil salinity [125].

In many cases, therefore, RS is not incorporated into a full SCORPAN model, but directly used as a predictor of the soil property of interest through generalized linear models, partial least squares regression (PLSR), geostatistics, machine learning, or other prediction tools (e.g., [126,127,128,129]). In this case, uncertainties and biases in the predictions should be noted, and care should be taken not to apply them outside their domain of validity. If the results are satisfactory, i.e., if RS data capture most of the variability of the soil attribute in a given area, it may not be necessary to incorporate them into a SCORPAN model (e.g., [123,124,130,131,132]).

Note that these considerations do not claim that working on the field and taking samples are not essential. Working on the field is a crucial step in DSM and RS for several reasons:

−

Lab measurements performed on samples taken into the field are essential to calibrate RS and DSM predictions, and they must be representative from both feature and spatial domains (e.g., [12,106,107,108,109,110,111,112,126,127,133,134,135]).

−

Using more high-resolution spatial covariates does not necessarily improve the prediction accuracy (e.g., [12,19,20,136]).

−

Validation should involve appropriate sampling in the field (e.g., [12,17,18,19,20,21,123,126,127,131,132,137]).

−

The lack of field data is often one of the main limiting factors of RS and DSM prediction performances (e.g., [13,21,138,139,140,141,142,143]).

−

Soil knowledge acquired in the field is useful to improve and assess the RS and DSM prediction performances (e.g., [17,19,20,21,127,138,139,144,145,146,147,148]).

3.2. Most Studied Soil Properties in French Remote Sensing Research

In order to have an overview of the most studied soil properties in French RS research, a dedicated bibliographic search was performed using the Clarivate Web of Science. The query returned 2426 articles, which shows a high productivity. Figure 1 shows the relative percentage of the main soil properties selected, and we decided to focus on those properties that exceed 5% of the total number of French research papers. We recognize that there is no guarantee that the search is exhaustive, given the query process. However, to our knowledge, this is the best effort of such a systematic literature search related to French soil RS research.

The most studied soil property is soil moisture (usually topsoil moisture), which has been studied either in the solar domain (e.g., [39,57,58,59,78]), the thermal infrared domain [59,60,78,79], or the microwave domain (e.g., [84,86,149,150,151,152,153,154,155,156,157,158,159]).

Soil moisture is one of the most studied soil properties in France and can be mainly related to the ESA SMOS (Soil Moisture and Ocean Salinity) mission. The origin of SMOS goes back to the early 1990s. As part of the development of passive microwave radiometers by CNES, experiments were conducted over crop fields in 1991 and 1993 in the remote sensing facilities of INRAE in Avignon. Wigneron et al. [160,161] demonstrated the possibility to simultaneously extract soil moisture (SM) and Vegetation Optical Depth (VOD), a parameter related to the vegetation water content, from multi-angular and bi-polarization passive microwave (MW) observations. Shortly after, Wigneron et al. [133,162] defined a retrieval algorithm based on the inversion of the L-MEB model (L-band Microwave Emission of the Biosphere). This algorithm, which is well suited for a spaceborne mission, was the core of a research proposal submitted to the ESA Call for Earth Explorer Missions. The proposal was accepted by ESA in 1999 (phase A), leading to the Soil Moisture and Ocean Salinity (SMOS) mission coordinated by Y. Kerr, its principal investigator [163]. The satellite was launched at the end of 2009 as part of the ESA Living Planet Programme. In parallel, many scientific advances have been made in the domain of active microwaves (radar) with the support of CNES. It is more difficult to extract SM from active MV sensors (i.e., scatterometers and Synthetic Aperture Radar (SAR)) than from passive sensors (i.e., radiometers). In active RS, measurements must take into account complex extinction effects due to vegetation structure and complex scattering effects due to soil roughness. However, while the spatial resolution of passive MW radiometers is limited to 35–50 km, that of active SAR systems can reach a few meters. Thus, the SM products retrieved from passive MW with low spatial resolutions mainly have applications in hydrology, climate studies, and meteorology [153], while the SM products retrieved from active systems can also find applications in agriculture [164].

The soil moisture retrieved by passive or active MW instruments corresponds to surface SM (1 to 5 cm in depth). Several studies have shown that root zone soil moisture (RZSM) can be inferred from time series of surface SM based on model assimilation or from simple temporal filtering methods [165].

Optical and thermal infrared (TIR) observations can also be used to map surface SM at high spatial resolution, but to date, the research is not sufficiently mature and no operational RS product exists [45,166,167,168]. Surface temperature measured by TIR sensors is related to evapotranspiration and can be used to estimate RZSM [169,170]. SM is an important parameter in itself, but it can also be used to estimate permanent soil properties. For example, temporal variations in SM can be used, through assimilation, to retrieve soil parameters such as soil hydraulic properties [134].

Soil roughness is an important characteristic that plays a key role in understanding/modeling different processes (e.g., slaking, runoff, infiltration, and erosion). It also has an important effect on the remote sensing signatures observed by satellites: it can be estimated because it is a parameter of interest in itself, or because its effects need to be corrected for a better estimation of other variables of interest (such as SM in the MW domain).

Soil roughness, or micro-topography, is an important influencing factor on the soil spectral response to a degree similar to that of “chromophores” [65], and is therefore a key parameter in modeling soil surface reflectance [56]. In addition, it is an indicator of the water infiltration potential, soil slaking, and erosion risk [98,171]. Several methods have been designed to measure soil surface roughness. Their main disadvantages were the time-consuming and tedious nature of the measurement in the field. Bretar et al. [172] and later Gilliot et al. [173] developed an easy-to-handle and fully automatic photogrammetric method to derive digital elevation models (DEMs) of the soil surface with millimeter resolution from multi-angular photographs captured in the field.

In the active MW domain, soil roughness has a significant effect that must be corrected to estimate other parameters of interest, such as SM and VOD. In the passive MW domain, soil roughness contributes to increased MW emissions from the soil. It is related to both micro-scale (~1 cm) and macro-scale (~100 m) topographic effects. The first global map of the observed roughness parameter from SMOS was produced by Parrens et al. [174]. The new generation of freely available radar sensors (e.g., Sentinel-1 since 2015) allows the scientific community to collect images with a revisit time (<6 days) and spatial resolution (~10 m) suitable for hydrological and agronomic applications, on local or regional scales [164]. At these scales, soil roughness mapping is mainly performed on bare agricultural soils. The radar signal increases with the roughness (root mean surface height, RMSH) according to an exponential law and becomes constant after a certain roughness threshold. Several studies show a faster saturation of the radar signal with the soil roughness when the wavelength and/or the incidence angle is small [93,175].

The third most important soil property is primarily related to soil organic carbon (SOC). This is not surprising considering that soils are the largest terrestrial organic C stocks and the question of their role in climate mitigation and adaptation is a hot topic [176,177,178,179]. Reviews on SOC mapping using RS are recent [180,181], especially using satellites [181]. Most studies rely on observation in the solar domain, at scales ranging from the field to regional, and finally, to national or even transnational scales. Most approaches have been conducted on bare, cultivated soils: the higher the scale, the lowers the accuracy. The prediction performance varies according to the type of soil [123], the date of acquisition [182], and the presence of crop residues on the surface [22,183]. As far as French research is concerned, RS studies have been conducted at local or regional scales and none of them have yet specifically addressed the use of RS at the national scale of metropolitan France. Studies at smaller scales have mainly relied on UAVs [173] and airborne imagery [184] at a very high spatial resolution (≤2.5 m), while at larger scales, most studies have relied on high-resolution satellite images: Hyperion [185], SPOT [186,187], Sentinel-2 [123,182], and/or Sentinel-2 with the joint use of Sentinel-1 [128,129]. The accuracies found for mainland France range from 2 g kg⁻¹ to 6 g kg⁻¹, leaving open the feasibility of monitoring SOC according to agricultural C-stocking practices. This feasibility, especially the influence of the composition of the soil surface, is the subject of studies at the French and European levels.

Soil color is an important soil property that is directly related to the reflectance of soil in the visible [46,51,80], which depends mainly on the so-called “chromophores” [108]: soil organic carbon content, Fe sesquioxydes, calcareous content, moisture, and to a lesser extent, soil texture and mineralogy. It is widely used in conventional mapping and soil taxonomy. The Munsell color coordinate system was related to visible bands by R. Escadafal [46,80], resulting in the cluster diagram of soil lines versus Munsell color. In our literature database, it was not always easy to distinguish soil color from soil brightness and soil albedo. While the term “soil color” is often used in soil science, the RS community that focuses on vegetation monitoring has investigated color-related soil surface properties such as soil brightness derived from the soil line concept or soil albedo, mainly derived from inversion of kernel-driven reflectance models on selected directional reflectance data [188].

Topsoil texture defined as the relative proportion of clay, silt, and sand is one of the most important characteristics of soils. It determines many of their physical and chemical properties and behaviors (e.g., water retention, mechanical behavior, infiltration capacity, nutrient and pollutant concentrations, friability, and trafficability). Attempts to map soil texture have shown that the best results were achieved in the visible, near infrared, and shortwave infrared from hyperspectral airborne images (e.g., [124,126,189,190]) or from multispectral satellite images [130,191,192,193]. These works have mainly focused on clay texture and/or clay content estimation due to the specific absorption of clay minerals around 2200 nm [194]. As observed for SOC estimation, the texture prediction performance of topsoil varies with the date of acquisition and the presence of crop residues on the surface [193].

Finally, in the early stages of the use of RS for soil mapping, RS was often used to delineate/identify soil types or soil classes (e.g., [61,63,64,65,66,67,68,69,70,71,72,73,75,76,195]). At the initiative of French geographers, in particular E. de Martonne [196], these approaches were initially based on the visual interpretation of aerial photographs, most often in photogrammetric mode for a better identification of morphological units [26,27,28]. We can introduce here the notion of soilscape, defined as “a landscape unit including a limited number of soil classes that are geographically distributed according to an identifiable pattern” [197]. A spatial layer of the soil factor can then be analyzed in terms of such identifiable patterns. Digital processing of single-date raster images using a per-pixel classifier [75,195,198] or “textural analysis” [74,76] has been proposed to mimic and replace visual interpretation. Some approaches have used digital clustering of soil profile data [199], through the DIMITRI algorithm [61,200] or hierarchical clustering [72] of semantic information about soil map units. The DIMITRI algorithm was later incorporated into the OASIS algorithm [201], whose objective is to identify repeated patterns in the image by calculating the “composition or class frequency vector” of a given square neighborhood in the image. A similar approach has been developed using the CLAPAS (CLAssement des PAysages et Segmentation) algorithm proposed by J.M. Robbez-Masson [202,203]. The main soilscapes of Lebanon were derived from a single-date Landsat image using OASIS [76], while SPOT was used to obtain large vineyard units in the southern Rhone Valley [74]. First applied to non-spectral data [197,204], CLAPAS was then applied to a MODIS time series at four dates to map the main landscape units of Brittany [205,206]. Pixel-based approaches have been especially applied in contexts where the soil is often bare, as in arid environments [195,198] or on Mediterranean vineyards with bare inter-rows [75]. In addition, a pixel-based approach relying on both a SPOT time series and morphometric layers (elevation and derivatives) to map South African terroirs was a step toward DSM formalism. It performed bootstrap selection of the training/validation areas in order to implement regression trees, and then a “Hyperstat” algorithm calculated the modal value in the output map stack and the assignment frequency as an indicator of mapping uncertainty [77].

Figure 2 shows the temporal evolution of the number of articles by soil properties. This trend relates to technological advances (new satellites and sensors) and scientific and political interest in specific soil properties. We observe a remarkable increase for all properties between the late 1990s and the early 2000s. This is likely due to the launch of Landsat 7 and its Enhanced Thematic Mapper Plus (ETM+) sensor, which allowed for more accurate assessments of the land surface [24]. In 2008, the United States Geological Survey (USGS) made all Landsat 7 data available to the public [207], and even free in 2009, which led to a large increase in downloads and probably explains the largest increase in publications after 2010. Studies on some soil properties have steadily increased (i.e., soil moisture), while some seem to have had a rather constant interest over the last twenty years (roughness, temperature, color, and albedo). In particular, although almost absent in 2007, SOC shows a strong increase in interest. This increase is arguably related both to the growing interest in SOC, especially driven by the IPCC [208], and to the launch of satellites that better capture a wide range of SOC-related spectral bands (e.g., [209,210]). Despite their relatively smaller numbers, soil texture studies show a notable increase that is also related to technological advances (imaging spectroscopy (e.g., [124,126])).

3.3. Trends in French Remote Sensing Research in Relation to Threats to Soil

In this section, we focus on articles that explicitly mention soil threats in our corpus of publications (Figure 1). Soil threats are classified according to the proposals made by the EU Soil Thematic Strategy [9,211] and the EU soil strategy for 2030 [10]. The relative importance of soil threats (Figure 3) is related both to their importance in the political agenda and to the feasibility of using RS products to monitor these threats.

The strong dominance of SOC is clearly related to the feasibility of reporting, monitoring, and verifying soil carbon stocks according to article 3.4. of the Kyoto Protocol (e.g., [212,213]). It also links to the COP21 Paris Agreement and the four per Mille initiative [177,214,215]. In addition to advances in temporal frequency, high spatial and spectral resolution RS imagery offers an unprecedented opportunity to study and monitor the spatiotemporal dynamics of SOC changes [111,181].

The second threat is the sealing of soils by urbanization, industrialization, road construction, etc. It corresponds to an almost irreversible loss of many soil functions and is particularly widespread in some regions of France (mainly peri-urban and coastal areas). The fine spatial resolution of the new RS sensors enables them to monitor the expansion of sealed areas over time. In France, the Land Cover Scientific Expertise Center (CES OSO, https://www.theia-land.fr/ accessed on 29 May 2023) automatically produces land-cover maps for metropolitan France from Sentinel-2A and Sentinel-2B data. These maps have a ground sampling distance (GSD) of 10 m. They are used to monitor the rate of soil sealing and/or to take into account the services rendered by soils in land-use planning (e.g., [216,217,218]).

The third threat is erosion, which is a serious problem in French agricultural soils (e.g., [219,220,221]). Soil crusting and slaking have been studied using RS for a long time (e.g., [98,222]). High-resolution imagery is an increasingly used technology for mapping erosion [223]. Soil surface changes can also be captured by DEM series from UAVs [224].

Examples of soil contamination are less numerous, often related to the impact of mining or industrial activities on trace elements (e.g., [225,226]), hydrocarbon and oil contamination (e.g., [227,228,229]), and some radionuclides [230]. Many studies are still based on hyperspectral, fluorescence, or photoluminescence measurements performed in the laboratory or by proximal RS (e.g., [226,231]). There are, however, a few studies using airborne hyperspectral or gamma-ray spectrometry data [227,228,230,232,233,234,235].

The fifth most important threat is soil salinization, which is a very serious problem for soil fertility. It mainly concerns arid and semi-arid regions and is often linked to the quality of irrigation water [4]. Mainland France is not a salinization hot spot, except for a few local areas such as the Rhone delta (Camargue). However, French research is very active in the field of soil salinization RS. This is partly due to the history of France and its relationships with certain countries. This is also related to the fact that salinization crusts are quite easily detectable, even in the visible, and that studies on soil salinity mapping and monitoring have been developing for a long time and are still ongoing (e.g., [99,101,198,236,237]).

Other soil threats are more rarely studied using RS. Biodiversity, compaction, and acidification are not directly related to RS imagery. Indicators of soil biodiversity or acidification can, however, be captured by RS through vegetation or land use. In this case, they are often combined with other covariates in DSM tools.

Figure 4 shows the evolution of the cumulative number of articles per soil threat.

The most striking trend is the almost exponential number of articles on soil carbon loss. As explained above, this trend relates to both policy issues and advances in the feasibility of monitoring SOC. Similarly, we observe a large increase in soil sealing studies, related to population growth, soil sealing, and land-use issues, as well as the opportunities provided by new high-resolution, time–frequency satellites. The increase in soil erosion and contamination studies also relates to global issues and advances offered by new sensors. In comparison, studies on soil salinization remain rather stable; biodiversity is slowly emerging, while acidification and compaction are almost negligible.

3.4. Remote Sensing of Soil Properties as Training Data for Digital Soil Mapping

Similar to laboratory measurements, RS predictions can be used as an “s” factor in the SCORPAN model (i.e., they are used as if they were “measurements” of the variable of interest). In this case, the main challenge is to incorporate the uncertainties associated with these inputs into the SCORPAN model. Note that this issue also applies to laboratory measurements, although it can be assumed that they are generally more accurate than those derived from RS data. Another important limitation of using RS data as direct inputs to calibrate DSM models is that soil reflectance only captures information about the surface layer, while soils are 3D objects. Only radar or gamma-ray spectrometry techniques are capable of observing a soil at depth.

If the prediction does not include a spatial component, and if there are a sufficient number of validation points, an easy and often-used way to verify the absence of a spatial component is to perform geostatistical analyses on the prediction residuals. If there is a spatial structure on the residuals, it is possible to improve the predictions by kriging or find that this trend might be related to another covariate that could be incorporated in a SCORPAN model.

It is important to keep in mind the consequences of uncertainty if these predictions are not only used to predict a value at a given time, but also intend to monitor changes in a soil characteristic. First, cumulative errors at two dates will affect the minimum detectable change (i.e., [238]). Second, it is debatable whether it is better to map soil attribute values at two dates and then infer changes from their differences, or whether it is better to directly map the changes themselves. Both from a spatial and soil process point of view, the spatial structures and controlling factors of a soil attribute at two different dates may not be the same as the changes between two dates. It should be noted that the above remarks are relevant to both for the “direct RS detection of soil properties” and DSM.

4. Incorporating RS as Covariates in DSM

4.1. Soil Property Maps Using Remote Sensing on Bare Soils as Covariates in DSM

Instead of using the soil spectral response to inverse a model prediction of a soil property, an alternative is to incorporate the map of that prediction as a covariate in a SCORPAN model [239,240,241]. These studies are limited by the spatial availability of bare soil images; the influence of perturbing factors such as soil roughness [157], moisture, and vegetation residues on the topsoil [102,242]; mismatches between soil sampling and remote sensing acquisition dates (e.g., [186]); and, in some cases, the low variability and range of some soil property values [123,126,182]. Increasing the frequency of acquisition dates (e.g., with S1 and S2 products) now enables the detection of more bare soil fields and area expansion by mosaicking several images [128] and reduces the influence of soil moisture and crop residues on the prediction accuracy (e.g., [128,129]). Figure 5, extracted from [243], shows an example of RS data on bare soil used as covariates in a DSM model to predict SOC in a region of central France. Interestingly, the covariates derived from RS imagery on bare soil are among the most important in SOC prediction.

In this example, the question of how to combine these soil predictions with a model adapted to usually vegetated areas remains to be solved, in order to obtain a full DSM map of soil properties over a large area including various land uses. This general question applies for nearly all RS data obtained from bare soils. If one wishes a full map of a landscape, RS data on bare soils are almost never available on the whole area, except for deserts.

4.2. Remote Sensing Data as Proxies of Soil Properties Controlling Factors in the SCORPAN DSM Model

4.2.1. Climate

On a global scale, climate is one of the major controlling factors for weathering and subsequent pedogenesis. Actually, the birth of pedology relates to Dokutchaev’s pioneering work on the zonal climatic distribution of soil types in Russia (e.g., [244,245,246,247]). At this scale, climate is also the major controlling factor for SOC contents and stocks (e.g., [177,248]), either through the effect of climate on vegetation-derived SOC inputs (e.g., in the hottest and driest regions) or through the effect of very low temperatures on lowering the rate of SOM mineralization (e.g., in the circumpolar and high-elevation regions).

Figure 6 illustrates the global effect of climate on the world’s carbon stock distribution downloaded from the website https://soilgrids.org/ accessed on 29 May 2023. The original work on DSM at the global scale is from Poggio et al. [249]. As explained in this article, climate is one of the most important driving factor of SOC stocks at the global scale. The map clearly highlights the effect of temperature in the Northern latitudes and in the mountains. There is also an effect of a very wet and hot climate, for example, in Indonesia.

At the scale of mainland France, latitude and altitude are clearly the two main controlling factors for SOC, together with land cover (e.g., [250,251,252]). The effect of altitude was also demonstrated for climatic gradients related to elevation or rainfall in smaller regions (e.g., [238,253,254]).

Rainfall intensity is another major controlling factor for soil weathering [255] and erosion [256]. Finally, climate change is one of the factors accelerating certain pedogenic processes and threats to soil [4].

4.2.2. Land-Cover and Vegetation Characteristics

Land-cover classification based on RS products is a relevant covariate for predicting many soil attributes, in particular SOC. In mainland France, studies have used land-cover maps as covariates to map SOC concentrations or stocks (e.g., [138,250,257,258,259]). Land-cover classification from RS products can be a relevant covariate for predicting many soil attributes. Other studies have also shown that SOC sequestration potential was also related to land cover (e.g., [216,260,261]). Land cover is also clearly related to soil pH as well as parent material [259]. RS products are one of the main tools used for assessing human pressure on soil contamination through the distance to urban or industrial areas and major roads (e.g., [262,263,264]).

Changes in land cover (and the transition duration from one land cover to another) are also direct controlling factors of erosion risks. In a review paper, King et al. [222] summarized soil erosion monitoring and modeling from RS. The first interesting piece of information accessible via RS is vegetation cover, as bare soils are generally the areas that contribute most to runoff (e.g., [265,266,267]). RS has also been used as an indicator of slaking [268]. Hill et al. [269] monitored land degradation, soil erosion, and desertification in a French Mediterranean region and demonstrated excellent accuracy in identifying soil degradation levels using airborne imaging spectroscopy (AVIRIS) and Landsat TM data. For erosion as well, the distributions of land use across the landscape and connectivity between fields are studied using RS and provide useful information on runoff concentration risks [270,271]. RS information on soil roughness, for example using radar data, can also significantly improve erosion prediction and modeling [93,222].

Finally, land cover may sometimes be correlated to some specific local soil properties. For example, in large, cultivated plains such as the Beauce region in France, small, forested plots are often related to very shallow soils [272].

The identification of bare soils uses vegetation indexes such as NDVI (Normalized Difference Vegetation Index). These were designed to assess the amount of green vegetation [273]. As such, they can also be used as proxies of organic carbon inputs to the soil when integrated over rather long periods of time (e.g., monthly or annual averages, [259,260]). In this case, the amount of green vegetation is supposed to be related to the carbon originating from roots and/or crop residues.

Vegetation indexes, such as land surface temperature (LST), as evapotranspiration factors, can also be an indicator of vegetation stress and the soil’s available water capacity (AWC) for plants during dry periods [274,275,276]. In general, the AWC can be calculated using pedotransfer functions (PTFs), established between water retention characteristics and some other spatialized information (e.g., texture, depth or rooting depth, bulk density, and coarse elements, [277,278]). Another way to obtain indirect information about the AWC is to assimilate RS data into vegetation functioning models (e.g., [279,280,281]). Fusion approaches have also been developed to combine local soil property measurements with RS data [282,283,284]. A recent review by Cousin et al. [285] provides an insightful overview of the ways to estimate the AWC using RS data as covariates for DSM AWC prediction from field to regional scales.

4.2.3. Relief, Topography, and Landforms

Relief, topography, and landforms are among the main controlling factors in the SCORPAN model [11,286]. The covariates used in DSM are generally computed derivatives of Digital Elevation Models (DEMs). For global mapping, Sanchez et al. [287] and Arrouays et al. [286] recommended using the Shuttle Radar Topographic Mission (SRTM) DEM at 90 m resolution. The main reason for recommending SRTM was its availability and reasonable resolution considering the number of pixels generated (about 18 billions). Of course, there are now more accurate DEMs available from other RS products, including Lidar for the most accurate ones (e.g., [288,289]). DEM derivatives have proven to be very efficient as covariates in DSM at different resolutions and area sizes, from field to national scales in France (e.g., [138,139,140,144,145,259,272,290,291,292,293]). The integration of DEM derivatives at various resolutions has also proven to be useful in capturing different scales of variability (e.g., [141,144,272]). Note that in some cases, elevation replaces temperature or rainfall/snowfall. It is clear that relief attributes directly reflect or control many factors in soil evolution (e.g., temperature, slope, soil depth, erosion and deposition, and alluvial plains). Conversely, some changes in soil surface topography can be useful covariates to help detect rill erosion [224], landslides [294,295], or the swelling/shrinkage of clay soils (e.g., [194,296,297]).

Note that, although French examples are not dominant in the literature (e.g., [76,129,144,201,202]), the automatic classification of landforms from DEM and other satellite imagery is one of the classical foundations of spatial segmentation, which is a technique often used for large-scale DSM (e.g., [142,298,299,300,301]).

4.2.4. Parent Material

Historically, information about parent material in DSM came primarily from digitizing geological or lithological maps, or soil maps in which the parent materials were described. This is still the case for many studies. However, RS products have proven their ability to detect and map certain parent materials through the spectral signature of outcrops, regolith, or upper soil, especially using airborne gamma-ray spectrometry. This has been shown in France in a large number of situations (e.g., [135,141,232,272,302,303,304,305]). Figure 7 shows four DSM predictions of topsoil texture [305] in a French department, using a large set of covariates and then adding a lithology map, or a gamma-ray spectrometry map, or both. It shows that gamma ray spectrometry captures some rather large spatial patterns linked to the parent material that were not captured by the lithological map.

Applications of RS to detect typical parent materials or mineralogy are quite common in tropical environments (e.g., [109,306]).

4.2.5. Age

According to Chen et al. [21], age is the least-used covariate for medium- and large-scale DSM. It is a rather complicated covariate to interpret because it mixes the duration of weathering and pedogenesis (both its duration and the various climatic and vegetation conditions in which soil formation occurred). Overall, age is quite difficult to infer directly from RS data. Where soils have developed from the underlying parent material, the age can be inferred from lithological or geological information (see Section 4.2.4). However, to our knowledge, there are no RS sensors that are capable of capturing age at the pedogenesis time span. For shorter periods, time can be captured through land-cover change history or other types of measurements. Although we did not find any soil studies using gamma-ray spectrometry or other tools to measure the erosion/sedimentation rates in French soils, some works suggest that this technique may have potential for detecting change related to nuclear accidents (e.g., [307,308]).

4.2.6. Soil Management Practices

In the “o” factor or the SCORPAN model, all organisms are included, although in most cases, “o” is represented only by vegetation indexes or biomass-related vegetation characteristics such as NPP. Among the organisms that have a major effect on soil properties and changes, one of them is obviously humans, through their impact on soil management practices. We refer here to a recent review by Bégué et al. [309], which emphasizes that the potential of RS to characterize cropping practices is not yet optimal, and that the new generation of high-resolution sensors (space, wavelength, and time) should allow for the improved characterization of cropping practices. In another review, Bégué et al. [310] stated that “less than 10% of the publications on remote sensing and agriculture actually focused on cropping practices” in the previous decade. In their 2018 review, they distinguished between three broad groups of cropping practices ([309], Figure 1):

• Crop succession (i.e., the sequence of crops or fallows in consecutive years at the field scale);
• Cropping patterns (e.g., temporal sequence and spatial arrangement of crop, fallows, and landscape features in a particular land area);
• Cropping techniques (e.g., irrigation, organic amendment, crop residue management, tillage, harvesting, duration of bare soil, grazing of moving cattle feed, and implemented on a piece of land).

It is beyond the scope of this study to detail all the sensors, indices, and data they listed in their review to detect all these soil management practices. Figure 8 shows the main cropping practices reviewed by Bégué et al. [309].

One of their main conclusions is “the majority of the studies are exploratory investigations, tested at a local scale with a high dependence on ground data, involving one remote sensing sensor at a time, and are constrained by local knowledge and conditions”. However, they pointed out the high potential related to the increasing availability of RS data suitable for fine resolution and time-scale monitoring (e.g., Sentinel 1 and 2, hyperspectral sensors, Lidar) and the combination of different satellite sensors [311,312]. Research should therefore be conducted to better map these soil management practices. One of the major difficulties was obviously related to the low spectral, spatial, and time resolution. Another outstanding question is what type of information derived from RS should be considered as a covariate for DSM? There is probably not a universal answer to this question; it may depend on the case, area, or soil property.

Nevertheless, we believe that this essential information will be incorporated into future DSM and will provide major advances in soil knowledge. First, it will better filter out some information that is often considered “noise” for some DSM predictions (e.g., soil roughness or dead vegetation residues for SOC). Second, it will provide critical information on the controlling factors of soil changes on the local and landscape scales. Such quantified information may include, but is not limited to, crop rotation [179], grassland management [313], spatial arrangement of crops and land use, linear features such as hedges or drips [314], soil water regimes and their consequences on SOC inputs and salinization, etc. Indeed, a large field of research is still open in order to find the right RS soil management and cropping practices to incorporate in the DSM covariates. A major difficulty could be the diversity of situations in which we will approach large-scale studies. A preliminary large-scale delineation of “agrarian” regions could be a solution, as suggested by Bellon et al. [315]. Furthermore, this typically “multi-scale” RS challenge may benefit from data fusion of different RS data at different wavelengths, resolutions, and revisit periods.

Finally, from the perspective of local soil management and security, this is really crucial, because it can fill the gap between DSM and local farmers and actors, as crop and soil management practices are indeed, a very local way that farmers can try to improve soil conditions [316,317].

5. Use of Remote Sensing to Design the Sampling Strategy for DSM

When RS products are available and relevant in a target area, and if we assume that the SCORPAN model for DSM is efficient, it seems logical to optimize the soil sampling design based on covariate space, and thus to use for designing the sampling, RS products alone or with other covariates. This may be especially useful if the number of soil measurements/observations is constrained. This has led some authors to propose that soil sampling is designed in terms of covariate space. One of the seminal articles in this technique is the constrained Latin Hypercube Sampling (cLHS) proposed by Minasny and McBratney [318]. The cLHS sampling in the presence of covariates is designed to select a limited number of sites in an area.

According to Minasny and McBratney [318], “This kind of sampling should produce a reasonably efficient way of sampling soil and its environment so that the range of conditions are encountered, ensuring a good chance of fitting relationships if they exist”. The last words of this quote, i.e., “if they exist”, are very important to consider. A strong assumption of cLHS is that the selected covariates represent the variability of the target variable. In other words, if a controlling variable is missing from the set of covariates, the sampling design has a high probability of missing part of the variability of the target property.

The main advantage is to optimize the number of samples and thus the sampling effort and cost to capture the main combinations of the covariate feature space. One of its weaknesses is that it is not well suited to some spatial modeling techniques, such as geostatistics and kriging [131]. This is because cLHS, in its process of selecting sampling points, maximizes the feature difference between sampling points. The residuals of the predictions obtained after modeling using a cLHS sampling for calibration, can, however, be analyzed by geostatistics or other validation techniques, on the condition that an independent validation set is available (see an example in Figure 9, left). Note that in the figure extracted from Minasny and McBratney [318], they just use the stratified sampling design as an example compared to cLHS.

Various adaptations of cLHS proposed by different researchers consider accessibility and cost, or better space filling (e.g., [319,320,321]). Space filling has the advantage of not missing situations that could be under-represented by clHS. If, for example, a large area is characterized by similar combinations in the feature space, it might be under-sampled, even if it contains a soil variability that is not captured by the covariates. Wadoux et al. [322] recently studied sampling design optimization for soil mapping with random forest, and concluded that cLHS sampling performed worse than other sampling designs. It should be noted that, in their study, the subsampling design was based on a preliminary selection of an existing set of sample sites that was already a mixed sampling design, i.e., subsampling a regular grid using a method close to cLHS, and with known results (LUCAS Soil, [323]). In other words, the spatial coverage of their learning full sample was already rather efficient, and some covariates had already been used to locate the samples in the grid.

Note that the idea of exploring the cLHS method using covariates on a dense set of already sampled points was also tested by Loiseau et al. [141], who concluded that it was efficient at highlighting situations that were under-sampled in the original soil legacy dataset. Briefly, they tried to predict topsoil particle size distribution using various DSM techniques and a set of covariates. When applying cLHS to the covariates, they found that the initial dataset missed some covariates combinations, which explained why some areas had a bigger uncertainty, or even a bias, in the predictions.

Another interesting characteristic of direct RS prediction is that it can be used as a dataset to simulate the spatial distribution of a soil property and thus be used to simulate different validation sampling strategies in spatial or feature space (i.e., [132]) and highlight their advantages and limitations. This kind of direct prediction, if accurate enough, can also help test different sampling strategies for the training dataset. This test on sampling strategies was conducted by Lagacherie et al. [131] in order to assess their impact on the accuracy of the prediction of topsoil clay content in Tunisia. In particular, they showed that stratified random sampling, which is often recommended for soil map validation [137], could lead to inconsistent results if the number of strata and/or the number of random samples taken in each strata were not large enough. For a recent synthesis about sampling in DSM, see Brus [137].

6. Limitations and Challenges

In this section, we focus on the main limitations of RS products for soil mapping and monitoring and some challenges to improving the quality and relevance of RS imagery for soil.

6.1. Using RS Products to Predict Soil Properties

The use of RS measurements as a substitute for soil analyses is a way to increase the density of soil information at a low cost. However, this surrogate often remains difficult to obtain with satisfactory accuracy and broad spatial coverage. With respect to accuracy, work remains to be conducted to better distinguish the effects of different factors on the spectral response of a soil property. A typical example can be taken for mapping SOC using the spectral signature of bare soil. A large number of limitations include, for example, the following:

• How do we discriminate bare soil from vegetated soil, e.g., what NDVI threshold should we use? Which unmixing method should we use for partly vegetated pixels?
• How do we distinguish bare soil from soil covered with dead vegetation residues? Additionally, what threshold should we use?
• How do we take into account the complex effects of different soil properties at a given time/location (e.g., vegetation, residues, moisture, roughness, SOC, and clay or lime content)?
• How do we extend the coverage of bare soils? Mosaicking is now well developed, but what time period should we take? What are the limitations related to mosaicking? How can we take into account the fact that different dates are usually associated with different situational conditions (e.g., moisture and roughness)?
• In many cropped regions, more and more soil will no longer be bare due to the implementation of cover crops and/or seeding under a vegetative cover.
• How do we take advantage of various RS products and sensors, depending on their resolution (in space, time, and spectral domain), remoteness and related perturbing factors?
• At the field or farm-scale, responding to these questions may imply using several RS sensors. For example, using UAV RS enables the selection of the most adapted dates to obtain bare soil imagery together with minimizing the perturbing factors of soil spectral data. This is, however, not feasible for broad-scale monitoring. The effect of some perturbing factors may also be studied using airborne imagery. Thus, studies comparing laboratory spectra, proximal sensing, UAV, or airborne data to satellite imagery, should be conducted in order to estimate/reduce the errors due to changes in spectral resolution/remoteness and to improve the criteria for selecting the relevant satellites data (e.g., [324,325,326]).

6.2. Most RS Products Only Capture Topsoil Information

Almost all DSM attempts show that the prediction performance decreases with depth (e.g., [21,259]). This is easily understandable because most RS products only capture information about the soil surface and/or terrain attributes at the time of acquisition. For example, it is physically impossible to detect deep soil attributes with visible–near infrared RS data. Some products, however, have the capacity to integrate a response over a less limited depth (radar and gamma-ray). Therefore, one challenge is to better incorporate different sensors that are capable of obtaining information at different depths into the DSM framework. Other information than bare soil spectra can be useful to map soil properties at depth. Gamma-ray spectrometry data, however, comes mainly from low elevation airborne sensors, which are not yet available in many parts of the world. Magnetic surveys are most often used for detecting soil pollution (e.g., [327,328,329]). They may also, however, reveal lithology, weathering, and water regime of soils (e.g., [330,331,332]). Yet, to our knowledge, there are no studies using airborne magnetic data as covariates in DSM.

Another challenge is to integrate some soil knowledge into the DSM. For example, it is well known that in many cases, soils are organized in topo-sequences or “catenas”. This is the case of erosion catenas from a plateau to a gentle slope. At the top of the eroded portion of the catena, one can obtain information about deeper horizons present on soils developed on the plateau, as these underlying horizons may have been exposed at the surface by erosion. The same inferential information can be incorporated into DSM in the case of steep slopes, where outcrops can provide information about the parent soil material. In this case, a sampling design based on topo-sequences and available high-resolution RS imagery can be very efficient (e.g., [333]).

This type of integration of soil knowledge into the DSM is not trivial because it requires a soil scientist to know the soil spatial organization and the scales involved in order to incorporate such information into the DSM. In some cases, these rules are described in conventional soil maps, and a challenge is to take them into account to disaggregate complex soilscape map units into more precise classes (e.g., [146,147,334,335,336]).

Finally, a solution often chosen to address the limitations of bare soil RS and to obtain information on some deep soil properties is to incorporate other RS data, such as vegetation indices or land-cover changes, as covariates in a DSM model (see Section 4.2, and further, Section 7.2.2). At the global scale, Padarian et al. [337] recently used a semi-mechanistic model combined with an MODIS time series of changes in land cover to model changes in the global SOC stocks from 2001 to 2020. Though this approach is still not perfect, because it does not take into account climate change and some inherent soil properties, it provided an estimate of the global losses of SOC stocks, mainly due to the conversion of forests to cultivated lands.

6.3. Relative Permanence of Soil Properties and Revisit Time for Soil Monitoring

In Section 2, we described the relative permanencies of soil properties. One question is how these relative permanencies relate to current and future satellite revisit times. Indeed, it is a great advance that satellite revisit times are becoming shorter, especially with some of the commercially available data products with multiple weekly flybys.

Some applications might expect us to be able to deliver soil information corresponding to this time span. The challenge is the number of cloud-free days due to weather conditions, especially in the northern or mountainous regions of France. Even when there are no clouds, we may be confronted with shadow effects, and the topographic image corrections may not be accurate enough to monitor mountainous regions.

Another challenge is dealing with the observed spatial and temporal heterogeneity that is inconsistent with our understanding of the relative permanence and natural variability of soil properties. Finally, how do we translate this into effective sampling designs for soil monitoring? At field or farm-scale levels, UAV imagery helps design more efficient field sampling designs for monitoring, especially for SOC. It can also increase the density of observations in space and time and help to select the relevant soil signals and filter/correct some perturbing factors. One challenge is to combine detailed surveys and RS products (from proximal sensing to UAV, airborne, and satellite data) in order to better calibrate the models obtained by hyperspectral and high-resolution satellite sensors and better design and optimize field sampling.

A final challenge is to make a clear distinction between what the changes in carbon storage and the changes in carbon sequestration are [179,215,260]. A promising way, at least for topsoil, could be to relate soil spectra to pools of SOC with different residence times. Several studies have shown the potential of laboratory Vis–NIR and MIR spectroscopy for an efficient estimation of SOC fractions, SOC monitoring, and modeling (e.g., [338,339,340]). To our knowledge, the feasibility of upscaling this kind of relationship to RS data needs to be assessed. This may increase the usefulness of RS to map and monitor SOC pools.

7. Main Progresses, Perspectives, and Prospects

7.1. The Increasing Availability of Remote Sensors and of Their Spatial, Spectral, and Temporal Resolution over Time

Taking the first Landsat satellite (1972) as a starting point, the last fifty years have been characterized by the launch of an increasing number of Earth observation satellites with different spatial coverages and resolutions, spectral bands, and revisiting times. Meanwhile, a large number of airborne spectral sensors have been tested on smaller spatial coverages in order to explore or simulate the potential of new sensors, to obtain faster images and shorten revisit times. These technological advances have improved the DSM in several ways:

The increase in spectral range and spectral band resolution has led to a drastic increase in the number of soil properties that can be predicted using RS (e.g., [43,107,341,342,343,344,345]).

The increase in spatial resolution allows the mapping of increasingly fine spatial distributions of soil properties and their controlling factors (e.g., [325,343,346]).

Shorter revisit periods allow for the better detection of changes in soil properties over shorter time periods and better capture of soil management effects (e.g., [347,348,349,350]).

Increasing the number of RS products over time allows both the design of composite products and a shift from static mapping to monitoring certain changes (e.g., [128,326,351]).

Merging data from RS products generally provides better predictions than using a single RS product (e.g., using a time series to increase both the prediction accuracy and map coverage [352], using synergistic information to retrieve soil information [353] or fusing products to increase both spatial and temporal resolution (e.g., [354,355,356,357])).

Several missions with high signal-to-noise ratios, high spatial resolution, but limited coverage per day have been launched recently: the Italian PRISMA (PRecursore IperSpettrale della Missione Applicativa, launch 22 March 2019, 30 m; [358,359]), the German EnMAP (Environmental Mapping and Analysis Program, launch 1 April 2022, 30 m, [360]), and the Japanese HISUI (Hyperspectral Imager Suite, launch 5 December 2019, 30 m; [361]). The Italian–Israeli SHALOM (Spaceborne Hyperspectral Applicative Land and Ocean Mission, launch 2024, 10 m [362]) is also expected to be launched soon. In addition, some global mapping missions are planned or under study for the coming years, with variable spatial coverage and different ground sampling distances, such as the former HyspIRI mission [363,364] and the Sentinel-10/CHIME (Copernicus Hyperspectral Imaging Mission for the Environment) satellite proposed as an ESA candidate mission [365]. This wide range of satellites will represent a major step toward multi-resolution and multi-precisions soil mapping. Virtual satellite constellations (LANDSAT/S2, PLANET, EarthDaily…) are promising tools for the future. One of the challenges is to harmonize the data because the sensors, acquisition conditions, and even the processing chains are different.

7.2. The Increasing Importance of RS Data in DSM

The trends we have observed and the future of RS in DSM can be schematically divided into the following two categories:

• The use of RS products as surrogates for in situ measurements.
• The incorporation of RS products as covariates for DSM.

In this section, we summarize the main conclusions that can be drawn from this review.

7.2.1. Use of RS Products as Surrogates for In Situ Measurements

As developed in Section 3.4., and already highlighted by Lagacherie and Gomez [240], in some cases, “the estimations of soil properties at each pixel of a hyperspectral image is assumed to be precise enough to be considered as a measurement of a soil property and to be used as such for feeding DSM models”. The huge advantage of these situations is that the spatial sampling is much denser than what is usually available for the soil data to train DMS models. A better coverage of spatial soil variability, especially short-scale variability, which is often not captured by existing soil data, can be expected.

Moreover, there are still large spatial gaps in the available soil information to run DSM models around the world [13]. Therefore, these myriad surrogates for in situ measurements pave the way for both filling gaps in learning data and revealing finer spatial structures. These surrogates have proven to be efficient for some soil properties in the case of airborne hyperspectral imagery (e.g., [124,126,132,240,282]). However, we need to consider that these surrogates have uncertainties, and that the relationships that enabled the estimation of these surrogates should not be applied outside of their domain of validity. This is a big challenge, as soil data are still missing in large parts of the world.

7.2.2. Incorporation of RS Products as Covariates in DSM

We have seen that RS products can be used as covariates in DSM models (see Section 4.2). This is primarily because many RS products not only reflect soil surface properties, but can provide spatiotemporal information about two categories of information that can be used in DSM.

RS products can provide information on the controlling factors of soil formation and changes. This is particularly the case for vegetation indexes, which provide indirect information about SOC inputs to soils. Conversely, soil moisture and temperature products can provide indirect information on the mineralization rate of SOC. Land cover and changes in land cover are related to different soil outputs and inputs that can influence pH, C and N cycles, nutrients, and possible sources of contamination. Airborne gamma-ray spectroscopy can provide information on the soil depth, presence of peats, soil weathering, and the nature and mineralogy of the parent material, or help better map clay content and/or mineralogy. The high-resolution DEMs derived from RS are essential to capture the different scales of soil processes and soil erosion/redistribution in the landscape at different scales. They are also very useful information for water redistribution and thus water regimes.

One big challenge with the multiplication of RS imagery is to conduct the right selection of covariates to avoid overfitting and respect the parsimony and relevance principles in DSM. Several ways of selecting covariates have been proposed and cited in this review. We must keep in mind that among the criteria used for selecting covariates and/or aggregating them in synthetic ones, we should take into account not only parsimony and redundancy, but also their meanings in terms of soil science. Conversely, another big challenge of DSM and data-driven modeling is to use them as knowledge discovery tools that may reveal or help with the understanding of new soil processes at various scales. In other words, RS and DSM are not only tools that can extend the spatially statistical relationships established on a set of leaning variates; they should also be able to discover new soil processes.

Although still underutilized, information on soil management practices at the field level is probably one of the most promising covariates with the development of hyperspectral sensors with high spatial and temporal resolution.