Seismic ground motion in urban environment is generally modified by the presence of buildings. This phenomenon, referred to as Site-City Interaction, is mainly due to the radiated waves emitted through the soil by vibrating structures. Therefore, the urban seismic wavefield is altered from the free field motion with zones of higher and lower energy, i.e. largest and smaller ground motion amplitude, respectively. Because of the complexity of the phenomenon, modelling and simulations of large dense cities might be performed with large computational costs. In this paper, a random urban wavefield is obtained through a hierarchical multi-scale approach that defines representative areas obtained by homogenization of cluster of buildings at local-scale. An equivalent single oscillator is hence derived for each representative areas by considering uncertain parameters such as the fundamental period. Therefore, a proposed model combining discrete model and wave propagation model is used to evaluate the interactions among the several representative areas at city-scale to derive the seismic ground motion map for urban areas in terms of 50% fractile peak acceleration. Monte Carlo Simulation study is performed for validation purposes. Finally, the proposed framework is applied to a real urban area using data from LIDAR maps.