We introduce a computational toolset, named GROmaρs, to obtain and compare time-averaged density maps from molecular dynamics (MD) simulations. GROmaρs efficiently computes density maps, by fast multi-Gaussian spreading of atomic densities onto a 3-dimensional grid. It complements existing map-based tools by enabling spatial inspection of atomic average localization during the simulations. Most importantly, it allows the comparison between computed and reference maps (e.g. experimental), through calculation of difference maps, and local and time-resolved global correlation. These comparison operations proved useful to quantitatively contrast perturbed and control simulation-data sets and to examine how much biomolecular systems resemble both synthetic and experimental density maps. This was especially advantageous for multi-molecule systems, in which standard comparisons, like RMSDs, are difficult to compute. In addition, GROmaρs incorporates absolute and relative spatial free energy estimates to provide an energetic picture of atomistic localization. This is an open-source GROMACS-based toolset, thus allowing for static or dynamic selection of atoms, or even coarse-grained beads, for the density calculation. Furthermore, masking of regions was implemented, to speed up calculations and to facilitate the comparison with experimental maps. Beyond map comparison, GROmaρs provides an straightforward method to detect solvent cavities and average charge distribution in biomolecular systems. We employed all these functionalities to inspect the localization of lipid and water molecules in aquaporin systems, the binding of cholesterol to the G Protein Coupled Chemokine Receptor Type 4, and the identification of permeation pathways through the dermicidin antimicrobial channel. Based on these examples, we anticipate a high applicability of GROmaρs for the analysis of MD simulations and their comparison with experimentally determined densities.