In this paper we present maps of K, eU, and eTh abundances of Elba Island (Italy) obtained with a multivariate spatial interpolation of airborne {gamma}-ray data using the constraints of the geologic map. The radiometric measurements were performed by a module of four NaI(Tl) crystals of 16 L mounted on an autogyro. We applied the collocated cokriging (CCoK) as a multivariate estimation method for interpolating the primary under-sampled airborne {gamma}-ray data considering the well-sampled geological information as ancillary variables. A random number has been assigned to each of 73 geological formations identified in the geological map at scale 1:10,000. The non-dependency of the estimated results from the random numbering process has been tested for three distinct models. The experimental cross-semivariograms constructed for radioelement-geology couples show well-defined co-variability structures for both direct and crossed variograms. The high statistical correlations among K, eU, and eTh measurements are confirmed also by the same maximum distance of spatial autocorrelation. Combining the smoothing effects of probabilistic interpolator and the abrupt discontinuities of the geological map, the results show a distinct correlation between the geological formation and radioactivity content. The contour of Mt. Capanne pluton can be distinguished by high K, eU and eTh abundances, while different degrees of radioactivity content identify the tectonic units. A clear anomaly of high K content in the Mt. Calamita promontory confirms the presence of felsic dykes and hydrothermal veins not reported in our geological map. Although we assign a unique number to each geological formation, the method shows that the internal variability of the radiometric data is not biased by the multivariate interpolation.