Spatial heterogeneity in soil properties has an impact on crop response. There is a growing demand for rapid and non-invasive acquisition of fine-scale information on soil and plant variation for site-specific management. Proximal (Electromagnetic Induction sensing (EMI), Ground Penetrating Radar imaging (GPR), hyperspectral spectroscopy (HS) and remote sensing (RS) can complement direct sampling. However, sensor data fusion techniques jointly analysing data from different sources, are still being developed. In this work, homogenous zones within an 1.5-ha arable field were delineated through an approach combining EMI, GPR, HS and RS. The multi-sensor data were split into four groups: 1) bulk electical conductivity (EC) from EMI data, 2) amplitude of GPR signal, 3) the first principal components related to five bands (green, yellow, red, rededge, near-infrared) of hyperspectral reflectance and 4) the vegetation indices (NDVI, NDRE and Near/Green) calculated from the remote sensing images. The data of each group were separately analysed and interpolated at the nodes of a 0.5 x 0.5 m-grid by using cokrigingor kriging. To obtain spatially contiguous clusters, a combined approach, based on multivariate geostatistics and a non parametric density function algorithm of clustering, was applied to the overall multivariate data set of the estimates. This approach proved to efficiently analyze data from different sources and to be suitable for delineation of homogeneous sub-field zones in site-specific management.