  We propose a novel approach to compute the cross section of near-threshold
$J/\psi$ and $\Upsilon$ production in electron-proton scattering at large
photon virtualities $Q^2$ based on an operator product expansion. We show that
the process can be used to extract the gluon part of the D-term gravitational
form factor of the proton. At the subleading level, it is also sensitive to the
trace anomaly effect of QCD.
