This paper introduces the Wavelet Phase Harmonics (WPH) statistics. They are interpretable low-dimensional statistics which describe 2D non-Gaussian density fields. These statistics are built from WPH moments, which have been recently introduced in data science and machine learning community. In this paper, we applied WPH statistics to projected matter density fields from the Quijote N-body simulations. We find by computing the Fisher information matrix, that the WPH statistics can place more stringent constraints on 5 cosmological parameters when compared to the combination of power-spectrum and bi-spectrum. We also use the WPH statistics to successfully generate from a maximum entropy model new 2D density fields that reproduce the PDF, mean, power-spectrum, bispectrum and the Minkowski functionals of the input density fields. While separate methods have been proven very efficient for parameter estimations and statistical syntheses for for large scale structure, WPH statistics are the first statistics that can both achieve a more stringent cosmological parameter constraint, and produce a sufficiently accurate simulation of the Universe while being interpretable.