This paper provides a new methodological framework to generate empirical ground shaking scenarios, designed for engineering applications and civil protection planning. The methodology is useful both to reconstruct the ground motion pattern of past events and to generate future shaking scenarios, in regions where strong-motion datasets from multiple events and multiple stations are available. The proposed methodology combines (1) an ad-hoc nonergodic ground motion model (GMM) with (2) a spatial correlation model for the source region-, site-, and path-systematic residual terms, and (3) a model of the remaining aleatory error to take into account for directivity effects. The associated variability is a function of the type of scenario generated (bedrock or site, past or future event) and it is minimal for source areas where several events have occurred and for sites where recordings are available. In order to develop the region-specific fully nonergodic GMM and to compute robust estimation of the residual terms, the approach is calibrated on a highly dense dataset compiled for the area of central Italy. Example tests demonstrate the validity of the approach, which allows to simulate acceleration response spectra at unsampled sites, as well as to capture peculiar physical features of ground motion patterns in the region. The proposed approach could be usefully adopted for data-driven simulations of ground shaking maps, as alternative or complementary tool to physic-based and stochasticbased approaches.