Gravity and magnetic imaging produce numerical images proportional to the density and magnetization source distribution through the upward continuation of potential field data. However, the inherent ambiguity and non-uniqueness in gravity or magnetic imaging restrict the reliability of the imaging model. We develop a joint iterative imaging framework for gravity and magnetic data based on Gramian constraints. The imaging fields of gravity and magnetic field are initially calculated by the depth from extreme points method. With the recovered magnetization and density distribution, the gradient directions of the Gramian function for these model parameters are calculated and utilized to regularize the joint imaging directions for gravity and magnetic fields. The introduction of Gramian constraints achieves mutual complementarity of the gravity and magnetic data and generates density and magnetization models with more distinct boundaries and improved structural coherence. The proposed approach is tested by two synthetic examples and applied to field data from the Galinge iron deposit in Qinghai, China. The real case results are verified by information from drillholes and physical properties measurements of ore and rock samples. Joint imaging provides an alternative to joint inversion with improved resolution and reliability of imaging models compared to separate imaging through the complementarity of gravity and magnetic data.