Several algorithms have been proposed to calculate the spatial entanglement spectrum from high order Renyi entropies. In this work we present an alternative approach for computing the entanglement spectrum with quantum Monte Carlo for both continuum and lattice Hamiltonians. This method provides direct access to the matrix elements of the spatially reduced density matrix and we determine an estimator that can be used in variational Monte Carlo as well as other Monte Carlo methods. The algorithm is based on using a generalization of the Swap operator, which can be extended to calculate a general class of density matrices that can include combinations of spin, space, particle and even momentum coordinates. We demonstrate the method by applying it to the Hydrogen and Nitrogen molecules and describe for the first time how the spatial entanglement spectrum encodes a covalent bond that includes all the many body correlations.