We construct the $S$-wave part of the electromagnetic vector annihilation current to $\mathcal{O}({\ensuremath{\alpha}}_{s}{v}^{2})$ on the lattice for heavy quarks whose dynamics are described by the NRQCD action, and $v$ is the nonrelativistic quark velocity inside the meson. The lattice vector current for $Q\overline{Q}$ annihilation is expressed as a linear combination of lattice operators with quantum numbers $L=0$, ${J}^{P}={1}^{\ensuremath{-}}$, and the coefficients are determined by matching this lattice current to the corresponding continuum current in QCD to $O({v}^{2})$ at one-loop. The annihilation channel gives a complex amplitude and a proper choice for the contours of integration is needed; a simple Wick rotation is not possible. In this way, and with a careful choice of subtraction functions in the numerical integration, the Coulomb-exchange and infrared singularities appearing in the amplitudes are successfully treated. The matching coefficients are given as a function of the heavy quark mass $Ma$ in lattice units. An automated vertex generation program written in Python is employed, allowing us to use a realistic NRQCD action and an improved gluon lattice action. A change in the definition of either action is easily accommodated in this procedure. The final result, when combined with lattice simulation results, describes the electromagnetic decays of heavy quarkonia, notably the $\ensuremath{\Upsilon}$ meson.