We propose two algorithms that use both models and datasets to estimate angular power spectra from channel covariance matrices in massive MIMO systems. The first algorithm is an iterative fixed-point method that solves a hierarchical problem. It uses model knowledge to narrow down candidate angular power spectra to a set that is consistent with a measured covariance matrix. Then, from this set, the algorithm selects the angular power spectrum with minimum distance to its expected value with respect to a Hilbertian metric learned from data. The second algorithm solves an alternative optimization problem with a single application of a solver for nonnegative least squares programs. By fusing information obtained from datasets and models, both algorithms can outperform existing approaches based on models, and they are also robust against environmental changes and small datasets.