We explore here the benefits of using constraints from seismic tomography in gravity data inversion and how inverted density distributions can be improved by doing so. The methodology is applied to a real field case in which we reconstruct the density structure of the Pyrenees along a southwest–northeast transect going from the Ebro basin in Spain to the Arzacq basin in France. We recover the distribution of densities by inverting gravity anomalies under constraints coming from seismic tomography. We initiate the inversion from a prior density model obtained by scaling a pre-existing compressional seismic velocity Vp model using a Nafe–Drake relationship: the Vp model resulting from a full-waveform inversion of teleseismic data. Gravity data inversions enforce structural similarities between Vp and density by minimizing the norm of the cross-gradient between the density and Vp models. We also compare models obtained from 2.5-D and 3-D inversions. Our results demonstrate that structural constraints allow us to better recover the density contrasts close to the surface and at depth, without degrading the gravity data misfit. The final density model provides valuable information on the geological structures and on the thermal state and composition of the western region of the Pyrenean lithosphere.