Pedometricians have spent a lot of effort on mapping soil types and basic soil properties. However, end-users typically need a more elaborate soil quality index for land management. Soil quality indices are typically derived from multiple individual soil properties by evaluating whether specific criteria are met. If this is based on individually mapped soil properties, then an important consequence is that cross-correlations between soil properties are ignored. This makes it impossible to quantify the uncertainties associated with the mapped indices. The objective of this study was to map a soil potential multifunctionality index for agriculture (Agri-SPMI) over a 12,125 km(2) study region located along the French Mediterranean coast to help urban planners preserve soils of the highest quality. The index considered the ability of soils to fulfil four functions under five land use scenarios. A binary map represented each soil function fulfilment for a given scenario. The final soil quality index map was the sum of the 20 binary maps. A regression cokriging model was developed to map the basic soil properties first individually from legacy soil data and spatial soil covariates using a random forest algorithm, and next, interpolate the residuals using cokriging and the linear model of coregionalisation. The mapping uncertainties of soil properties were propagated by calculating the soil quality index over 300 stochastic simulations of soil properties derived from the linear models of coregionalisation. Results showed a poor prediction accuracy of the quality index, mainly because some soil properties were poorly predicted (notably available water capacity and coarse fragments) and used in combination with extreme thresholds that defined land suitability. Overall, the uncertainty was correctly quantified because the stochastic simulations reproduced the width of the observed distribution well, but the shapes of the distributions differed considerably from those of the observations. We envisage some ways for improvement, such as creating probability maps instead of the mean from simulations, and changing the prediction support from point to area.