Goslar district in Lower-Saxony, Germany, was one of the most important mining and metal production areas in the world. The long time of mining, mineral processing and smelting activities in the area have enriched the forest soil with many heavy metals. Spatial distribution of soil heavy metals Pb, Zn, and As was investigated using geostatistical methods and geographic information system (GIS). Geostatistical Analyst extension of ArcGIS was used to explore the raw data, transform to normality and calculate the experimental variogram and its fitted model. Ordinary Kriging was applied to predict the spatial distribution of soil heavy metals. A GIS model was built to optimize the spatial distribution of heavy metals by selecting new soil samples in sites of maximum Standard Prediction Error. Comparison between predicted surfaces before and after taking the new samples was studied. Hotspot locations of polluted soil were investigated by creating specific geostatistical structures for each intensive contaminated site. These specific predicted surfaces of hotspot locations are merged into the general predicted surface of the same soil heavy metal.