Properly modeling and predicting the scattering response of a metasurface is a particularly challenging task. This has been shown to be especially difficult if the metasurface supports both local and nonlocal interactions, in the form of lattice coupling effects, multipolar contributions, or bianisotropic responses. So far, existing methods for addressing this problem often lack generality, being limited to normal incidence or applicable only to dipolar responses under oblique incidence. More general models, such as Green's function, T-matrix, or generalized susceptibility approaches, do exist; however, their application becomes challenging, particularly when analyzing the angular scattering response of complex unit cells or inhomogeneous background media. We overcome these limitations by providing a rigorous and comprehensive formalism that accommodates both oblique incidence and the presence of different superstrate and substrate. This is achieved by extending our existing metasurface modeling framework to account for nonlocal and multipolar contributions up to the any desired order and properly accounting for the scattering effects due to an inhomogeneous background medium. Additionally, our method is based on exact spherical multipole decomposition, which intrinsically accounts for toroidal contributions. We demonstrate the effectiveness of our approach by modeling the response of several dielectric and plasmonic metasurfaces that exhibit sharp spectral features including bound states in the continuum. Overall, our formalism yields excellent agreement with full-wave simulations.