Soil contamination by heavy metals is of particular concern, due to the direct negative impact on crop yield, food quality and human health. Although the conventional approach to monitor heavy metals relies on field sampling and lab analysis, the proliferation in the use of portable spectrometers has reduced the cost and time of investigation. However, discrepancies in spectral data from different spectrometers increase the modeling time and undermine the model accuracy for spatial mapping. This study, therefore, took advantage of the readily accessible Landsat 7 data to predict and map the spatiotemporal distribution of ten heavy metals (i.e., Sb, Pb, Ni, Mn, Hg, Cu, Cr, Co, Cd and As) over a 640 km2 area in Belgium. The Land Use/Cover Area Frame Survey (LUCAS) database of a region in north-eastern Belgium was used to retrieve variation in heavy metals concentrations over time and space, using the Landsat 7 imagery for four single dates in 2009, 2013, 2016 and 2020. Three regression methods, namely, partial least squares regression (PLSR), random forest (RF) and support vector machine (SVM) were used to model and predict the heavy metal concentrations for 2009. By comparing these models unbiasedly, the best model was selected for predicting and mapping the heavy metal distributions for 2013, 2016 and 2020. RF turned out to be the optimal model for 2009 with a coefficient of determination of prediction (R2P) and residual prediction deviation of prediction (RPDP) ranging from 0.62 to 0.92, and 1.23 to 2.79, respectively. The measured heavy metal distributions along the river floodplains, at the highlands and in the lowlands, were generally high, compared to their RF spatiotemporal predictions, which decreased over time. Increasing moisture contents in the floodplains adjacent to the river channels and the lowlands were the primary contributors to the reduction in the satellite reflectance spectra. However, topsoil erosion from rainfall, snowmelt as well as wind into the lowlands could have influenced the reduction in heavy metal spatiotemporal predicted values over time in the highlands. The spatiotemporal prediction maps produced for the heavy metals for the four different years revealed a good spatial similarity and consistency with the measured maps for 2009, which indicates their stability over the years.