Data collected by space probes and Earth-based observations allow to investigate multiple aspects of planetary bodies in the Solar System. Typical spatial data utilized for geological and geophysical global scale investigations include imagery basemaps, digital elevation models, feature density maps, geochemical, magnetic and gravity anomaly maps. Different datasets often present local correlations, and their simultaneous use can help unveil the wide diversity of geological and geophysical processes that contribute to the planetary evolution. A geostatistical approach can be applied to quantify the level of correlation and anticorrelation between the different datasets. We developed a method based on a linear regression analysis of multiple planetary data, producing georeferenced maps that can be imported into a Geographic Information System (GIS) environment. We show how this geostatistical method can be used to quantify the level of correlation between different spatial data, and to investigate the geological and geophysical processes of planetary bodies. We present the results of the geostatistical method for the Mercury science case, using datasets returned by the MESSENGER mission.