|Abstract:|| This work offers a novel methodological framework to address the problem of generating data-driven earthquake shaking fields at different vibration periods, which are key to support decision making and civil protection planning. We propose to analyse the entire profiles of spectral accelerations and project their information content to unsampled locations in the system, based on the theory of Object Oriented Spatial Statistics (O2S2). The proposed methodology combines a non-ergodic ground motion model (GMM) with a fully functional model for the residual term, the latter consisting of (i) the spatially-varying systematic effects due to source, site and path, and (ii) the remaining aleatory error. The proposed methodology allows to generate multiple shaking scenarios conditioned on the data, jointly and consistently for all the vibration periods, overcoming the intrinsic limitations of existing multivariate approaches to the problem. The approach is tested on a vast dataset of ground motion records collected in the study-area of the Po Plain (Northern Italy), for which a region-specific fully non-ergodic GMM was previously calibrated. Our validation tests demonstrate the potentiality of the approach, which is capable to effectively simulate spectral acceleration profiles, while keeping the ability to capture the main physical features of ground motion patterns in the region.