This contribution deals with the design and analysis of multiscale numerical methods for modeling floods in complex urban environments, where the challenge arises from the stark contrast between the scale of the flow domain (at the district or city level) and the smaller-scale features, such as buildings or walls, that must be accurately represented. Additionally, the irregular geometry of urban domains leads to solutions with numerous corner singularities. Starting with a linear diffusion problem in a domain containing a large number of polygonal perforations, representing structural features of urban areas, we propose a low-dimensional approximation space based on a coarse polygonal partitioning of the domain. Similar to other multiscale numerical methods, such as the Multiscale Finite Element Method or the Multiscale Hybrid-Mixed Method, this coarse space, is spanned by locally discrete harmonic basis functions. Along the subdomain boundaries, the basis functions are piecewise polynomial. Due to its energy minimization property, the coarse space naturally incorporates singular basis functions. This allows us to derive an original H^1 error estimate that remains robust with respect to both the geometry of the perforations and the singularities induced by the geometry. To tackle urban flood simulations, we develop a mass-lumped Galerkin discretization of the nonlinear Diffusive Wave model using the coarse space described earlier. Numerical experiments demonstrate that the method achieves reasonable accuracy on coarse grids while offering a substantial speedup, up to a factor of 100, compared to fine-scale simulations.