For probability density estimation on Riemannian manifolds, many methods have been proposed, parametric or non-parametric. However, whether they can be implemented for general manifolds is in question. Here we propose a kernel-based method for density estimation and sampling of singular measures supported on Riemannian submanifolds. Based on the theory of Pelletier (2005), the kernels are functions of the Riemannian geodesic distance on the submanifold, which we estimate from a neighborhood graph of data. To improve the estimate of geodesic distances, we densify the graph by numerical retraction. With an estimated kernel bandwidth, we can sample from the model using the same retraction. Our method applies to both constraint manifolds and data manifolds, and is more efficient compared with diffusion-based methods and eigenfunction projection.