Scientific research often starts with data collection. However, many researchers pay insufficient attention to this first step in their research. The author, researcher at Wageningen University and Research, often had to conclude that the data collected by fellow researchers were suboptimal, or in some cases even unsuitable for their aim. One reason is that sampling is frequently overlooked in statistics courses. Another reason is the lack of practical textbooks on sampling. Numerous books have been published on the statistical analysis and modelling of data using R, but to date no book has been published in this series on how these data can best be collected. This book fills this gap. Spatial Sampling with R presents an overview of sampling designs for spatial sample survey and monitoring. It shows how to implement the sampling designs and how to estimate (sub)population- and space-time parameters in R. Key features Describes classical, basic sampling designs for spatial survey, as well as recently developed, advanced sampling designs and estimators Presents probability sampling designs for estimating parameters for a (sub)population, as well as non-probability sampling designs for mapping Gives comprehensive overview of model-assisted estimators Covers Bayesian approach to sampling design Illustrates sampling designs with surveys of soil organic carbon, above-ground biomass, air temperature, opium poppy Explains integration of wall-to-wall data sets (e.g. remote sensing images) and sample data Data and R code available on github Exercises added making the book suitable as a textbook for students The target group of this book are researchers and practitioners of sample surveys, as well as students in environmental, ecological, agricultural science or any other science in which knowledge about a population of interest is collected through spatial sampling. This book helps to implement proper sampling designs, tailored to their problems at hand, so that valuable data are collected that can be used to answer the research questions.

Spatial soil applications frequently involve binomial variables. If relevant environmental covariates are available, using a Bayesian generalized linear model (BGLM) might be a solution for mapping such discrete soil properties. The geostatistical extension, a Bayesian generalized linear geostatistical model (BGLGM), adds spatial dependence and is thus potentially better equipped. The objective of this work was to evaluate whether it pays off to extend from a BGLM to a BGLGM for mapping binary soil properties, evaluated in terms of prediction accuracy and modelling complexity. As motivating example, we mapped the presence/absence of the Pleistocene sand layer within 120 cm from the land surface in the Dutch province of Flevoland, using the BGLGM implementation in the R-package geoRglm. We found that BGLGM yields considerably better statistical validation metrics compared to a BGLM, especially with - as in our case - a large (n = 1,000) observation sample and few relevant covariates available. Also, the inferred posterior BGLGM parameters enable the quantification of spatial relationships. However, calibrating and applying a BGLGM is quite demanding with respect to the minimal required sample size, tuning the algorithm, and computational costs. We replaced manual tuning by an automated tuning algorithm (which eases implementing applications) and found a sample composition that delivers meaningful results within 50 h calculation time. With the gained insights and shared scripts spatial soil practitioners and researchers can - for their specific cases - evaluate if using BGLGM is feasible and if the extra gain is worth the extra effort. Highlights Does adding spatial correlation to a Bayesian GLM for mapping a binary soil variable pay off? We aim to make spatial Bayesian hierarchical modelling accessible for pedometricians. Most hierarchical models work well when enough observations are provided, even without covariates. Including spatial correlation might sometimes be worth the extra effort and computational costs.

For many decades, soil scientists have produced spatial estimates of soil properties using statistical and non-statistical mapping models. Commonly in soil mapping studies the map quality is assessed through pairwise comparison of observed and predicted values of a soil property, from which statistical indices summarizing the quality of the entire map are computed. Often these indices are based on average error and correlation statistics. In this study, we recommend a more appropriate and effective method of map evaluation by means of Taylor and solar diagrams. Taylor and solar diagrams are summary diagrams exploiting the relationship between statistical indices to visualize differentiable aspects of map quality into a single plot. An important advantage over current map quality evaluation is that map quality can be assessed from the combined effect of a few statistical quantities, not just on the basis of a single index or list of indices. We illustrate the use of common statistical indices and their combination into summary diagrams with a simulation study and two applications on soil data. In the simulation study nine maps with known statistical properties are produced and evaluated with tables and summary diagrams. In the first case study with soil data, change in the quality of a large-scale topsoil organic carbon map is tracked for a number of permutations in the mapping model parameters, whereas in the second case study several maps of topsoil organic carbon content for the same area, made by various statistical and non-statistical models, are compared and evaluated. We consider that in all cases better insights in map quality are obtained with summary diagrams, instead of using a single index or an extensive list of indices. This underpins the importance of using integrated summary graphics to communicate on quantitative map quality so as to avoid excessive trust that a single map quality index may suggest.

Mapping of environmental variables often relies on map accuracy assessment through cross-validation with the data used for calibrating the underlying mapping model. When the data points are spatially clustered, conventional cross-validation leads to optimistically biased estimates of map accuracy. Several papers have promoted spatial cross-validation as a means to tackle this over-optimism. Many of these papers blame spatial autocorrelation as the cause of the bias and propagate the widespread misconception that spatial proximity of calibration points to validation points invalidates classical statistical validation of maps. We present and evaluate alternative cross-validation approaches for assessing map accuracy from clustered sample data. The first method uses inverse sampling-intensity weighting to correct for selection bias. Sampling-intensity is estimated by a two-dimensional kernel approach. The two other approaches are model-based methods rooted in geostatistics, where the first assumes homogeneity of residual variance over the study area whilst the second accounts for heteroscedasticity as a function of the sampling intensity. The methods were tested and compared against conventional k-fold cross-validation and blocked spatial cross-validation to estimate map accuracy metrics of above-ground biomass and soil organic carbon stock maps covering western Europe. Results acquired over 100 realizations of five sampling designs ranging from non-clustered to strongly clustered confirmed that inverse sampling-intensity weighting and the heteroscedastic model-based method had smaller bias than conventional and spatial cross-validation for all but the most strongly clustered design. For the strongly clustered design where large portions of the maps were predicted by extrapolation, blocked spatial cross-validation was closest to the reference map accuracy metrics, but still biased. For such cases, extrapolation is best avoided by additional sampling or limitation of the prediction area. Weighted cross-validation is recommended for moderately clustered samples, while conventional random cross-validation suits fairly regularly spread samples.

For decades scientists have produced maps of biological, ecological and environmental variables. These studies commonly evaluate the map accuracy through cross-validation with the data used for calibrating the underlying mapping model. Recent studies, however, have argued that cross-validation statistics of most mapping studies are optimistically biased. They attribute these overoptimistic results to a supposed serious methodological flaw in standard cross-validation methods, namely that these methods ignore spatial autocorrelation in the data. They argue that spatial cross-validation should be used instead, and contend that standard cross-validation methods are inherently invalid in a geospatial context because of the autocorrelation present in most spatial data. Here we argue that these studies propagate a widespread misconception of statistical validation of maps. We explain that unbiased estimates of map accuracy indices can be obtained by probability sampling and design-based inference and illustrate this with a numerical experiment on large-scale above-ground biomass mapping. In our experiment, standard cross-validation (i.e., ignoring autocorrelation) led to smaller bias than spatial cross-validation. Standard cross-validation was deficient in case of a strongly clustered dataset that had large differences in sampling density, but less so than spatial cross-validation. We conclude that spatial cross-validation methods have no theoretical underpinning and should not be used for assessing map accuracy, while standard cross-validation is deficient in case of clustered data. Model-free, design-unbiased and valid accuracy assessment is achieved with probability sampling and design-based inference. It is valid without the need to explicitly incorporate or adjust for spatial autocorrelation and perfectly suited for the validation of large scale biological, ecological and environmental maps.

It is commonly accepted that an estimated soil variogram can be transferred to another similar area for deriving the tolerable spacing of a sampling grid or, more generally, the sample size, given a requirement on the quality of the soil property map of the recipient area. The quality of the derived tolerable grid spacing depends on how similar the population variograms of the donor area and recipient area are. In practice we are uncertain about the variograms of both areas due to sampling errors. Ideally, the uncertainty about the variogram of the donor area is accounted for in deriving the tolerable grid spacing. To assess the transferability, we should also account for uncertainty in the estimated variogram of the recipient area. In this study the transferability of variograms of soil pH, P, Mg and K is analysed for three grassland fields in Ireland, which are similar in soil‐forming factors. One field served as donor area, the other two as recipient area. For all three fields and for each soil property, 500 variograms were sampled from the posterior distribution of the variogram parameters. Results showed that the estimated variogram parameters of the recipient fields differed largely from those of the transferred variograms. The ranges of estimated mean kriging variance values for the various grid spacings, as obtained with the two sets of variograms (one set of the donor field, one set of the recipient field), did not overlap. Even after scaling the transferred variogram with an estimate of the variance of the recipient field, the transferred variogram was of no use for determining the tolerable grid spacing. The difference in the variograms can possibly be explained by the difference in historical land use.

Several misconceptions about the design-based approach for sampling and statistical inference, based on classical sampling theory, seem to be quite persistent. These misconceptions are the result of confusion about basic statistical concepts such as independence, expectation, and bias and variance of estimators or predictors. These concepts have a different meaning in the design-based and model-based approach, because they consider different sources of randomness. Also, a population mean is still often confused with a model mean, and a population variance with a model-variance, leading to invalid formulas for the variance of an estimator of the population mean. In this paper the fundamental differences between these two approaches are illustrated with simulations, so that hopefully more pedometricians get a better understanding of this subject. An overview is presented of how in the design-based approach we can make use of knowledge of the spatial structure of the study variable. In the second part, new developments in both the design-based and model-based approach are described that try to combine the strengths of the two approaches.

If a map is constructed through prediction with a statistical or non‐statistical model, the sampling design used for selecting the sample on which the model is fitted plays a key role in the final map accuracy. Several sampling designs are available for selecting these calibration samples. Commonly, sampling designs for mapping are compared in real‐world case studies by selecting just one sample for each of the sampling designs under study. In this study, we show that sampling designs for mapping are better compared on the basis of the distribution of the map quality indices over repeated selection of the calibration sample. In practice this is only feasible by subsampling a large dataset representing the population of interest, or by selecting calibration samples from a map depicting the study variable. This is illustrated with two real‐world case studies. In the first case study a quantitative variable, soil organic carbon, is mapped by kriging with an external drift in France, whereas in the second case a categorical variable, land cover, is mapped by random forest in a region in France. The performance of two sampling designs for mapping are compared: simple random sampling and conditioned Latin hypercube sampling, at various sample sizes. We show that in both case studies the sampling distributions of map quality indices obtained with the two sampling design types, for a given sample size, show large variation and largely overlap. This shows that when comparing sampling designs for mapping on the basis of a single sample selected per design, there is a serious risk of an incidental result.

Area-to-point kriging (ATPK) is a geostatistical method for creating high-resolution raster maps using data of the variable of interest with a much lower resolution. The data set of areal means is often considerably smaller ($$<\,50 $$ observations) than data sets conventionally dealt with in geostatistical analyses. In contemporary ATPK methods, uncertainty in the variogram parameters is not accounted for in the prediction; this issue can be overcome by applying ATPK in a Bayesian framework. Commonly in Bayesian statistics, posterior distributions of model parameters and posterior predictive distributions are approximated by Markov chain Monte Carlo sampling from the posterior, which can be computationally expensive. Therefore, a partly analytical solution is implemented in this paper, in order to (i) explore the impact of the prior distribution on predictions and prediction variances, (ii) investigate whether certain aspects of uncertainty can be disregarded, simplifying the necessary computations, and (iii) test the impact of various model misspecifications. Several approaches using simulated data, aggregated real-world point data, and a case study on aggregated crop yields in Burkina Faso are compared. The prior distribution is found to have minimal impact on the disaggregated predictions. In most cases with known short-range behaviour, an approach that disregards uncertainty in the variogram distance parameter gives a reasonable assessment of prediction uncertainty. However, some severe effects of model misspecification in terms of overly conservative or optimistic prediction uncertainties are found, highlighting the importance of model choice or integration into ATPK.

This study investigates sampling design for mapping soil classes based on multiple environmental features associated with the soil classes. Two types of sampling design for calibrating the prediction models are compared: conditioned Latin hypercube sampling (CLHS) and feature space coverage sampling (FSCS). Simple random sampling (SRS), which does not utilize the environmental features, is added as a reference design. The sample sizes used are 20, 30, 40, 50, 75, and 100 points, and at each sample size 100 sample sets were drawn using each of the three types of design. Each of these sample sets was then used to calibrate three prediction models: random forest (RF), individual predictive soil mapping (iPSM), and multinomial logistic regression (MLR). These sampling designs were compared based on the overall accuracy of predicted soil class maps obtained by these three prediction methods. The comparison was conducted in two study areas: Ammertal (Germany) and Raffelson (USA). For each of these two areas a detailed legacy soil class map is available. These soil class maps were used as references in a simulation study for the comparison. Results of both study areas show that on average FSCS outperforms CLHS and SRS for all three prediction methods. The difference in estimated medians of overall accuracy with CLHS and SRS was marginal. Moreover, the variation in overall accuracy among sample sets of the same size was considerably smaller for FSCS than that for CLHS. These results in the two study areas suggest that FSCS is a more effective sampling design.

Machine learning techniques are widely employed to generate digital soil maps. The map accuracy is partly determined by the number and spatial locations of the measurements used to calibrate the machine learning model. However, determining the optimal sampling design for mapping with machine learning techniques has not yet been considered in detail in digital soil mapping studies. In this paper, we investigate sampling design optimization for soil mapping with random forest. A design is optimized using spatial simulated annealing by minimizing the mean squared prediction error (MSE). We applied this approach to mapping soil organic carbon for a part of Europe using subsamples of the LUCAS dataset. The optimized subsamples are used as input for the random forest machine learning model, using a large set of readily available environmental data as covariates. We also predicted the same soil property using subsamples selected by simple random sampling, conditioned Latin Hypercube sampling (cLHS), spatial coverage sampling and feature space coverage sampling. Distributions of the estimated population MSEs are obtained through repeated random splitting of the LUCAS dataset, serving as the population of interest, into subsets used for validation, testing and selection of calibration samples, and repeated selection of calibration samples with the various sampling designs. The differences between the medians of the MSE distributions were tested for significance using the non-parametric Mann-Whitney test. The process was repeated for different sample sizes. We also analyzed the spread of the optimized designs in both geographic and feature space to reveal their characteristics. Results show that optimization of the sampling design by minimizing the MSE is worthwhile for small sample sizes. However, an important disadvantage of sampling design optimization using MSE is that it requires known values of the soil property at all locations and as a consequence is only feasible for subsampling an existing dataset. For larger sample sizes, the effect of using an MSE optimized design diminishes. In this case, we recommend to use a sample spread uniformly in the feature (i.e. covariate) space of the most important random forest covariates. The results also show that for our case study, cLHS sampling performs worse than the other sampling designs for mapping with random forest. We stress that comparison of sampling designs for calibration by splitting the data just once is very sensitive to the data split that one happens to use if the validation set is small.