National-scale soil datasets exhibit variation over widely disparate spatial scales. Geological variation is an important source at the coarsest scale, and one in which exhaustive information is commonly available, in geological maps. Superimposed on this is continuous spatial variation caused by factors such as relief, vegetation and diffuse ‘background’ pollution. Further variation is caused by locally distinct factors such as point pollution from industrial sites or occasional geological anomalies. In this paper, we propose a single statistical model to encompass all of these effects which we describe as ‘geological variation’, ‘continuous spatial variation’ and ‘local anomalies'. In our model, the geological and continuous spatial variation are described, respectively, by the fixed and random effects of a linear mixed model (LMM) and the local anomalies lead to observations which are spatial outliers with respect to the LMM. We fit the model to a survey of 1887 observations of cadmium concentration in soil (Cd) collected on an incomplete regular grid across the French metropolitan territory (550 000 km2) and use it to predict Cd across France. We find that (i) it is not possible to fit a valid model—in terms of cross-validation statistics—of Cd variation unless the effects of local anomalies are identified and separated from the larger-scale processes; (ii) the LMM is not valid if the outliers are merely discarded but a valid model does result if the outliers are winsorized. On the basis of these findings we suggest a practical robust algorithm for national-scale spatial analysis.