State density distribution, in contrast to worst-case reachability, can be leveraged for safety-related problems to better quantify the likelihood of the risk for potentially hazardous situations. We developed a data-driven method to compute the density distribution of reachable states for nonlinear and even black-box systems. Our approach can estimate the set of all possible future states as well as their density. Moreover, we could perform online safety verification with probability ranges for unsafe behaviors to occur. We show that our approach can learn the density distribution of the reachable set more accurately with less data and quantify risks less conservatively and flexibly compared with worst-case analysis. We also study the use of such an approach in combination with model predictive control for verifiable safe path planning under uncertainties.