Bayesian networks (BNs) have the ability to perform inference on uncertain variables given evidence on observed quantities, which makes them relevant mathematical tools for the updating of ground-motion fields based on strong-motion records or macroseismic observations. Therefore, the present article investigates the use of BN models of spatially correlated Gaussian random fields as an accurate and scalable method for the generation of ground-motion maps. The proposed BN model is based on continuous Gaussian variables, as opposed to discrete variables as in previous formulations, and it is built to account for cross-correlated ground-motion parameters as well as macroseismic observations. This approach is validated with respect to the analytical solution (i.e., conditional multivariate normal distributions), and it is also compared with the U.S. Geological Survey ShakeMap method, thus demonstrating a better ability to model jointly the interevent and intraevent error terms of ground-motion models. The scalability of the approach, that is, its capacity to be applied to large grids, is ensured by a grid subdivision strategy, which appears to be computationally efficient and accurate within an error rate of a fraction of percent. Finally, the BN implementation is demonstrated on a real-world example (the 2016 M-w 6.2 Kumamoto, Japan, foreshock), where vector-valued shake maps of cross-correlated intensity measures are generated, along with the integration of macroseismic observations.