This report describes the potential and functionality of software for spatial analysis, prediction and stochastic simulation of continuous soil properties using data from the Dutch Soil Information System (BIS). A geostatistical framework and R codes were developed. The geostatistical model of a soil property has a deterministic component representing the mean value within a soil category, and a stochastic component of standardized residuals. The standardized residuals are interpolated or simulated based on the simple kriging system. The software was tested in four case studies: exchangeable soil pH, clay content, organic matter content and Mean Spring Water table depth (MSW). It is concluded that the geostatistical framework and R codes developed in this study enable to predict values of continuous soil properties spatially, and to quantify the inaccuracy of these predictions. The inaccuracy of a spatial prediction at a certain location is quantified by the kriging variance, which can be interpreted as an indication of the uncertainty about the true value.