Soil surveys with line-scanning platforms appear to have great advantages over the traditional methods used to collect soil information for the development of field-scale soil mapping and applications. These carry VNIR (visible and near infrared) spectrometers and have been used in recent years extensively for the assessment of soil fertility at the field scale, and the delineation of site-specific management zones (MZ). A challenging feature of VNIR applications in precision agriculture (PA) is the massiveness of the derived datasets that contain point predictions of soil properties, and the interpolation techniques involved in incorporating these data into site-specific management plans. In this study, fixed-rank kriging (FRK) geostatistical interpolation, which is a flexible, non-stationary spatial interpolation method especially suited to handling huge datasets, was applied to massive VNIR soil scanner data for the production of useful, smooth interpolated maps, appropriate for the delineation of site-specific MZ maps. Moreover, auxiliary Sentinel-2 data-based biophysical parameters NDVI (normalized difference vegetation index) and fAPAR (fraction of photosynthetically active radiation absorbed by the canopy) were included as covariates to improve the filtering performance of the interpolator and the ability to generate uniform patterns of spatial variation from which it is easier to receive a meaningful interpretation in PA applications. Results from the VNIR prediction dataset obtained from a pivot-irrigated field in Albacete, southeastern Spain, during 2019, have shown that FRK variants outperform ordinary kriging in terms of filtering capacity, by doubling the noise removal metrics while keeping the computation cost reasonably low. Such features, along with the capacity to handle a large volume of spatial information, nominate the method as ideal for PA applications with massive proximal and remote sensing datasets.