The use of a highly localized plasmonic field has enabled us to achieve sub-nanometer resolution of Raman images for single molecules. The inhomogeneous spatial distribution of plasmonic field has become an important factor that controls the interaction between the light and the molecule. We present here a gauge invariant interaction Hamiltonian (GIIH) to take into account the non-uniformity of the electromagnetic field distribution in the non-relativistic regime. The theory has been implemented for both resonant and nonresonant Raman processes within the sum-over-state framework. It removes the gauge origin dependence in the phenomenologically modified interaction Hamiltonian (PMIH) employed in previous studies. Our calculations show that, in most resonant cases, the Raman images from GIIH are similar to those from PMIH when the origin is set to the nuclear charge center of the molecule. In the case of nonresonant Raman images, distinct differences can be found from two different approaches, while GIIH calculations provide more details and phase information of the images. Furthermore, the results from GIIH calculations are more stable with respect to the computational parameters. Our results not only help to correctly simulate the resonant and nonresonant Raman images of single molecules but also lay the foundation for developing gauge invariant theory for other linear and nonlinear optical processes under the excitation of non-uniform electromagnetic field.