We compute a light front wave function for heavy vector mesons based on long distance matrix elements constrained by decay width analyses in the Non Relativistic QCD framework. Our approach provides a systematic expansion of the wave function in quark velocity. The first relativistic correction included in our calculation is found to be significant, and crucial for a good description of the HERA exclusive $\mathrm{J}/\psi$ production data. When looking at cross section ratios between nuclear and proton targets, the wave function dependence does not cancel out exactly. In particular the fully non-relativistic limit is found not to be a reliable approximation even in this ratio. The important role of the Melosh rotation to express the rest frame wave function on the light front is illustrated.