Disaggregating a regional-extent digital soil map using Bayesian area-to-point regression kriging for farm-scale soil carbon assessment

Most soil management activities are implemented at farm scale, yet digital soil maps are commonly available at regional or national scale. Disaggregating these regional and/or national maps is applicable for farm-scale tasks, particularly in data-poor or limited situations. Although disaggregation is a frequently discussed topic in recent digital soil mapping literature, the uncertainty of the disaggregation process is not often discussed. Underestimation of inferential or predictive uncertainty in statistical modelling leads to inaccurate statistical summaries and overconfident decisions. The use of Bayesian inference allows for quantifying the uncertainty associated with the disaggregation process. In this study, a framework of Bayesian area-to-point regression kriging (ATPRK) is proposed for downscaling soil attributes, in particular, maps of soil organic carbon. An estimation of point support variograms from block-supported data was carried out using the Monte Carlo integration via the Metropolis–Hastings algorithm. A regional soil carbon map with a resolution of 100 m (block support) was disaggregated to 10 m (point support) information for a farm in northern New South Wales (NSW), Australia. The derived point support variogram has a higher partial sill and nugget, while the range and parameters do not deviate much from the block support data. The disaggregated fine-scale map (point support with a grid spacing of 10 m) using Bayesian ATPRK had an 87 % concordance correlation with the original coarse-scale map. The uncertainty estimates of the disaggregation process were given by a 95 % confidence interval (CI) limit. Narrow CI limits indicate that the disaggregation process gives a fair approximation of the mean soil organic carbon (SOC) content of the study site. The Bayesian ATPRK approach was compared with dissever, which is a regression-based disaggregation algorithm. The disaggregated maps generated by dissever had 96 % concordance correlation with the coarse-scale map. Dissever achieves this higher concordance correlation through an iteration process, while Bayesian ATPRK is a one-step process. The two disaggregated products were validated with 127 independent topsoil carbon observations. The validation concordance correlation coefficient for Bayesian ATPRK disaggregation was 23 %, while downscaled maps generated from dissever had 18 % concordance correlation coefficient (CCC). The advantages and limitations of both disaggregation algorithms are discussed.

1 Introduction

The proliferation of digital soil mapping (DSM) and modelling techniques during the past decade has made vast amounts of digital soil information available (Lagacherie, 2008; Grunwald, 2009; Minasny and McBratney, 2016). For example, the GlobalSoilMap project (Arrouays et al., 2014) aspired towards delivering a global coverage of soil attributes at 100 m resolution for all six standard soil depths. The Soil and Landscape Grid of Australia (SLGA; Grundy et al., 2015), the maps of topsoil organic carbon content across Europe (de Brogniez et al., 2014), the national digital soil maps of France (Mulder et al., 2016), and the carbon map of Nigeria (Akpa et al., 2016) are examples of large-extent digital soil mapping. However, these national or regional digital maps do not capture the sufficient spatial variability of soil for farm-scale decision-making (Malone et al., 2017). Yet, these maps can be converted into useful farm-scale information using a combination of fine-scale environmental covariates and geostatistical methods. This process is called downscaling or disaggregation (Malone et al., 2017).

Disaggregation is the process of detailing information collected at a small or coarse scale towards a larger or finer scale (Cheng, 2008). Disaggregating techniques can be classified into two categories, namely geostatistical and regression techniques. Area-to-point kriging (ATPK; Kyriakidis, 2004) and its variants, such as area-to-point regression kriging (ATPRK), are the most common geostatistical techniques applied to disaggregate soil resource information. Incorporating both point and areal data in the kriging system is a variant of ATPK presented by Goovaerts (2011). Firstly, that study demonstrated how ATPK mapped the variability within soil units and, secondly, introduced a general formulation for incorporating area-to-area, area-to-point, and point-to-point covariance in the kriging system. Replacing the areal mean with summary statistics in the kriging system is another improved ATPK technique suggested by Brus et al. (2014). Román Dobarco et al. (2016) extended this approach using summary statistics in the ATPRK system. Furthermore, a Bayesian regression kriging model was presented by Reid et al. (2013). They introduced the Gaussian process model as a novel image-fusion technique for the panchromatic sharpening of low-resolution satellite images. They used images from unmanned aerial vehicles flown over agricultural fields to demonstrate their model.

Regression approaches are also a popular approach for downscaling soil attributes. For example, Malone et al. (2012), formalised a statistical downscaling algorithm. In addition, dissever is based on the Liu and Pu (2008) algorithm for downscaling continuous earth resources with a case study of downscaling soil organic carbon content. This method was refined by Roudier et al. (2017), allowing the user to choose between a range of regression models apart from the originally contained generalised additive models (GAMs). The Disaggregation and Harmonisation of Soil Map Units through Resampled Classification Trees (DSMART) approach is another popular regression-based approach proposed by Odgers et al. (2014) to downscale soil map units to soil series information. However, it is different to the discussed method since there is no fine gridding involved.

Disaggregation is a change in the support problem (Cressie, 1991; Gotway and Young, 2002). Spatial support is the geometrical size, shape orientation, and spatial unit of an observation or a prediction. In the disaggregation process, when the support changes the associated statistical and spatial properties of the variable also change, adding a major source of uncertainty (Truong et al., 2014). ATPK is a theoretically based method that uses the classical kriging principles to predict point support using block support information (Kyriakidis, 2004). Incorporating fine-scale covariate data in the kriging system (ATPRK) is superior to ATPK as it allows for within-class variability, spatial autocorrelation, and also ancillary covariates representing the Soils, Climate, Organisms, Parent material, Age and (N) space or spatial position (SCORPAN) factors (McBratney et al., 2003; Kerry et al., 2012).

In disaggregating soil information, quantifying the uncertainties of this process is not often discussed. Underestimation of inferential or predictive uncertainty in statistical modelling can lead to inaccurate statistical summaries and overconfident decisions (Draper, 1995). Therefore, it is important to account for the underlying parameter uncertainty associated with the disaggregation process. The use of Bayesian inference in the ATPRK system enables one to quantify the uncertainty associated with inferring point support parameters from the block-supported system, in addition to estimating the product uncertainty (Diggle et al., 1998).

To the best of our knowledge, Bayesian ATPRK has not been used for downscaling soil information. In this paper, we present a Bayesian ATPRK approach for downscaling soil attributes. Specifically, the study focuses on downscaling regional soil organic carbon (SOC) maps to high-resolution farm-scale maps. SOC is arguably a key soil property due to its conferring benefits to the soil's physical and chemical properties and its potential for storing atmospheric carbon. Accurate estimations of SOC stocks at farm scale are important for precision agricultural needs. In addition, the fine-scale maps can be used for designing sampling strategies for the quantification of SOC stocks (de Gruijter et al., 2016) required by carbon sequestration auditing process, such as the carbon farming initiatives in Australia (Malone et al., 2017).

A regional soil carbon map with a resolution of 100 m (block support) was disaggregated into 10 m (point support) information for a farm in the Spring Ridge region in northern NSW. Application of the Bayesian ATPRK approach is presented and demonstrated within the study. In addition, a comparative analysis with the dissever algorithm, a regression-based downscaling method, is performed to assess the performance of Bayesian ATPRK in the context of the spatial downscaling of digital soil carbon mapping.

2 Theoretical background

Assuming we are observing a continuous Gaussian process, we denote the process as S(p) for locations p∈D, where D is the region of interest. For block support (B) data, where B⊂D, we presume that the observations are block averages as follows:

$$\begin{array}{}\text{(1)}& S\left(B\right)=\mathrm{|}B{\mathrm{|}}^{-\mathrm{1}}\underset{B}{\int }S\left(p\right)\mathrm{d}p,\end{array}$$

where |B| denotes the area of B.

The underlying Gaussian process S can be expressed as follows:

$$\begin{array}{}\text{(2)}& \mathit{S}=\sum _{k=\mathrm{0}}^l\mathbf{m}\mathit{\beta }+\mathit{\epsilon },\end{array}$$

where m is the design matrix of the trend function and β is the vector of coefficients; ε is the stochastic residual with a mean of zero and covariance function C(θ), where θ represent the parameters of the covariance function. Since the block-supported observations S_B and point support observations S_p are jointly Gaussian, the prediction problem can be introduced as follows:

$$\begin{array}{}\text{(3)}& \left[\begin{array}{c}{\mathit{S}}_{\mathrm{p}}\\ {\mathit{S}}_{\mathrm{B}}\end{array}\right]\sim N\left(\left[\begin{array}{c}\mathbf{m}\\ \mathbf{M}\end{array}\right]\mathit{\beta },\left[\begin{array}{ll}{\mathbf{C}}_{\mathrm{pp}}& {\mathbf{C}}_{\mathrm{pB}}\\ {\mathbf{C}}_{\text{Bp}}& {\mathbf{C}}_{\mathrm{BB}}\end{array}\right]\right),\end{array}$$

where C_pp is the variance–covariance matrix of S_p, C_BB is the variance–covariance matrix of S_B, and C_pB and C_Bp are the variance–covariance matrices between S_p and S_B and vice versa. m is the design matrix for fine-scale covariates of S_p.

Since the joint distribution is Gaussian, the optimal predictions at unsampled point locations can be given as follows:

$$\begin{array}{}\text{(4)}& \widehat{{\mathit{S}}_{\mathbf{p}}}={\mathbf{m}}^{\prime }\mathit{\beta }+{\mathbf{C}}_{\mathrm{pB}}^{\prime }{\mathbf{C}}_{\mathrm{BB}}\left({\mathit{S}}_{\mathrm{B}}-\mathbf{M}\mathit{\beta }\right),\end{array}$$

where $\widehat{{\mathit{S}}_{\mathbf{p}}}$ is the vector of predicted soil carbon at points.

And the variance–covariance matrix for the prediction error is the following:

$$\begin{array}{}\text{(5)}& \text{var}\left({\mathit{S}}_{\mathrm{p}}\mathrm{|}{\mathit{S}}_{\mathrm{B}}\right)={\mathbf{C}}_{\mathrm{PP}}-{\mathbf{C}}_{\mathrm{pB}}{\mathbf{C}}_{\mathrm{BB}}^{-\mathrm{1}}{{\mathbf{C}}_{\mathrm{pB}}}^{\prime }.\end{array}$$

This shows that the inference of S_p from S_B is straightforward and resembles a simple kriging formulation (Cressie and Wikle, 2011). However, the estimation of a point support variogram from the block support data is not a trivial task. Usually, the estimation is done via deregulation or deconvolution of the empirical variogram. Goovaerts (2008), Gotway and Young (2007), and Pardo-Igúzquiza and Atkinson (2007) proposed an iterative numerical procedure for the deconvolution of the block support variogram for irregular block support. Wang et al. (2015) suggested an integration approach, while Nagle et al. (2011) used a maximum likelihood approach. Gelfand et al. (2001) used the Markov chain Monte Carlo (MCMC) integration to estimate the point support variogram from block support data where the parameters of the point support variogram were given as non-informative priors. Truong et al. (2014) used a similar approach, but the parameters of the point support variogram were given as informative priors. In this paper, the variogram parameters were inferred directly from the data via Bayesian estimation with an objective of inferencing the uncertainty of the disaggregation process.

2.1 Bayesian inference

The likelihood function of the Gaussian process S can be given as $f\left(\mathit{S}\right)\mathrm{|},\mathit{\theta })=N(\mathbf{m}\mathit{\beta },\mathbf{C}\left(\mathit{\theta }\right))$. C is a function of the Euclidean distance h, with the vector of parameters θ. Then, the likelihood of predictive distribution at point locations can be given as follows:

$$\begin{array}{}\text{(6)}& f\left({\mathit{S}}_{\mathrm{p}}\mathrm{|}{\mathit{S}}_{\mathrm{B}}\right)=\int f({\mathit{S}}_{\mathrm{p}}\mathrm{|}{\mathit{S}}_{\mathrm{B}};\mathit{\beta },\mathit{\theta })f(\mathit{\beta },\mathit{\theta }\mathrm{|}{\mathit{S}}_{\mathrm{B}})\mathrm{d}\mathit{\beta }\mathrm{d}\mathit{\theta },\end{array}$$

where ${\mathit{S}}_{\mathrm{p}}\mathrm{|}{\mathit{S}}_{\mathrm{B}},\mathit{\beta },\mathit{\theta }$ is distributed as follows:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& N\left(\mathbf{M}\mathit{\beta }+{\mathbf{C}}_{\mathrm{Bp}}^{\prime }\left(\mathit{\theta }\right){\mathbf{C}}_{\mathrm{pp}}^{-\mathrm{1}}\left(\mathit{\theta }\right)({\mathit{S}}_{\mathrm{p}}-\mathbf{m}\mathit{\beta })\right),\\ & {\mathbf{C}}_{\mathrm{BB}}\left(\mathit{\theta }\right)-{\mathbf{C}}_{\mathrm{Bp}}^{\prime }\left(\mathit{\theta }\right){\mathbf{C}}_{\mathrm{pp}}^{-\mathrm{1}}\left(\mathit{\theta }\right){\mathbf{C}}_{\mathrm{Bp}}\left(\mathit{\theta }\right)).\end{array}\end{array}$$

Each entry of the distribution can be estimated through a Monte Carlo integration. Then, we can replace Eq. (3) with $\widehat{f}\left(({\mathit{S}}_{\mathrm{p}}{\mathit{S}}_{\mathrm{B}}{)}^{\prime }\mathrm{|}\mathit{\beta },\mathit{\theta }\right)$, where “hat” denotes a Monte Carlo integration. Then it is apparent, from the following, that:

$$\begin{array}{}\text{(8)}& \widehat{f}\left(({\mathit{S}}_{\mathrm{p}}{\mathit{S}}_{\mathrm{B}}{)}^{\prime }\mathrm{|}\mathit{\beta },\mathit{\theta }\right)=\widehat{f}\left(({\mathit{S}}_{\mathrm{p}}\widehat{{\mathit{S}}_{\mathrm{B}}}{)}^{\prime }\right)\mathrm{|}\mathit{\beta },\mathit{\theta }.\end{array}$$

To predict at ${\widehat{\mathit{S}}}_{\mathrm{p}}$ we require $f({\widehat{\mathit{S}}}_{\mathrm{p}},{\mathit{S}}_{\mathrm{B}})$, which is given by Eq. (3). Using $\widehat{f}\left(\right({\mathit{S}}_{\mathrm{p}},{\mathit{S}}_{\mathrm{B}}{)}^{\prime }\mathrm{|}\mathit{\beta },\mathit{\theta })$ and Eq. (8), we can obtain $f\left(\right({\mathit{S}}_{\mathrm{p}}\widehat{{\mathit{S}}_{\mathrm{B}}}{)}^{\prime })\mathrm{|}\mathit{\beta },\mathit{\theta }$ to sample ${\widehat{\mathit{S}}}_{\mathrm{p}}$.

Then, given the priors [β,θ], the Bayesian model is specified and a Monte Carlo fitting procedure can be carried out to maximise the following log likelihood function. A complete description of deducing the Bayesian area-to-point inference can be found in Gelfand et al. (2001).

$$\begin{array}{}\text{(9)}& \begin{array}{rl}L(\mathit{\beta },\mathit{\theta }\mathrm{|}{\mathit{S}}_{\mathrm{B}})=& {\displaystyle \frac{n}{\mathrm{2}}}\log \left(\mathrm{2}\mathit{\pi }\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\log \mathrm{|}{\mathbf{C}}_{\mathrm{BB}}\mathrm{|}\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{S}}_{\mathrm{B}}-\mathbf{M}\mathit{\beta }{)}^{\mathrm{T}}{\mathbf{C}}_{\mathrm{BB}}^{-\mathrm{1}}({\mathit{S}}_{\mathrm{B}}-\mathbf{M}\mathit{\beta }),\end{array}\end{array}$$

where C_BB is block-to-block covariance and can be approximately calculated as follows:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathbf{C}}_{\mathrm{BB}}(i,j)=\\ & \left\{\begin{array}{l}\frac{c_{\mathrm{1}}}{K}+\left(\frac{\mathrm{2}}{K(K-\mathrm{1})}\right){\sum }_{k=\mathrm{1}}^K{\sum }_{l>k}^K\left(c_{\mathrm{1}}+c_{\mathrm{0}}-{\mathit{\gamma }}_{\mathrm{p}}(s_k-s_l))\right),i\ne j\\ \frac{\mathrm{1}}{K^{\mathrm{2}}}\left(c_{\mathrm{1}}+c_{\mathrm{0}}-{\mathit{\gamma }}_{\mathrm{p}}(s_k-s_l))\right),i=j,\end{array}\right.\end{array}\end{array}$$

where i,j index are block-supported observations, k,l are index discretised points within a block, and K is the number of discretisation points per block (Truong et al., 2014). c₁, c₀ are the partial sill and the nugget component of the semivariogram γ_p.

Figure 1 The study area (left; sources: Esri, DigitalGlobe, GeoEye, i-cubed, USDA-FSA, USGS, AEX, Getmapping, AeroGRID, IGN, IGP, Swisstopo, and the GIS user community) and the coarse-scale topsoil soil organic carbon map (right).

3 Case study

3.1 Study area and coarse-scale soil carbon map

We illustrate the methodology by disaggregating a regional topsoil organic carbon map of the University of Sydney E.J. Holtsbaum Agricultural Research Station (Fig. 1). The landholding, also known as “Nowley”, is in the Spring Ridge district of the central/northwestern slopes of New South Wales (NSW), Australia. Its area is approximately 2300 ha and is managed as a mixed farming enterprise of cattle grazing throughout the year and wheat, barley, and canola in winter and sorghum and sunflower in the summer. Cropping is performed on soils derived from basalt parent materials. The average annual rainfall is 600 mm.

The available map is a NSW topsoil carbon map, which was generated using a local regression kriging technique (Somarathna et al., 2016). These maps are analogues of the nationally available soil carbon maps of Australia, which are available in separate layers corresponding to the depth intervals definition in the GlobalSoilMap product specifications (Arrouays et al., 2014). The standard soil depth layers are 0–5, 5–15, 15–30, 30–60, and 60–100 cm.

This NSW soil carbon map has a grid spacing of 100 m × 100 m, with block support of the same size. Hence, the map spatial resolution and the support are equal, which in turn means that the pixel values that make up the map represent a predicted mean for the area of each pixel. Commonly, the pixel values of a digital soil map are on point support, with the pixel value representing the value at a point that is usually at the centre of each pixel (Malone et al., 2013). Therefore, block kriging was used convert the maps to point-supported maps before the analysis.

The topsoil of the study area was represented by the depth interval of 0–7.5 cm, as this was the depth interval used in previous work conducted at this site for soil carbon auditing (de Gruijter et al., 2016). The respective coarse (NSW) map for the 0–7.5 cm depth interval was predicted by using mass-preserving splines (Malone et al., 2009) of the 0–5 and 5–15 cm standardised maps. This map will hereafter be referred to as the coarse-scale map. It is this map that subsequent descriptions of spatial downscaling will refer to. The experimental variogram of the block-supported coarse-scale map is given in Fig. 2.

Figure 2 Block support experimental variogram of the coarse-scale (100 m resolution) soil carbon map.

3.2 Environmental covariates

Elevation, topographic wetness index (TWI), gamma radiometric thorium, gamma radiometric potassium, Landsat band 4, and Landsat band 7 were used as environmental covariates in the spatial downscaling model. Landsat data have a resolution of 30 m. Since the coarse map was disaggregated into a 10 m resolution, the Landsat data were pan sharpened using the panchromatic band (15 m) and then resampled to 10 m.

The elevation data had been collected during a ground survey using an all-terrain vehicle with a 20 m swath. Then the data were processed using a local block kriging approach to produce elevation data at a 10 m resolution. Based on the developed digital elevation model (DEM), the TWI was derived. The gamma radiometric data were derived from aerial surveys (Minty et al., 2009), which were resampled to 10 m. A more detailed description of data collection and post-processing can be found in Malone et al. (2017).

3.3 Model parameter inferencing using the MCMC simulation

Suppose our observed soil carbon data (S) follows a second-order stationary Gaussian process which is characterised by the mean (trend) and the variogram where ${\mathit{S}}_{\mathrm{p}}\sim N({\mathit{S}}_{\mathrm{p}},{\mathit{\gamma }}_{\mathrm{p}})$ and where ${\mathit{\gamma }}_{\mathrm{p}}=(h,\mathit{\theta })$ is a function of Euclidian distance h and θ vector of parameters of Matérn variogram (Eq. 14). The block averages S_B have the same spatial mean as S_p but a different spatial structure (Gotway and Young, 2002). The block support covariance C_BB is also a function of h, with a parameter vector of θ, but has a different support.

$$\begin{array}{}\text{(11)}& C_{ij}=c_{\mathrm{0}}{\mathit{\delta }}_{ij}+c_{\mathrm{1}}\left[{\displaystyle \frac{\mathrm{1}}{{\mathrm{2}}^{\mathit{\upsilon }-\mathrm{1}}\mathrm{\Gamma }\left(\mathit{\upsilon }\right)}}{\left({\displaystyle \frac{h}{r}}\right)}^{\mathit{\upsilon }}{K}_{\mathit{\upsilon }}\left({\displaystyle \frac{h}{r}}\right)\right],\end{array}$$

where C_ij is the covariance between observation i and j, h represents the separation distance between i and j, δ_ij denotes the Kronecker delta (δ_ij=1 if i=j and δ_ij=0 when I≠j), c₀+c₁ signifies the sill variance, and K_v is the modified Bessel function of the second kind of order υ. Γ is the gamma function, r denotes the distance or range parameter and υ is the spatial smoothness. Calculating C_BB requires discretisation, and it was done using a regular grid pattern, with equal 10 m spacing between discretising points within a block. Then C_BB was calculated according to Eq. (11).

The priors were given as informative priors and considered as normally distributed. The coarse-scale soil carbon map was regressed against the selected covariates to estimate the approximate prior values for the coefficient of the linear trend (β). Gaussian priors were assigned for both the spatial trend (β) and covariance component (θ) of the spatial model. MCMC integration was carried out using the Metropolis–Hastings algorithm. Metropolis–Hastings is a computationally efficient algorithm that allows sampling from conditional distributions (Diggle et al., 1998). Calibration of the model was done using two parallel chains of MCMC with 100 000 iterations.

3.4 Spatial prediction

The mean 95 % upper and lower CI limits of the spatial model parameters were calculated using the values of the last 10 000 MCMC iterations. These calibrated parameters were used to predict the SOC content at the fine scale. The predictions were made onto the 10 m × 10 m points based on Eq. (4).

Figure 3 Schematic representation of dissever disaggregation technique.

3.5 Disaggregation using dissever

The proposed Bayesian model was also compared with the dissever regression-based disaggregation approach. Dissever is a regression-based downscaling approach proposed by Malone et al. (2012). Dissever operates in two steps, namely initialisation and iteration. At the initialisation step, the coarse map is resampled using the nearest neighbour resampling technique to match the resolution of available fine-scale covariates. Then the resampled map is regressed against the fine-scale covariates. Next, the regressed fine map is upscaled through averaging, and then the mass balance deviation (deviation factor – DF) from the coarse map is calculated. The regressed fine map is corrected using the DF, and then the algorithm progresses to the iteration step. In the iteration step, regression and upscaling continue until the mass balance deviation between the upscaled map and the original map is less than a set minimum (threshold). The following flow diagram (Fig. 3) is an illustrative description of the dissever algorithm. This differs from Bayesian ATPRK since it is not focused on addressing the change in support problem, and it is a fine-gridding exercise. A detailed and comprehensive account of dissever can be found in Malone et al. (2012; 2017).

Following the procedure mentioned above, the coarse-scale carbon map of the area was disaggregated into 10 m using dissever. The environmental covariates used in the Bayesian approach were similarly made available for dissever.

3.6 Independent validation

3.6.1 Validation data

The accuracy of the model was tested using an independent dataset, which consists of 127 soil samples from 0–7.5 cm soil cores (cross section area of 41 cm² and height of 7.5 cm). The samples are considered to have point support. These data were collected over two separate soil surveys during 2014 and 2015, using stratified random sampling. Stratification was based on the SOC prediction fields generated from digital soil mapping. A detailed description about the stratified random sampling process can be found in de Gruijter et al. (2015). The carbon content of the collected soil cores was measured using dry combustion by an Elementar Vario Max CNS macro elemental analyser (Elementar Analysensysteme GmbH, Langenselbold, Germany). The measured carbon contents are in carbon concentrations (g of C per 100 g of soil).

3.6.2 Validation statistics

The model accuracy was compared using mean squared error (MSE) and Lin's concordance correlation coefficient (CCC; Lin, 1989) between the disaggregated values and observed independent values at sampling points.

CCC is given by the following:

$$\begin{array}{}\text{(12)}& \text{CCC}={\displaystyle \frac{\mathrm{2}\mathit{\rho }{\mathit{\sigma }}_x{\mathit{\sigma }}_y}{{\mathit{\sigma }}_x^{\mathrm{2}}+{\mathit{\sigma }}_y^{\mathrm{2}}+{\left({\mathit{\mu }}_x-{\mathit{\mu }}_y\right)}^{\mathrm{2}}}},\end{array}$$

where ρ is the correlation coefficient, σ_x,σ_y are the variances of observed (x) and predicted (y) values, and μ_x,μ_y are the respective means. CCC is scaled between −1 and 1, with the latter implying perfect agreement and the former implying perfect reverse agreement.

MSE is given by the following:

$$\begin{array}{}\text{(13)}& \text{MSE}={\displaystyle \frac{\sum _{i=\mathrm{1}}^{n^{\prime }}{\left({\mathit{S}}_{\mathrm{p}}-{\widehat{\mathit{S}}}_{{\mathrm{p}}^{\prime }}\right)}^{\mathrm{2}}}{n^{\prime }}},\end{array}$$

where S_p is the observed values of SOC at independent sampling points, ${\widehat{\mathit{S}}}_{{\mathrm{p}}^{\prime }}$ is predicted SOC values at the fine scale, and n^′ is the number of independent SOC samples at a fine scale. MSE is an indication of the accuracy of predictions, while CCC indicates both accuracy and precision of the disaggregation exercise.

Figure 4 Correlations and posterior distributions of linear trend parameters of the spatial model.

Table 1 1. Posterior parameters of the linear trend model.

β₀, β₁, β₂, β₃, β₄, β₅, and β₆ are the coefficients of the linear trend for intercept, elevation, TWI, gamma radiometric thorium, gamma radiometric potassium, and Landsat band 4 and Landsat band 7 parameters, respectively.

4 Results and discussion

The use of the linear mixed model (LMM) as the spatial prediction model allows integration of fine-scale covariate information and spatial correlation in a single model. Unlike regression kriging in which regression and residual kriging are used in two phases of modelling, LMM is a one-step process. We used Bayesian inference techniques to estimate the parameters of LMM to deliver uncertainty estimates of the model parameters. This section provides a description of the posterior parameters of the spatial model followed by a comparison of spatially disaggregated maps using Bayesian ATPRK and dissever.

4.1 Bayesian parameter estimation of the spatial model

Figure 4 illustrates the posterior statistics of the linear trend parameters. The upper diagonal panel displays the correlation coefficients among the covariates, the lower diagonal depicts the scatter plots, while the diagonal shows the posterior distributions. The diagonal histograms are well defined, which indicates that all parameters of the linear trend are relatively well defined. The correlations among covariates are minimal, which indicates that there is no multicollinearity between the covariates.

Figure 5 displays the posterior distribution of Matérn variogram parameters and their correlations. Posterior values of the range parameter are similar to the prior values. However, the partial sill and nugget parameters are different from their priors. The posterior partial sill is higher than the prior sill, which agrees with the fact that the coarse-scale variability of soil carbon is lower compared to the finer scale (Fig. 5b).

Tables 1 and 2 provide the estimates of the parameters of the spatial model along with the 95 % confidence limits, which is an indication of the uncertainty of the parameters. The estimated posterior parameters have narrow CI limits, indicating the estimated mean provides a good approximation of the soil carbon content of the study area.

4.2 Spatial disaggregation using Bayesian ATPRK

Figure 6 displays the Bayesian ATPRK disaggregated soil carbon content of Nowley at 10 m, along with the 95 % confidence interval limits. Unlike other methods which only give a single estimate of mean soil carbon content, the Bayesian technique provides a range and/or interval of the most probable value with a 95 % confidence. These estimates are more reliable as they are generated using numerous repeated samplings. Accuracy and confidence in estimates are important parameters for planning, designing, and implementing soil management. Therefore, Bayesian ATPRK provides a comprehensive output of soil carbon estimates at the farm scale.

Figure 5 (a) Correlations and posterior distributions of variogram parameters (b) and the point and block support variograms.

Figure 6 Bayesian ATPRK disaggregated map with 95 % confidence limits.

4.3 Comparison of Bayesian ATPRK and dissever disaggregation methods

The disaggregated maps from both techniques were upscaled using a spatial overlay to assess the accuracy of disaggregation. CCC values were calculated for both techniques. Figure 7 shows the agreement between the cell values of the coarse map and the upscaled disaggregated maps.

Figure 7 Correlation plots between upscaled disaggregated values from Bayesian ATPRK, dissever, and original coarse-scale block values.

Concordance correlation between the original coarse values and upscaled product were 87 % and 96 % for Bayesian ATPRK and dissever, respectively. Dissever iterates until the agreement between upscaled and original coarse cell reaches a set minimum, while Bayesian ATPRK output is a result of one simulation considering both the correlation with the covariates and covariance structure of the residuals. This explains the stronger concordance associated with dissever.

Figure 8 illustrates the comparison between the dissever and Bayesian ATPRK disaggregated product with the coarse-scale source map. A closer examination reveals that Bayesian ATPRK product resembles the spatial pattern of SOC more closely than the dissever product, while dissever seemingly over-predicts the SOC content.

Table 2 Posterior parameters of the point support Matérn variogram model. CI refers to the confidence interval.

Figure 8 Comparing the disaggregated map with original coarse-scale SOC map.

Conversely, dissever is computationally more efficient compared to Bayesian ATPRK. Calculating the block-to-point and block-to-block distance matrixes is a computationally expensive procedure. The whole study area cannot be computed in a single run using a desktop computer. Therefore, the study area has to be split into tiles, and predictions were carried out for each of the tiles. This resulted in abrupt changes between two tiles in the Bayesian ATPRK disaggregated map.

4.4 Validation of disaggregated maps

Independent validation results between disaggregated maps using dissever and Bayesian ATPRK indicate that Bayesian ATPRK technique is slightly more accurate compared to dissever-produced output in general.

The CCC value of Bayesian ATPRK was 23 %, while that of dissever was 18 %. Greater CCC indicates more accurate and precise predictions. The mean squared error (MSE %) of the Bayesian ATPRK product was 0.45 while MSE for dissever was 0.49, which indicates that Bayesian ATPRK disaggregation is slightly more accurate than dissever disaggregation.

Although the validation CCC is low, it still is a sufficient approximation in terms of SOC as the SOC prediction accuracy generally remains low. The coarse map is a product of regression for the whole region, which also has a lower concordance correlation with the sampled data. Through the disaggregation process, the errors tend to propagate, widening the gap between the actual and predicted SOC content.

Overall, we have demonstrated using Bayesian ATPRK for downscaling coarse-scale rasters to more informative fine-scale rasters by using available fine-scale covariates. The study proves that Bayesian ATPRK provides more accurate results with uncertainty estimates of the disaggregation process. There are a few limitations to this study. Firstly, we did not deal with the uncertainty from all sources of measurement errors. The coarse-scale map always carries a certain degree of error and/or uncertainty. Also, there is an uncertainty associated with the splining of the maps. Both of these uncertainties can be treated as measurement errors and filtered by adding as a noise variance along the diagonal of the covariance matrix C_BB. Another way of correcting the bias is to use the measured data in the spatial model to reflect the actual ground conditions. Malone et al. (2017) demonstrated that the inclusion of the ground observations improves the accuracy of the disaggregated maps. A further limitation of our method is that the disaggregated map produced from Bayesian ATPRK has a block effect due to the final map being a collection of tiles produced separately due to the high computing demands of the model. This tiling effect can be eliminated by using high-performing and parallel-computing solutions to generate the fine-scale image in one go. Also, low-ranking representations (Wikle, 2010) of the covariance can also be used as a possible solution to the computational burden of the process.

Furthermore, although we have used the proposed Bayesian ATPRK in the DSM domain, this method can also be used to downscale coarse-scale satellite data. Even though there is a wide range of fine-resolution satellite images available in the present context, there is still a scarcity of fine-resolution data for certain parts of the world. Spatial downscaling is a solution to this data sparsity. There is a range of disaggregation and/or pan-sharpening techniques available in the literature. However, very few techniques can tackle the uncertainty of source images and the uncertainty of the disaggregation process. Unlike other methods, Bayesian ATPRK can explicitly help to tackle observation errors of the satellite images.

5 Conclusions and recommendations

The study demonstrates that Bayesian techniques can be effectively used to address the change in the support problems. Monte Carlo integration provides the estimates of point-supported covariance structures using block-supported data along with uncertainty estimates of the parameter inference. This study used informative priors of covariance parameters, as the soil carbon data were available for the study site. In the absence of prior information, priors can be given as non-informative priors.

Bayesian ATPRK provides an effective way of disaggregating coarse-scale SOC maps. The output of this study facilitates policy making and soil management activities by providing the maps with confidence limits, which lead to a wider understanding of the SOC content of the area. However, the method is computationally intensive, and the use of high-performing computing techniques is recommended for faster and smoother predictions.

Furthermore, the noise of coarse-scale observations can be filtered through incorporating the noise variance in the diagonal of C_BB. The use of point data (if available) in the spatial model can improve the accuracy of the disaggregation. Finally, the use of both techniques is recommended for lower uncertainty and improved accuracy of the disaggregation process.