# The utility of remotely-sensed vegetative and terrain covariates at different spatial resolutions in modelling soil and watertable depth (for digital soil mapping)

James Taylor, Frédéric Jacob, Mauricio Galleguillos, Laurent Prevot, N. Guix,

Philippe Lagacherie

► To cite this version:

James Taylor, Frédéric Jacob, Mauricio Galleguillos, Laurent Prevot, N. Guix, et al. The utility of remotely-sensed vegetative and terrain covariates at different spatial resolutions in modelling soil and watertable depth (for digital soil mapping). Geoderma, 2013, 193-194, pp.83-93. 10.1016/j.geoderma.2012.09.009. hal-02650587

HAL Id: hal-02650587

https://hal.inrae.fr/hal-02650587

Submitted on 4 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The utility of remotely-sensed vegetative and terrain covariates at different spatial resolutions in modelling soil and watertable depth (for digital soil mapping)

J.A. Taylor<sup>a,*</sup>, F. Jacob<sup>b</sup>, M. Galleguillos<sup>a</sup>, L. Prévot<sup>a</sup>, N. Guix<sup>c</sup>, P. Lagacherie<sup>a</sup>

Digital soil modelling and mapping is reliant on the availability and utility of easily derived and accessible co-variates. In this paper, the value of covariates derived from a time-series of remotely-sensed ASTER satellite imagery and digital elevation models were evaluated for modelling two soil attributes — soil depth and watertable depth. Modelling was performed at two resolutions: a fine resolution (15 m pixels) that relates to the resolution of the ASTER Visible-NIR bands, and a larger resolution (90 m pixels) that relates to the resolution of the thermal bands of the ASTER imagery. Upscaling to a larger pixel size and downscaling to a smaller pixel size were performed to adjust the covariates where necessary. A regression tree approach was used to model soil depth and watertable depth, recorded as a binary 'deep' or 'shallow' response, using the ASTER imagery-derived covariates and digital terrain attributes (DTAs). Modelling was performed at a single spatial resolution (15 or 90 m pixels) using the imagery-derived covariates only, the DTAs only, or a mixture of both. A multi-resolution model was also generated, by using both imagery-derived covariates and DTAs at both resolutions. When mixed with the DTAs, the imagery-derived covariates helped explain the uncertainty (variance) in the soil depth data but not in the watertable depth. The ASTER-derived down-scaled evapotranspiration-based covariates were of particular significance in the soil depth modelling. Watertable depth was best explained by models that used DTAs at a smaller pixel size. Information on veg-etative growth was neither superior nor complementary to information on terrain for modelling watertable depth. Using a multi-resolution model significantly improved the modelling of soil depth but not of watertable depth. The effect of covariate and modelling resolution on model performance is discussed within the context of the GlobalSoilMap.net project.

## 1. Introduction

The prediction of subsoil attributes is generally more difficult than topsoil attributes because many of the environmental covariates used in digital soil modelling and mapping (DSMM) only generate a topsoil or surface response (e.g. Visible and Near-infrared (Vis-NIR) imagery and gamma radiometry). Meanwhile, sensors that are capable of subsoil penetration (e.g. electro-magnetic induction (EMI) and ground-penetrating radar (GPR) sensors) often produce a signal that is difficult to interpret or to deconstruct into individual subsoil attributes and/or different subsoil layer responses. However, in agricultural production situations, mapping subsoil properties, particularly subsoil constraints, is as important if not more so, than mapping topsoil properties. Indeed, information relating to soil depth, depth to an impeding layer or depth to a watertable is critical as this, together with

texture, determines the potential available soil moisture that drives crop yield potential and run-off/infiltration partitioning.

In this article, the issue of identifying relevant covariates and spatial resolution for modelling soil depth and watertable depth is investigated. There are several possible approaches to modelling and mapping soil depth and watertable depth. Most commonly, both are modelled directly from elevation and digital terrain attributes (DTAs) (e.g. Finke et al., 2004). Several geophysical sensors (e.g. EMI, GPR) are able to generate subsoil responses that can be related to soil depth under certain conditions (Bramley et al., 2000). Such sensing technologies have been used for subsoil modelling/mapping, either with or without terrain information (see Grunwald, 2009 for examples). An alternate solution to deriving a direct subsoil response is to map the topsoil response(s) and then use a soil inference system to predict subsoil properties from the topsoil response(s). However, to the authors' knowledge, this has not been done for soil depth or watertable depth. Another alternative, and a focus of this work, is to incorporate a vegetation-soil inference system into the modelling/mapping. The type and vigour of vegetation

<footnote><sup>a</sup> INRA, UMR LISAH, Bât. 24, 2 Pl. Pierre Viala, 34060, Montpellier, France</footnote>

<footnote><sup>b</sup> IRD, UMR LISAH, Bât. 24, 2 Pl. Pierre Viala, 34060, Montpellier, France</footnote>

<footnote><sup>c</sup> VetAgro Sup, Clermont Université, BP 10448, F-63000, Clermont-Ferrand, France</footnote>

<footnote>* Corresponding author. Tel.: +1 716 792 2800.
E-mail address: james.taylor@cornell.edu (J.A. Taylor).</footnote>
growth at a given site will be driven, amongst others, by the soil properties, and predominantly subsurface soil properties, at the site. It should therefore be possible to infer subsoil information from the plant response, particularly in monocultures.

Vegetation type and vigour can be determined from remotely-sensed multi- and hyper-spectral imagery. Whilst the use of vegetative information in soil mapping is not new (McBratney et al., 2003), the use of vegetative covariates, such as vegetative indices derived from imagery band information, in quantitative digital approaches has been somewhat limited to date. Vegetative covariates have generally only been derived from the visible and near-infrared regions of the electro-magnetic spectrum (EMS) and are usually only used empirically within models. A survey of recent DSMM papers (Grunwald, 2009) found that 16.7% used a vegetation response from remotely sensed imagery in the digital soil model. However none of these approaches tried to model or map either soil depth or soil moisture content. The recent review of Mulder et al. (2011) also found no direct application of Vis-NIR imagery to the mapping of soil depth or moisture. Mulder et al. (2011) observed that 'in regional studies NDVI (normalised difference vegetation index) data have been related to soil type patterns rather than to specific soil properties'.

Since 1999, high-resolution (<100 m pixels) images in the Thermal-IR (TIR) section of the EMS have been available to researchers via the ASTER and LANDSAT 7 satellite sensors. These TIR bands have been used to directly map surface properties over bare soils (Goosens et al., 1999; Martínez-Montoya et al., 2010). However there have been no reports of the use of these TIR bands in a vegetation-soil inference system in the DSMM literature to our knowledge. TIR imagery is able to provide information regarding the temperature of a vegetation canopy. If crop and climatic conditions are known, then the image TIR band observation(s) can be translated into a measure of crop water stress at the time of image acquisition. Recent works have explored the possibilities of combining thermal infrared measurement with surface energy balance modelling to retrieve plant water stress and relate it to soil moisture or soil hydrodynamic properties (Boulet et al., 2009; Crow et al., 2008; Olioso et al., 2005). Meanwhile, Galleguillos et al. (2011a,b) demonstrated it was possible, over Mediterranean vineyards, to map evapotranspiration (ET) from ASTER Thermal-IR observations by using simplified energy balance models.

In this study, vegetative covariates derived from a time-series of images will be used. Previous DSMM studies have generally focused on single images or interpreting images within a time-series on an individual image basis (Carre and Girard, 2002; Liu et al., 2008; Meirik et al., 2010; Rivero et al., 2007). This is problematic for development of a digital soil model as the observations from an individual image can only be considered as a snap shot in time. It does not indicate how conditions evolve in a dynamic environment. Modelling and mapping with these 'snapshots' makes repeatability difficult over time or in a different spatial context. Using covariates derived from a standardised time-series captures trends in the imagery rather than focusing on the actual observation values, which should produce more robust covariates for both temporal and spatial extrapolation of the results. Further, a time-series of remotely-sensed ET images throughout the season is expected to provide information on soil moisture that in turn relates to soil depth and watertable access. The same applies for vegetative indices that relate to canopy extension and vegetation growth, although their relationship with soil moisture is not as direct as ET.

The primary intention of this study is to explore how well different combinations of covariates at two distinct spatial resolutions explain the variability in soil and watertable depth in a catchment. The objective is not to generate digital soil maps of the catchment, but to gain a better understanding of the utility of potential covariates in DSMM. This will be initially undertaken by modelling soil depth and watertable depth with only imagery-derived covariates or only DTAs. Subsequently, models with both imagery-derived covariates

and DTAs will be generated to identify if there is an advantage to selecting amongst the two sets of covariates. Since the original imagery used was obtained at two different pixel sizes (15 m and 90 m square pixels), the models will be run at these two different scales, with the relevant imagery either upscaled to a larger pixel size or downscaled to a smaller pixel size where required. Two digital elevation models (DEMs) were also available at approximately the same resolution as the imagery (10 m and 80 m pixels) so that the DTAs can be calculated at the relevant scale without the need for up/downscaling. Since one of the models will be performed at the same scale (90 m) as the current GlobalSoilMap.net initiative (www.globalsoilmap.net), this investigation will also provide some information on the implications for DSMM at this scale and a smaller pixel size.

## 2. Methods

### 2.1. Site description

The study area is the Peyne watershed (43.49°N, 3.37°E) located near Pezenas in the Languedoc-Roussillon region, southern France. This 65 km² watershed is mainly under vineyards (~70% of the watershed area), which are predominantly non-irrigated (~90% of vineyards). The remaining area includes other crops, forest, native scrubland and urban zones.

The watershed is very variable in lithology and pedology. The Peyne watershed is underlain by heterogeneous Miocene marine and lacustrine sediments, i.e. marl, limestones and calcareous sandstones, which were partly overlain by successive alluvial deposits ranging from Pliocene to Holocene. Recent volcanic activity and the localised transport of colluvial material along slopes add to the pedological complexity of the region. Elevation in the watershed ranges from 20 to 230 m above sea level (Guix-Hébrard et al., 2007).

### 2.2. Soil and watertable depth data

Information on soil depth and depth to the watertable during the 2003-04 seasons was available at 41 sites from a previous study (Guix-Hébrard et al., 2007). The depth of interrogation for soil was limited to 2.5 m with 20 sites having profiles deeper than the depth of interrogation. Of the other 21 sites, the soil depth was <2 m (range [0.6-1.8 m]). Therefore there was a distinct separation between the two groups and the soil depth was divided into a binary categorical variable of 'deep' (>2.5 m) and 'shallow' (<2 m). Similarly, the watertable depth data was converted into a binary categorical variable of either 'deep' or 'shallow' watertables. Since the depth to watertable is temporally variable (Finke et al., 2004), the depth at the start of the growing season (April/May) and at veraison (end of July) were considered in segregating the watertable depths. 'Deep' watertables were considered to have a depth >3 m at the start of the growing season and a depth >3.9 m at veraison. All other sites were characterised as having a 'shallow' watertable. The agreement between the categorical soil and watertable depths, as obtained by using a contingency table analysis (JMP 8.0.1, SAS Institute Inc., Cary, NC, USA), was not strong ($\kappa$ = 0.36, $n$ = 41).

### 2.3. Imagery

A time-series of ASTER images was collected over the catchment throughout 2007 and 2008. However, due to both the revisit frequency of the TERRA satellite and the weather conditions (i.e. cloud cover); this was an irregular time-series. The available imagery was subset into two groups. The 2007 analysis was confined to a May 10th and September 15th image that represented vine status at the start and end of the growing season. In 2008, the analysis was limited to five images acquired over a six week period from June 22 to July 31. This period covers the key phenological stage of veraison, and is a period
where water stress in vines is increasing as the hot dry Mediterranean summer takes effect (Taylor et al., 2010). The ASTER images contain bands in the Visible (Vis), Near-Infrared (NIR) and Thermal-Infrared (TIR) sections of the EMS. The Vis-NIR imagery were obtained at 15 m pixels whilst the TIR where collected at 90 m pixels. All images were assessed for accuracy of geo-location and geo-rectified if necessary.

### 2.3.1. NDVI

Normalised Difference Vegetation Index (Eq. (1)) (Rouse et al., 1973; Tucker, 1979) was calculated using the ASTER Red and NIR bands on the original pixel resolution (15 m). The individual red and NIR bands from the imagery were upscaled to 90 m pixels (the same resolution as the TIR imagery) by applying an averaging filter over a 6*6 pixel window. The upscaled NDVI was then calculated from the upscaled bands.

$$NDVI = (NIR - R)/(NIR + R). \qquad (1)$$

### 2.3.2. Evapotranspiration

The ASTER TIR imagery was converted into evapotranspiration (ET) estimations following the process outlined in Galleguillos et al. (2011a,b). The S-SEBI model (Roerink et al., 2000) permitted the derivation of daily ET (ETd) from the ASTER imagery, where the latter captures contrasts driven by evaporation processes. The approach is twofold. Evaporative fraction (Λ), the ratio of latent heat flux to available energy, is computed from the differences between pixel temperature (Ts) and maximum (Tmax) and minimum (Tmin) temperatures within the corresponding albedo class. Next, assuming the instantaneous Λ at satellite overpass is equal to the daily Λ, and neglecting daily soil heat flux, ETd is derived by extrapolating, at the daily scale, instantaneous net radiation at satellite overpass Rni, through the ratio Cd = Rnd/Rni, where Rnd is daily net radiation. This process is expressed in Eq. (2) (L is the latent heat of evaporation) and further details are available in Galleguillos et al. (2011a).

$$ET_d = \Lambda \frac{R_{nd}}{L} = \frac{T_{max} - T_s}{T_{max} - T_{min}} \frac{Cd_i R_{ni}}{L}. \qquad (2)$$

## 2.4. Downscaling of ET data

Due to the small-scale of vineyard size in the study area (mean size is ~1 ha), many of the 90 m pixels are mixed pixels that incorporate urban and other agricultural effects. To better isolate the vine response, the ET layers were processed using a linear downscaling technique to the same resolution as the original NDVI layers (15 m pixels). A statistical linear downscaling approach was selected as this is simpler and shown to be as effective as a scale-invariant physical modelling approach for downscaling of thermal imagery (Liu and Pu, 2008).

A high quality vector (polygon) landuse map exists for the catchment and was updated for the survey period. The landuse map was simplified into 6 classes (vines, urban, bare soil, cultivated fields, orchards and forest), converted into a 5 m raster layer and intersected with the desired output raster (15 m pixels). Each 15 m pixel therefore contained information on the percentage of each landuse type within the pixel (with nine 5 m pixels per 15 m pixel). The linear downscaling method outlined in Liu and Pu (2008), which uses an initialisation and iteration process, was then adapted and applied to the ET layers.

### 2.4.1. Linear downscaling theory

The ET of the 90 m pixel is a function of the landuse within the pixel. If $f_k(i)$ is the fraction of landuse k within pixel i then the ET of the 90 m pixel is defined as

$$ET(i) = \sum_{k=1}^{K} f_k(i) ET_k(i) \qquad (3)$$

where $ET_k(i)$ is the ET of landuse k in pixel i.

Each 90 m pixel i is composed of 36 15 m subpixels. Each of these subpixels (j) has a unique evapotranspiration ET(ij) that depends on the percentage of each landuse type (derived previously) within the 15 m subpixel. Therefore, ET(ij) can be expressed as:

$$ET(i_j) = \sum_{k=1}^{K} f_k(i_j) ET_k(i) \qquad (4)$$

where $f_k(i_j)$ is the fraction of landuse k in subpixel i_j such that $\sum_{j=1}^{j} ET(i_j) = ET(i)$ for pixel i.

The ET of landuse k in the original pixel i, $ET_k(i)$, can be estimated using an iterative process (Liu and Pu, 2008; Malone et al., 2012).

### 2.4.2. Initiation stage

At the initialization stage, the ET of each subpixel i_j ($\hat{ET}^{l-1}(i_j)$) is set equal to the ET of its corresponding larger pixel i (ET(i)) ($i = 1,...,n$ and $j = 1,...,36$).

### 2.4.3. Iteration stage

At the l-th iteration, a linear regression model between the ET estimate ($\hat{ET}^{l-1}(i_j)$), obtained from the previous iteration ($l-1$), and the land cover fractions ($f_k(i_j)$) is fitted across the 36 subpixels. Note that an intercept term is not used in the regression model in order to avoid co-linearity as the land cover factions sum to one. The model estimate of the ET is obtained as

$$\hat{ET}^{l-1}(i_j) = \hat{ET}_i f_1(i_j) + \cdots + \hat{ET}_k f_k(i_j) \qquad (5)$$

where $\hat{ET}_1, ..., \hat{ET}_k$ are estimates of regression coefficients for the landuse types.

Given the constraint in Eq. (4), $\hat{ET}^{l-1}(i_j)$ is updated to $\hat{ET}^l(i_j)$ using the following equation.

$$\hat{ET}^l(i_j) = \hat{ET}^{l-1}(i_j) + E_T(i) - \frac{1}{36} \sum_{j=1}^{36} \hat{ET}^{l-1}(i_j). \qquad (6)$$

The iteration procedure proceeds until a stop point is reached, in this case when the regression coefficient of determination between two continuous iterations is less than a pre-defined threshold (0.01). At this point, $ET_k$ is returned for each landuse class within the original mixed pixel (i). If landuse k is absent from the original mixed pixel then $ET_k=0$. The important point is that the downscaling extracts the vine (and other landuse) response within the 90 m pixel. Therefore the vine response can be mapped from mixed pixels as well as pure (100% vineyard) pixels. See Liu and Pu (2008) for further details.

## 2.5. Temporal imagery covariates

The process described in §2.4 yielded layers of ET and NDVI in 2007 and 2008 at both 15 and 90 m pixel resolutions. A numeric subscript is used to differentiate between the layers i.e. ET15 and ET90 for ET at 15 and 90 m pixels respectively. The ET15 corresponds to the landuse coefficient derived for vines in §2.4 whilst the ET90 is the 90 m 'mixed' pixel value.

In 2007 a difference was calculated between the early season (May 10th) and late season (September 15th) layers for both covariates at both resolutions. This yielded 4 covariates; $\Delta ET_{15}$, $\Delta ET_{90}$, $\Delta NDVI_{15}$ and $\Delta NDVI_{90}$. These difference covariates provide information on the change in the vines photosynthetically active biomass (NDVI) and in the plant water status/functioning (ET) over the season. Vines with access to subsurface moisture (watertables) should experience relatively little change compared to vines that lack watertable
access at the end of the season (whether or not the watertable is absent or intermittent during the season) (Guix-Hébrard et al., 2007).

For the 2008 time-series, the temporal mean (μ) and standard deviation (σ) of *ET* and *NDVI* were calculated across the 5 seasonal layers at both resolutions (μ*ET*<sub>15</sub>, μ*ET*<sub>90</sub>, μ*NDVI*<sub>15</sub>, μ*NDVI*<sub>90</sub>, σ*ET*<sub>15</sub>, σ*ET*<sub>90</sub>, σ*NDVI*<sub>15</sub> and σ*NDVI*<sub>90</sub>). This provided some general statistics on the *ET* and *NDVI* observation at each pixel, i.e. whether the observation at a location was high or low, variable or constant over time etc. In addition, information on the temporal trend in the *ET* and *NDVI* responses were derived. Since *ET* may be affected by short term climatic variations, the data was scaled before the temporal trend analysis was performed. In this way the temporal change of a pixel was assessed relative to other pixels. For each layer (at both 15 and 90 m), the observations (*ET* or *NDVI*) were scaled within the layer by subtracting the mean of all the observations across the layer from each pixel observation and dividing by the standard deviation of the observations across all layer pixels. For each pixel, the temporal trend in the scaled *ET* and *NDVI* observations was then characterized using linear regression against time (the day of the year when the imagery was collected). The gradient (m) of the linear fit was recorded for each pixel. Similar to the difference covariates (Δ*ET*<sub>15</sub>, Δ*ET*<sub>90</sub>, Δ*NDVI*<sub>15</sub> and Δ*NDVI*<sub>90</sub>) from the 2007 data, these covariates (μ, σ and m) provide information on how the vines are coping as evaporative demand increases. The μ covariates give an overall picture of vine performance (vigorous, well watered or not) and should identify locations with constant access to soil moisture (via very deep soils or watertable access). The σ and m covariates indicate how the vine situation is changing during this crucial phenological period relative to the other vines in the study area. This may help distinguish vines on shallow soils (with no watertable access), which are already exhibiting stress leading into this period (veraison), and deeper soils with good water holding capacity (but no watertable access) that become depleted over this period and reveal a significant change in vegetative vigour and plant water status.

Ideally it would be preferable to have all the imagery-derived covariates in the same year. However, the analysis was constrained by the available ASTER imagery. In 2007, a mid-season time-series was not available and similarly, in 2008, it was not possible to obtain both an early and late season image to compute the difference. Previous work has shown that the vegetative response (pattern) in perennial vineyards is temporally stable in this region (Kazmierski et al., 2011), which permits the mixing of covariates derived from different years.

By convention, the vector of covariates associated with a layer is denoted in bold. Thus all the covariates derived from the *ET*<sub>15</sub> layer (Δ, μ, σ, m) are denoted as *ET*<sub>15</sub>, to distinguish it from the original layer (denoted as *ET*<sub>15</sub>).

## 2.6. Digital terrain attributes

There were two digital elevation models available for the survey area. A larger pixel size DEM captured at 3" (~80 m pixel<sup>1</sup>) by the Shuttle Radar Topographic Mission (SRTM) (obtained from http://srtm.csi.cgiar.org) and a smaller (10 m) pixel size DEM obtained from an aerial LIDAR survey (Bailly et al., 2008). The DEMs were analysed using CLASS Spatial Analyst 2.0.1 (Teng et al., 2005; Teng et al., 2008) to extract primary, secondary and tertiary digital terrain attributes. The CLASS Spatial Analyst software is designed as a hydrological tool to operate on catchment areas. Defining catchment boundaries helps to eliminate some of the boundary effect issues in secondary and tertiary terrain analysis. The programme also calculates Multi-resolution Valley Bottom Flatness (MrVBF) (Gallant and

Dowling, 2003). This has been shown to be related to watertable depths in Australian conditions (Ransley et al., 2007). The CLASS Spatial Analyst was run according to the user guide (Teng et al., 2005) for both DEMs to derive slope (S), aspect (A), curvature (C), planar curvature (PlanC), profile curvature (ProfC), flow derivatives, the compound topographic wetness index (CTI) and MrVBF covariates at both 10 and 80 m resolutions. For the calculation of flow derivatives and the CTI, the D∞ option of the software was used. Recent work has shown that the effectiveness of DTAs depends on the neighbourhood size used in the calculation as well as on the DEM pixel size (Roecker and Thompson, 2010). However, the effect of neighbourhood size is not included in this analysis to avoid too many effects in the analysis. The neighbourhood used is the default settings from within the CLASS Spatial Analyst.

As for the imagery-derived covariates, a numeric subscript is used to differentiate between the same terrain covariate at different resolutions e.g. MrVBF<sub>10</sub> and MrVBF<sub>80</sub>. Again, the vector of covariates associated with a particular DEM are denoted in bold i.e. DEM<sub>10</sub> and DEM<sub>80</sub>.

There were several DTAs generated by the CLASS software that were not included for subsequent analysis. The flow direction and accumulation were omitted as they are strongly correlated to the CTI, which is derived from the flow attributes. Elevation was also omitted. As a covariate, elevation is an absolute value. However, it is the relative location, not the absolute value, of the elevation in the landscape that is important for understanding the pedology and hydrology at a site. The MrVBF DTA is a surrogate for relative elevation. Aspect was also omitted as this does not contribute directly to soil formation. Before omission, modelling was performed with the omitted covariates to verify that they were either not selected or, when selected, they were correlated but not causative and nonsensical in pedological terms in explaining soil or watertable depth.

The schematic diagram in Fig. 3 provides an overview of the steps taken in pre-processing the imagery, deriving the temporal image covariates and generating the digital terrain.

## 2.7. Statistical analysis

The ASTER imagery-derived covariates and DTA-derived covariates were extracted at the 41 survey points. In total there were 14 covariates; 8 imagery-derived covariates (μ*ET*, σ*ET*, μ*NDVI*, σ*NDVI*, m*ET*, m*NDVI*, Δ*ET* and Δ*NDVI*) and 6 DTAs, (S, C, ProfC, PlanC, CTI and MrVBF) at both resolutions. This gave a total of 28 different covariates. To investigate how these covariates related to the two dependent variables (soil depth and water table depth), a binary recursive partition algorithm (regression tree) was used. The partitioning was unsupervised and the point of splitting was determined by selecting the split that maximised the likelihood ratio chi-square statistic (SAS Institute Inc., 2007). Splitting was halted when all individuals had been assigned to a group (node) that was either homogeneous or contained insufficient individuals to form an effective split. The minimum number of individuals in a node was set at 5, which corresponded to 12% of the dataset. Considering a lower number for minimum individuals may produce unrepresentative nodes. Covariates remained within the algorithm at all times, thus the same covariate could be selected at several splitting iterations within the regression tree.

For categorical analyses, the negative loglikelihood (−lnLLH) can be used in a similar manner to the sum of squares for continuous covariates to determine a goodness of fit for a model (denoted as r<sup>2</sup>). This is a measure of the amount of total uncertainty that can be attributed to the model and is calculated as 1 − [(−lnLLH Model)/ (−lnLLH for Corrected Total)] (SAS Institute Inc., 2007). The −lnLLH for Corrected Total is derived by fitting a model where the probabilities are estimated by fixed rates for each response level. This provides the background uncertainty for a model with no effects and is analogous to the Total Sum of Squares (Model SS + Residual SS)

<footnote n="1">¹ In general 3" relates to a 90 m pixel, however, at the latitude of this study the resulting image obtained from CGIAR (SRTM_37_4) had a pixel size of 80.1 m.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Maps of the derived imagery-based covariates. Covariates at different pixel sizes are paired and presented on a common legend. Labels are provided on the legend. In each pairing the smaller pixel layer (15 m) is on the left and larger pixel layer (90 m) is on the right. Gaps in the *ET*₁₅ layers are due to areas where there were no vines according to the landuse map i.e. the vine coefficient after downscaling = 0. The white area in the lower centre of all layers masks a cloud artefact.</figcaption></figure>

in ANOVA with continuous data. A leave-one-out cross validation (LOOCV) was performed for each model. The LOOCV was also assessed using an $r^2$ value based on the —InLLHs. All analysis was performed using the Partition platform in JMP 8.0.1 (SAS Institute Inc., Cary, NC, USA).

## 2.8. Models

The intention of the analysis is to investigate the utility of both ASTER imagery-derived covariates and DTA-derived covariates in modelling soil and watertable depth, and to investigate if there is an effect of pixel size (associated with the different pixel sizes of the

ASTER imagery bands) in the modelling. To achieve this, the following combinations of independent covariates were used as inputs to the regression tree algorithm.

Model 1 Only terrain covariates at 10 m pixel size.

Model 2 Only terrain covariates at 80 m pixel size.

Model 3 Only image covariates at 15 m pixel size.

Model 4 Only image covariates at 90 m pixel size.

Model 5 Terrain and Image covariates at the smallest pixel size (10/15 m).

Model 6 Terrain and Image covariates at the largest pixel size (80/90 m).

Model 7 All available covariates.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Maps of the DTAs used as covariates. Covariates at different pixel sizes are paired and presented on a common legend. Labels are provided on the legend. In each pairing the smaller pixel layer (10 m) is on the left and the larger pixel layer (80 m) is on the right. The CTI and MrVBF are tertiary terrain attributes and were derived after the DEM was constrained to the watershed to avoid boundary effects. The white outline indicates the survey area.</figcaption></figure>

For the discussion, the subscript *s* or *w* is appended to the model number to specify when the model being discussed relates to a specific dependent variable (*s*—soil depth and *w*—watertable depth).

# 3. Results and discussion

## 3.1. Effect of pixel size on the imagery and terrain covariate maps

Maps of the various imagery-derived covariates and DTAs are presented in Figs. 1 and 2 respectively. Each covariate is presented at both a small and large pixel size on a common legend for comparison. As expected, the **NDVI**<sub>15</sub> maps have more fine detail than the **NDVI**<sub>90</sub> maps. The filtered **ET**<sub>15</sub> maps do not show more detail than the **ET**<sub>90</sub> maps but do show a spatial smoothing of the response due to the down-scaling. The terrain maps also show similar patterns between **DEM**<sub>10</sub> and **DEM**<sub>80</sub> covariates with much greater detail in the **DEM**<sub>10</sub> covariates e.g. *Curve*<sub>10</sub> vs. *Curve*<sub>80</sub>. The exception is the tertiary CTI covariate. At a large pixel size, the CTI<sub>80</sub> follows the general stream flow pattern in the catchment. However the lower connectivity in the CTI<sub>10</sub> map means that the stream order is not as clear.

## 3.2. Overall modelling results

A summary of the results of running Models 1–7 on the soil and watertable depth responses is given in Table 1. Fits ($r^2$) are given

for both the overall model (fitted with all the data and no validation) and the LOOCV. Also recorded are the number of splits in the regression tree and the covariates selected by the algorithm. In general the selected covariates are arranged in the order selected by the partitioning algorithm; however repetitions of covariate selections are not indicated. Given that a covariate could be selected at multiple splits in the regression tree, the number of selected covariates does not always equal the number of splits. For example, Model 1<sub>s</sub> of soil depth has 5 splits but only 3 covariates as *Curve*<sub>10</sub> and *Slope*<sub>10</sub> are repeated.

Model fitness was assessed primarily from the LOOCV fits. For soil depth, the best result was obtained from Model 7<sub>s</sub>, which included all available covariates at all scales. The models that used only DTAs (Models 1 and 2) or only imagery (Models 3 and 4) at a single scale generally performed poorly compared to models that incorporated both DTAs and imagery-derived covariates.

For the watertable depth modelling, the best result was obtained from Models 5<sub>w</sub> and 7<sub>w</sub>, which were identical and used the small pixel imagery-derived covariates and small pixel DTAs. Using only imagery-derived covariates, at either resolution (Models 3<sub>w</sub> and 4<sub>w</sub>), produced the worst results. The larger pixel-size **DTA**<sub>80</sub> covariates (Model 2<sub>w</sub>) also performed poorly. The best models all utilised the MrVBF<sub>10</sub> covariate in the primary split (Models 1<sub>w</sub>, 5<sub>w</sub> and 7<sub>w</sub>).

The best regression tree models were able to explain more of the uncertainty in soil depth (LOOCV $r^2$=0.89) than watertable
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Schematic illustrating (A) the pre-processing of the ASTER imagery into evapotranspiration (ET) and normalised difference vegetative index (NDVI) layers at two different pixel sizes (15 and 90 m); (B) the generation of the temporal imagery covariates from the ET and NDVI layers; and (C) the generation of primary, secondary and tertiary digital terrain attributes using the CLASS Spatial analyst. Bold black outlined text indicates input layers, dotted grey circles indicate intermediate steps and grey outlined text indicates derived covariate layers.</figcaption></figure>

(LOOCV r² = 0.63) (Fig. 4). The best model for watertable depth (Models 5<sub>w</sub> and 7<sub>w</sub>) was more parsimonious (3 splits) than the best soil depth model (Model 7<sub>s</sub>, 5 splits). In both cases, the best models incorporate both small pixel DTAs and small pixel temporal imagery-derived covariates in the partitioning process.

### 3.3. Selection of DTAs and effect of pixel size (Models 1–2)

The results indicate that decreasing the DEM pixel size from 90 m to 10 m pixels produced more effective DTAs for modelling both soil depth and watertable depth. By decreasing the pixel size there was
<table><caption>**Table 1**

Fits from the overall model and a leave-one-out cross validation (LOOCV) analysis and the number of partitions for a series of regression tree models of soil and watertable depth that use imagery-derived covariates or terrain covariates (DTAs) or a mixture of both imagery and terrain covariates.</caption><thead><tr><th>Dependent covariate</th><th>Model<sup>a</sup></th><th>Model r<sup>2</sup></th><th>LOOCV r<sup>2</sup></th><th>No. splits</th><th>Selected independent covariates</th></tr></thead><tbody><tr><td>Soil</td><td>Model 1</td><td>0.58</td><td>0.46</td><td>5</td><td>Curve<sub>10</sub>, MrVBF<sub>10</sub>, Slope<sub>10</sub></td></tr><tr><td>Soil</td><td>Model 2</td><td>0.50</td><td>0.32</td><td>7</td><td>CTI<sub>80</sub>, ProfC<sub>80</sub>, Slope<sub>80</sub></td></tr><tr><td>Soil</td><td>Model 3</td><td>0.58</td><td>0.45</td><td>5</td><td>σET<sub>15</sub>, σNDVI<sub>15</sub>, μNDVI<sub>15</sub>, μET<sub>15</sub></td></tr><tr><td>Soil</td><td>Model 4</td><td>0.51</td><td>0.41</td><td>4</td><td>mET<sub>90</sub>, μET<sub>90</sub>, μNDVI<sub>90</sub></td></tr><tr><td>Soil</td><td>Model 5</td><td>0.75</td><td>0.64</td><td>5</td><td>σET<sub>15</sub>, PlanC<sub>10</sub>, MrVBF<sub>10</sub>, mNDVI<sub>15</sub></td></tr><tr><td>Soil</td><td>Model 6</td><td>0.51</td><td>0.41</td><td>4</td><td>mET<sub>90</sub>, μET<sub>90</sub>, μNDVI<sub>90</sub></td></tr><tr><td>Soil</td><td>Model 7</td><td>0.91</td><td>0.89</td><td>5</td><td>σET<sub>15</sub>, PlanC<sub>10</sub>, Curve<sub>80</sub>, MrVBF<sub>10</sub>, mNDVI<sub>15</sub></td></tr><tr><td>Watertable</td><td>Model 1</td><td>0.66</td><td>0.57</td><td>4</td><td>MrVBF<sub>10</sub>, PlanC<sub>10</sub>, Slope<sub>10</sub></td></tr><tr><td>Watertable</td><td>Model 2</td><td>0.58</td><td>0.46</td><td>5</td><td>CTI<sub>80</sub>, PlanC<sub>80</sub></td></tr><tr><td>Watertable</td><td>Model 3</td><td>0.45</td><td>0.28</td><td>5</td><td>μNDVI<sub>15</sub>, μET<sub>15</sub>, mNDVI<sub>15</sub>, σNDVI<sub>15</sub></td></tr><tr><td>Watertable</td><td>Model 4</td><td>0.42</td><td>0.28</td><td>5</td><td>σET<sub>90</sub>, μET<sub>90</sub>, μNDVI<sub>90</sub>, mET<sub>90</sub></td></tr><tr><td>Watertable</td><td>Model 5</td><td>0.71</td><td>0.61</td><td>3</td><td>MrVBF<sub>10</sub>, μET<sub>15</sub>, Slope<sub>10</sub></td></tr><tr><td>Watertable</td><td>Model 6</td><td>0.69</td><td>0.59</td><td>5</td><td>CTI<sub>80</sub>, mNDVI<sub>90</sub>, Slope<sub>80</sub>, μNDVI<sub>90</sub>, μET<sub>90</sub></td></tr><tr><td>Watertable</td><td>Model 7</td><td>0.71</td><td>0.61</td><td>3</td><td>MrVBF<sub>10</sub>, μET<sub>15</sub>, Slope<sub>10</sub></td></tr></tbody></table>

a corresponding increase of 11–14% in the amount of uncertainty (r²) explained in both the soil depth and watertable depth modelling. The DTAs at a given pixel size were more effective at partitioning the uncertainty in watertable depth than soil depth.

The MrVBF₁₀ covariate was the most effective covariate for the first partitioning between shallow and deep watertables. However, the equivalent larger pixel covariate (MrVBF₈₀) was not selected in the modelling (Models 2<sub>w</sub> and 6<sub>w</sub>) despite similarities in the two maps (Fig. 2). With the 80 m DEM, the CTI and curvature covariates were preferred to MrVBF. The Curve₈₀ covariate was also selected in the best soil depth model (Model 7<sub>s</sub>) over Curve₁₀, however the Curve₁₀ and associated PlanC₁₀ and ProfC₁₀ covariates where often present in soil depth models that incorporated the DEM₁₀ covariates (for example Models 5<sub>s</sub> and 7<sub>s</sub>). It is possible that the preference for the Curve₈₀ covariate over Curve₁₀ is due to artefacts arising from the choice of neighbourhood size in the calculation (Roecker and

Thompson, 2010). Further studies in this area should investigate these possible effects.

### 3.4. Selection of temporal image covariates and effect of pixel size (Model 3–4)

There was little difference between the fit of Models 3 and 4 to soil depth and no difference in their fits to watertable depth. The pixel size had no effect on the amount of uncertainty explained, although it did affect the covariates selected. When modelling soil depth, ET covariates were selected as the first partition for both models although the type of covariates chosen differed (σ for Model 3<sub>s</sub> and m for Model 4<sub>s</sub>).

The imagery-derived covariates explained more of the uncertainty in soil depth, at both pixel sizes, than they did for watertable depth. This is the opposite result to the DTA models (Models 1 and 2). It

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Regression trees of the best performed models for both dependent variables — Model 7<sub>s</sub> for soil depth (left) and the identical Models 5<sub>w</sub> and 7<sub>w</sub> for watertable depth (right). Count indicates the total number of data at each node whilst the bar provides an indication of the proportion of shallow (grey) and deep (black) data at each node. The covariate and the value used for each split are also shown.</figcaption></figure>

<footnote>a Models are specified in §2.8.</footnote>
appears that information on vegetation development is not a good indicator of depth to the watertable; rather it indicates access to the watertable. The imagery is unable to differentiate between conditions where vines are accessing the watertable due to a shallow watertable or conditions where access is achieved by a deep watertable and deep rooting depth.

### 3.5. Effect of modelling with both DTAs and imagery covariates (Models 5–6)

When the DTAs and temporal imagery-derived covariates were mixed for a given pixel size (Models 5 and 6), the gain compared to the best results from the previous models (Models 1–4) differed for the two dependent variables. In the case of soil depth, Model 5s, which used the small pixel covariates, explained 18% more of the uncertainty in the soil depth data than Model 1s (the best result from Models 1–4). Model 6s was identical to Model 4s indicating that at a 90 m pixel resolution the DTAs provide no additional value to the imagery data.

The watertable depth models that used both imagery-derived covariates and DTAs produced similar results regardless of the pixel size ($r^2$ of 0.61 and 0.59 for Models 5w and 6w respectively) and these were only slightly better than Model 1w ($r^2$ = 0.57). Both Models 5w and 6w used a mixture of imagery-derived covariates and DTAs. Model 5w was the most parsimonious with 3 splits compared to 5 splits for Model 6w and 4 splits for Model 1w.

### 3.6. Multi-scale covariate modelling. (Model 7)

For soil depth, the Model 7s fit was superior to all other models ($LOOCV \ r^2$ = 0.89) and explained at least 25% more of the uncertainty than the next best model (Model 5s). The difference between Model 7s and 5s is the inclusion of the Curve80 covariate in Model 7s. This larger pixel indicator of curvature fills a gap that the small pixel imagery and DTAs cannot in Model 5s. The results from the regression tree analysis indicate that multi-scale DTAs and imagery-derived covariates are complementary, and a mixed approach yields a superior non-linear model. The choice of multi-scale DTAs in the best performed soil depth model would seem to support previous observations (Roecker and Thompson, 2010; Smith et al., 2006) that multi-scale terrain characterization is preferable for modelling and prediction.

In the soil depth models, the $\sigma ET_{15}$ and $mNDVI_{15}$ covariates were selected in the best models that mixed the terrain and imagery-derived covariates (Models 5s and 7s) rather than the mean response ($\mu ET$ and $\mu NDVI$), which were selected in the poorer performing Model 6s. The $\sigma ET_{15}$ and $mNDVI_{15}$ covariates are describing the change in the vines vegetative response relative to other vines in the image and are able to identify vines with either poor or good water supply as a result of either soil depth or watertable depth or both. These imagery observations may be independent of the underlying geo-hydrology and are able to identify areas where a vines access to groundwater does not correspond to the general rules of pedogenesis and geo-morphology, for example a shallow soil over a shallow watertable in a relatively low-lying position in the landscape. Since this is a complex pedo-landscape, there are several areas where soil depth is independent of watertable (Guillaume Coulouma, UMR LISAH, INRA, Montpellier, pers. comm.). This complexity is hypothesised to be the reason that the vine response provides little value as a covariate for modelling watertable depth.

For watertable depth modelling there was very little advantage in using multi-scale covariates. Model 7w produced similar results to Models 5w and 1w. Of primary interest in Model 7w for watertable depth is the choice of the scaled $\mu ET_{15}$ covariate (discussed further below). Overall the results for modelling watertable depth were not surprising. With the exception of perched watertables, the depth function for watertable is driven primarily by landscape and

geohydrological conditions with fluctuations imposed by the flux of groundwater in the system (Guix, 2005). These processes are generally well understood and fairly universal in application i.e. they are not specific to particular regions (Tassinari, 1998). For this reason the alternative DTAs, such as MrVBF, also appear to work well in this area even though the pedo-landscape is very different from the pedo-landscapes where the attribute was developed and has been previously applied (Ransley et al., 2007). Information on vegetative growth is neither superior nor complementary to information on terrain factors for modelling watertable depth.

In general, the temporal imagery-derived covariates were identified as relevant discriminators although generally the small pixel covariates were preferred (ET<sub>15</sub>, NDVI<sub>15</sub>). A high $\sigma ET_{15}$ value clearly defines points where shallow soils exist and is the first split in all models where the ET<sub>15</sub> covariates are available. In the best watertable depth model (Model 7w, Fig. 4b) a low $\mu ET_{15}$ identifies points with a deep watertable (poor vine access to ground water) in areas where shallow watertables are expected (high MrVBF). Similarly, in the last split of the best soil depth model (Fig. 4a), the $mNDVI_{15}$ response is differentiating soil depth on points that are not located on the valley floor. In these areas (slopes and hill crests), locations with a deep soil profile can be determined by a positive trend in the NDVI response ($mNDVI_{15}$ > 0.008) around veraison (relative to other locations in the catchment).

The relevance of the ET<sub>15</sub> covariates indicates an advantage to the downscaling of the ET data from 90 to 15 m pixels. The original larger pixel imagery (ET<sub>90</sub>) is known to contain 'mixed' pixels. Soil sampling was only performed within vineyards; however the response over the 90 m pixels is often 'mixed' (multiple landuses). The downscaling approach applied here attempted to filter the vineyard response, and the results indicate that it has produced more effective covariates for modelling.

### 3.7. Implications for the GlobalSoilMap.net project

The models performed here indicate that using covariates at multiple scales (Model 7) confers benefits to the soil depth modelling process. However multi-scale covariates complicate the movement from a digital soil modelling to a digital soil mapping application. McBratney (1998) provides an introduction to some concepts for model development in multi-scale situations. McBratney (1998) suggests that either a) the inputs be scaled (either up or down) to the relevant model/output scale before processing, b) the model be run at the preferred input scale (where some inputs may need to be scaled) and then the output rescaled to the relevant support or c) the inputs are added at their original scale and scaling is performed within the model to produce an output on a different support. The results here would indicate that some version of option (c) is preferable, especially for soil depth. In the context of the GlobalSoilMap.net project, multi-scale inputs are likely and from the results here preferable, however only a single scale output (90 m pixel size) is desired. Digital soil depth models therefore need to be able to both predict and scale the input data. Currently, most approaches to DSM have followed approach (a) in scaling the inputs to the desired support (usually 90 m pixels) prior to modelling/mapping.

Whilst filtering the imagery is beneficial, it cannot be ascertained from this study whether or not the pixel sizes used here are the optimum. It may be preferable to upscale/downscale to an intermediate value, such as 30, 45 or 60 m pixels. This is certainly an area for further investigation. The availability of either a 1" SRTM (~30 m) DEM for some regions of the world or a downscaled 30 m grid for others (e.g. the 30 m DEM of Australia being generated by the Australian Government (http://www.ga.gov.au/topographic-mapping/digital-elevation-data.html)) may present opportunities at this pixel size. It is certainly of interest to understand how the terrain covariates behave at these intermediate pixel sizes.
The approach of using vegetative covariates takes advantage of (and is reliant on) the canopy development within a perennial rain-fed monoculture system to reflect subsoil properties. Therefore, the application of this approach will likely be most effective in other regions where a perennial rain-fed monoculture exists. This may be an agricultural (e.g. other viticultural or horticultural zones) or a natural environment (e.g. coniferous forests). However, the potential application to annual (monoculture) systems should not be dismissed. In many cases it is in these areas where high quality digital soil maps are particularly needed to optimise production potential. Provided that the variance associated with management (sowing date, fertilisation and crop development) can be minimised, the canopy response across an annual system may reflect subsoil conditions. This may reflect a subsoil constraint(s), i.e. the effective rooting depth, rather than the actual soil depth. At a farm level, inter-annual temporal *NDVI* layers have been used to map subsoil constraints in annual wheat (*Triticum aestivum*) systems (Dang et al., 2011). The challenge remains to upscale this to a regional scale that is compatible with regional DSMM.

Finally, it is important to recall that regression trees and partitioning algorithms are known to have limitations when used for prediction. The intention here is to use them for data mining, not for prediction. A partitioning approach to a prediction may or may not be the preferred method and the determination of this is outside the scope of this paper. When partitioning algorithms are incorporated into prediction methodologies then there is usually some form of boot-strapping or jack-knifing, for example the random forest algorithm (Breiman, 2001), is preferred.

## 4. Conclusions

The regression tree analysis showed some clear trends for modelling soil depth and water table depth using remotely-sensed vegetative and digital terrain attributes. The incorporation of vegetative covariates derived from the thermal wavebands of satellite imagery was useful in modelling soil depth but not watertable depth. Terrain attributes from a small pixel DEM proved to be the best covariates for modelling both soil and watertable depths. There was a distinct advantage in modelling soil depth using both the vegetative and terrain covariates at multiple scales compared to single scale models. The value of the vegetative covariates in modelling soil depth is based on the ability of the plant to reflect the soil conditions in which it is growing. This illustrates that some soil attributes are able to be modelled through a vegetation-soil inference system. The success of the multi-scale model in this exploratory analysis indicates that further work to create multi-scale predictive models for DSMM is likely to be beneficial for some soil attributes.

## Acknowledgements

James Taylor's position with INRA is funded through the Agropolis Foundation and their financial support is greatly appreciated. ASTER data was kindly provided by the Jet Propulsion Laboratory. The obtaining of ASTER imagery derived covariates was possible; thanks to grants from the Chilean National Commission for Scientific and Technological Research, the French National Institute for Agricultural Research and the French National Program for Remote Sensing.

## References

* Baily, J., Lagacherie, P., Millier, C., Puech, C., Kosuth, P., 2008. Agrarian landscapes linear features detection from LIDAR: application to artificial drainage networks. International Journal of Remote Sensing 29 (12), 3489–3508.

* Boulet, G., Mougenot, B., Ben Abdelouahab, T., 2009. An evaporation test based on thermal infra-red remote-sensing to select appropriate soil hydraulic properties. Journal of Hydrology 376 (3–4), 589–598.

* Bramley, R.G.V., Proffitt, A.P.B., Corner, R.J., Evans, T.D., 2000. Variation in grape yield and soil depth in two contrasting Australian vineyards. In: Adams, J.A., Metherell, A.K. (Eds.), Soil 2000: New Horizons for a New Century. Australian and New Zealand Second Joint Soils Conference. Volume 2: Oral papers. Lincoln University, New Zealand Society of Soil Science, pp. 29–30.

* Breiman, L., 2001. Random forests. Machine Learning 45 (1), 5–32.

* Carre, F., Girard, M.C., 2002. Quantitative mapping of soil types based on regression kriging of taxonomic distances with landform and land cover attributes. Geoderma 110 (3–4), 241–263. http://dx.doi.org/10.1016/S0016-7061(02)00233-1.

* Crow, W.T., Kustas, W.P., Prueger, J.H., 2008. Monitoring root-zone soil moisture through the assimilation of a thermal remote sensing-based soil moisture proxy into a water balance model. Remote Sensing of Environment 112 (4), 1268–1281.

* Dang, Y.P., Pringle, M.J., Schmidt, M., Dalal, R.C., Apan, A., 2011. Identifying the spatial variability of soil constraints using multi-year remote sensing. Field Crops Research 123 (3), 248–258.

* Finke, P.A., Brus, D.J., Bierkens, M.F.P., Hoogland, T., Knotters, M., de Vries, F., 2004. Mapping groundwater dynamics using multiple sources of exhaustive high resolution data. Geoderma 123, 23–39.

* Gallant, J.C., Dowling, T.L., 2003. A multiresolution index of valley bottom flatness for mapping depositional areas. Water Resources Research 39 (12), 1347–1359.

* Galleguillos, M., Jacob, F., Prévot, L., Lagacherie, P., Liang, S., 2011a. Spatializing vineyard hydric status within heterogeneous Mediterranean watershed from high spatial resolution optical remote sensing. IEEE Geoscience and Remote Sensing Letters 8 (1), 168–172.

* Galleguillos, M., Jacob, F., Prévot, L., French, A., Lagacherie, P., 2011b. Comparison of two temperature differencing methods for estimating daily evapotranspiration with ASTER data over a Mediterranean watershed. Remote Sensing of Environment 115 (6), 1326–1340.

* Goosens, R., Alavi Panah, S.K., De Dapper, M., Kissyar, O., 1999. The use of thermal band of Landsat TM for the study of soil salinity in Iran (Ardakan area) and Egypt (Ismailia Province). Proc. Int. Conf. on Geoinformatics for Natural Resource Assessment, Monitoring and Management, Indian Institute of Remote Sensing, NRSA, Dehradun, India, pp. 454–459.

* Grunwald, S., 2009. Multi-criteria characterization of recent digital soil mapping and modeling approaches. Geoderma 152, 195–207.

* Guix, N., 2005. Prise en compte des nappes superficielles pour spatialiser le bilan hydrique à l'échelle sub-régionale Cas de la vigne en basse léeve de la Peyne - Hérault, France. Mémoire de thèse de Doctorat de l'université de Montpellier 2.

* Guix-Hébrard, N., Voltz, M., Trambouze, W., Garnier, F., Gaudillère, J.P., Lagacherie, P., 2007. Influence of watertable depths on the variation of grapevine water status at the landscape scale. European Journal of Agronomy 27, 187–196.

* Kazmierski, M., Glemas, P., Rousseau, J., Tisseyre, B., 2011. Temporal stability of within-field patterns of ndvi in non irrigated mediterranean vineyards. Journal International des Sciences de la Vigne et du Vin 45 (2), 61–73.

* Liu, D., Pu, R., 2008. Downscaling thermal infrared radiance for subpixel land surface temperature retrieval. Sensors 8, 2695–2706.

* Liu, J., Pattey, E., Nolin, M.C., Miller, J.R., Ka, O., 2008. Mapping within-field soil drainage using remote sensing, DEM and apparent soil electrical conductivity. Geoderma 143 (3–4), 261–272. http://dx.doi.org/10.1016/j.geoderma.2007.11.011.

* Malone, B.P., McBratney, A.B., Minasny, B., Wheeler, I., 2012. A general method for downscaling earth resource information. Computational Geosciences 41, 119–125. http://dx.doi.org/10.1016/j.cageo.2011.08.021.

* Martínez-Montoya, J.F., Herrero, J., Casterad, M.A., 2010. Mapping categories of gypseous lands in Mexico and Spain using Landsat imagery. Journal of Arid Environments 74 (8), 978–986. http://dx.doi.org/10.1016/j.jaridenv.2010.01.011.

* McBratney, A.B., 1998. Some considerations on methods for spatially aggregating and disaggregating soil information. Nutrient Cycling in Agroecosystems 50 (1–3), 51–62.

* McBratney, A.B., Mendonça Santos, M.L., Minasny, B., 2003. On digital soil mapping. Geoderma 117 (1–2), 3–52. http://dx.doi.org/10.1016/S0016-7061(03)00223-4.

* Meirik, E., Frazier, B., Brown, D., Roberts, P., Rupp, R., 2010. ASTER-based vegetation map to improve soil modeling in remote areas. In: Boettinger, J.L., Howell, D.W., Moore, A.C., Hartemink, A.E., Kienast-Brown, S. (Eds.), Digital Soil Mapping, Progress in Soil Science, Volume 2, Part II, Section A. Springer, Berlin, pp. 113–122. http://dx.doi.org/10.1007/978-90-481-8863-5_10.

* Mulder, V.L., de Bruin, S., Schaepman, M.E., Mayr, T.R., 2011. The use of remote sensing in soil and terrain mapping — a review. Geoderma 162 (1–2), 1–19. http://dx.doi.org/10.1016/j.geoderma.2010.12.018.

* Olioso, A., Inoue, Y., Ortega-Farias, S., Demarty, J., Wigneron, J.-P., Braud, I., Jacob, F., Lecharpentier, P., Ottlé, C., Calvet, J.-C., Brisson, N., 2005. Future directions for advanced evapotranspiration modelling: assimilation of remote sensing data into crop simulation models and SVAT models. Irrigation and Drainage Systems 19 (3–4), 355–376.

* Ransley, T., Tottenham, R., Baskaran, S., Brodie, R., 2007. Development of Method to Map Potential Stream-Aquifer Connectivity: A Case Study in the Border Rivers Catchment. Bureau of Rural Sciences, Canberra, 27 pp.

* Rivero, R.G., Grunwald, S., Bruland, G.L., 2007. Incorporation of spectral data into multivariate geostatistical models to map soil phosphorus variability in a Florida wetland. Geoderma 140 (4), 428–443. http://dx.doi.org/10.1016/j.geoderma.2007.04.026.

* Roecker, S.M., Thompson, J.A., 2010. Scale effects on terrain attribute calculation and their use as environmental covariates for digital soil mapping. In: Boettinger, J.L., Howell, D.W., Moore, A.C., Hartemink, A.E., Kienast-Brown, S. (Eds.), Digital Soil Mapping, Progress in Soil Science, Part II, Section A. Springer, Berlin, pp. 55–66. http://dx.doi.org/10.1007/978-90-481-8863-5_5.

* Roerink, G.J., Su, Z., Menenti, M., 2000. S-SEBI: a simple remote sensing algorithm to estimate the surface energy balance. Physics and Chemistry of the Earth-Part B: Hydrology, Oceans and Atmosphere 25 (2), 147–157.

* Rouse, J.W., Haas, R.H., Schell, J.A., Deering, D.W., 1973. Monitoring vegetation systems in the Great Plains with ERTS. Proceedings of the Third ERTS Symposium, Paper A 20, NASA SP-351, pp. 309–317.
SAS Institute Inc., 2007. JMP Statistics and Graphics Guide, Release 7. SAS Institute Inc.,
Cary, NC, USA. ISBN: 978-1-59994-408-1.

Smith, M.P., Zhu, A.X., Burt, J.E., Stiles, C., 2006. The effects of DEM resolution and
neighbourhood size on digital soil survey. Geoderma 137, 58–69.

Tassinari, C., 1998. Contribution 'a l'étude des sols méditerranéens (Languedoc, France): caractérisation de leurs nappes et des indicateurs morphologiques associés. Thèse de Doctorat de l'Ecole Nationale Supérieure Agronomique de Montpellier, 155p.

Taylor, J.A., Acevedo-Opazo, C., Ojeda, H., Tisseyre, B., 2010. Identification and significance of sources of spatial variation in grapevine water status. Australian Journal of Grape and Wine Research 16 (1), 218–226. http://dx.doi.org/10.1111/j.1755-0238.2009.00066.x.

Teng, J., Vaze, N.K., Tuteja, N.K., Gallant, J., 2005. CLASS spatial analyst: A GIS based tool for distributed hydrological modelling. In: Zerger, A., Argent, R.M. (Eds.), MODSIM 2005; International Congress on Modelling and Simulation. Modelling and Simulation Society of Australia and New Zealand. ISBN: 0-9758400-2-9, pp. 1485–1491.

Teng, J., Vaze, J., Tuteja, N.K., Gallant, J.C., 2008. A GIS-based tool for spatial and distributed hydrological modelling: CLASS spatial analyst. Transactions in GIS 12 (2), 209–225.

Tucker, C.J., 1979. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment 8, 127–150.