Past seismic events worldwide demonstrated that damage and death toll depend on both the strong ground motion (i.e., source effects) and the local site effects. The variability of earthquake ground motion distribution is caused by the local stratigraphic and/or topographic setting and buried morphologies (e.g., irregular sub-interface between soft and stiff soils) that can give rise to amplification and resonances with respect to the ground motion expected at the reference site. Therefore, local site conditions can affect an area with damage related to the full collapse or loss in functionality of facilities, roads, pipelines, and other lifelines. To this concern, the near-real-time prediction of ground motion variation over large areas is a crucial issue to support the rescue and operational interventions. A machine learning approach was adopted to produce ground motion prediction maps considering both stratigraphic and morphological conditions. A set of about 16 000 accelerometric data points and about 46 000 geological and geophysical data points was retrieved from Italian and European databases. The intensity measures of interest were estimated based on nine input proxies. The adopted machine learning regression model (i.e., Gaussian process regression) allows for improving both the precision and the accuracy in the estimation of the intensity measures with respect to the available near-real-time prediction methods (i.e., ground motion prediction equation and ShakeMaps). In addition, maps with a 50 m × 50 m resolution were generated, providing a ground motion variability in agreement with the results of advanced numerical simulations based on detailed subsoil models.