The spatial interpolation of punctual radiometric data for the realization of natural radioactivity maps poses several challenges associated to the integration of information referred to different measurement techniques, measurement errors, detectors’ field of view and morphological features. Indeed, the elaboration of a unique cartographic product with an appropriate descriptive legend from laboratory, in-situ and airborne gamma-ray spectroscopy data cannot be pursued without critically dealing with some delicate issues. Gamma ray surveys allow for monitoring the spatial distribution of terrestrial radioelements (K, U and Th) and in turn provide valuable insights on geological mapping, structural geology and soil surveying. The high efficiency and relatively good spectral resolution make Sodium Iodide (NaI) detectors particularly suitable for real-time insitu and airborne measurements, respectively characterized by a spatial footprint of ∼m^2 and ∼10^5 m^2. On the other hand, Hyper Pure Germanium (HPGe) detectors, thanks to the optimal spectral resolution combined with low radiation background, are ideal for achieving very low uncertainties in gamma-ray measurements performed in laboratory on rock and soil samples. We dealt with these problematic aspects proposing operative solutions regarding the statistical treatment of analyzed datasets, the heterogeneous experimental uncertainties and the spatial resolution of the measurements. The results obtained from the rigorous study of statistical distributions and the spatial correlation were integrated based on appropriated geostatistical interpolators. Taking on the challenge to treat heterogeneous input uncertainties data, the degree of confidence associated with two different gamma-ray techniques is considered, giving value to the spatial data represented in the map. Multivariate spatial interpolation enhances the estimation of radioelements distribution taking advantage of the correlation existing between the under-sampled gamma-ray measurements and the continuous distributions of geological formations. The described methods were validated through several surveys that cover approximately 50000 km^2 of the Italian territory: specific case studies will be presented and discussed
