Am I missing something here?
Well, you (among others) seem to be missing that to measure "momentum" is defined through the application of the gradient of the translation operator $\nabla \hat T_{\mathbf r}[~] := \frac{d}{d \mathbf r_{\mathcal S} }[~]$
to what's given through observational data (e.g. concerning a particular object $A$ under consideration)
by members of a suitable (inertial) system $\mathcal S$.
See also for instance this question: ("Momentum as Generator of Translations"; PSE/q/45067) , and questions linked there.
Arguably, sensibly, this operator should be applied to the quantity which directly characterizes the object $A$ under consideration: its duration (a.k.a. its "proper time") $\Delta \tau_A$ throughout the course of a (suitable, short) experimental trial.
Expressing explicitly $$ \Delta \tau_A := \sqrt{(\Delta \tau_{\mathcal S}[~A~])^2 - (\Delta \mathbf r_{\mathcal S}[~A~] / c)^2} $$
therefore $$ \eqalign{
\mathbf p_{\mathcal S}[~A~] & \simeq \frac{d}{d \mathbf r_{\mathcal S} }[~\sqrt{(\Delta \tau_{\mathcal S}[~A~])^2 - (\Delta \mathbf r_{\mathcal S}[~A~] / c)^2}~] \\
& = \frac{-\Delta \mathbf r_{\mathcal S}[~A~] / c^2}{\sqrt{(\Delta \tau_{\mathcal S}[~A~])^2 - (\Delta \mathbf r_{\mathcal S}[~A~] / c)^2}} \\
& = -\vec \beta_{\mathcal S}[~A~]~\gamma_{\mathcal S}[~A~] / c
}.$$
Certain proportionality coefficients can and should be factored in:
the "invariant mass" $m_A$, as an extensive quantity characterizing object $A$ (in the trial under consideration), and
$c^2$ for backward compatibility (in terms of dimensionality) with earlier, preliminary conceptions of "momentum" (e.g. in Newtonian physics), and
a minus sign.
Together, the proportionality coefficient $m_A~c^2$ constitutes the "invariant energy" ("rest energy", "center-of-mass energy") characterizing object $A$ (in the trial under consideration).
Consequently: $$ \mathbf p_{\mathcal S}[~A~] = m_A~c~\vec \beta_{\mathcal S}[~A~]~\gamma_{\mathcal S}[~A~] := m_A~\mathbf v_{\mathcal S}[~A~]~\gamma_{\mathcal S}[~A~].$$
Finally, expressing the relevant geometric data ($\Delta \tau_A$) together with suitable propotionality coefficients as a "phase" $\phi := \text{Exp}[~-i~\frac{m_A~c^2}{\hbar}~\Delta \tau_A~]$,
the appropriate operator to evaluate the momentum of object $A$ with respect to system $\mathcal S$ is: $$ \mathbf{ \hat p}_{\mathcal S} := -i~\hbar~\frac{d}{d \mathbf r_{\mathcal S} }[~].$$
Similarly, of course, for the evaluation of $A$'s energy with respect to system $\mathcal S$ by applying the operator $\hat E_{\mathcal S} := i~\hbar~\frac{d}{d \mathbf \tau_{\mathcal S} }[~]$:
$$ \eqalign{
E_{\mathcal S}[~A~] & := m_A~c^2~\frac{d}{d \mathbf \tau_{\mathcal S} }[~\sqrt{(\Delta \tau_{\mathcal S}[~A~])^2 - (\Delta \mathbf r_{\mathcal S}[~A~] / c)^2}~] \\ & = m_A~c^2~\frac{\Delta \tau_{\mathcal S}[~A~]}{\sqrt{(\Delta \tau_{\mathcal S}[~A~])^2 - (\Delta \mathbf r_{\mathcal S}[~A~] / c)^2}} \\ & = m_A~c^2~\gamma_{\mathcal S}[~A~] },$$
such that $A$'s "invariant energy" can also be expressed as
$$ m_A~c^2 = \sqrt{ (E_{\mathcal S}[~A~])^2 - c^2~(\mathbf p_{\mathcal S}[~A~])^2 }.$$