Why is the "canonical momentum" for the Dirac equation not defined in terms of the "gauge covariant derivative"? The canonical momentum is always used to add an EM field to the Schrödinger/Pauli/Dirac equations. Why does one not use the gauge covariant derivative? As far as I can see, the difference is a factor i in front of the vector potential. I know I'm combining two seemingly unrelated things, but they seem very similar, an the covariant form seems much "better" with respect to the inherent gauge freedom in the EM field. I can also see that with the canonical momentum form, the equations remain unchanged after an EM and a QM (phase) gauge transformation. Suffice to say my field theory knowledge is not that impressive.
 A: The identification goes as follows:
$$ \text{Kin. Mom.}~=~ \text{Can. Mom.} ~-~\text{Charge} \times \text{Gauge Pot.} $$
$$ \updownarrow $$
$$  m\hat{v}_{\mu} ~=~ \hat{p}_{\mu} - qA_{\mu}(\hat{x})$$
$$ \updownarrow $$
$$ \frac{\hbar}{i} D_{\mu} ~=~ \frac{\hbar}{i}\partial_{\mu} - qA_{\mu}(x) $$
$$ \updownarrow $$
$$  D_{\mu} ~=~ \partial_{\mu} -\frac{i}{\hbar} qA_{\mu}(x) $$
$$ \updownarrow $$
$$ \text{Cov. Der.}~=~ \text{Par. Der.} ~-~\frac{i}{\hbar}\text{Charge} \times \text{Gauge Pot.} $$
The imaginary unit $i$ is needed, e.g. because the derivative is an anti-hermitian operator (recall the usual integration-by-part proof), while the  momentum is required to be a hermitian operator in quantum mechanics.
A: (Qmechanic has already given the answer. However, since I spent some time writing the answer below so I am anyway posting it.)
Consider charged particle with charge $q$, and (nonzero) rest mass $m$ moving in a spacetime with coordinates $(x^0,x^1,\ldots ,x^{n-1})$. When there is no electromagnetic field then the action of particle is given as 
$$\tag {1} S=-mc\int \sqrt {\eta_{\mu\nu}\dot x^\mu(\lambda)\dot x^\nu(\lambda)}\,\mathrm d\lambda$$ 
Where $\lambda$ is the parameter along the trajectory $x^{\mu}(\lambda)$ of the particle and $\dot x^{\mu}$ means $\partial x^{\mu}(\lambda)/\partial\lambda$. Note that the action is Lorentz invariant. When there is a $U(1)$ gauge field $-iA_{\mu}\,\mathrm dx^{\mu}$ then we can add one more Lorentz invariant term to this action to generalize it as :-
$$\tag{2}S=-mc\int \sqrt {\eta_{\mu\nu}\dot x^\mu (\lambda)\dot x^\nu (\lambda)}\,\mathrm d\lambda-(q/c)\int \eta_{\mu\nu}A^{\mu}\dot x^{\mu}(\lambda)\,\mathrm d\lambda$$
Now in order to proceed its convenient to work in a particular inertial frame, and look at things from the viewpoint of inertial observer corresponding to that frame. In such a frame we can take $x^0(\lambda)/c=t=\lambda$. Above integral was from some point $\lambda_{0}$ to some $\lambda_{1}$. Now it becomes an integral from $t_0=x^0(\lambda_0)/c$ to $t_1=x^0(\lambda_1)/c$ and can be written as 
$$\tag{3} S=-mc^2\int \sqrt {1-v^2/c^2}dt-\int \left(qA^{0}-\frac {q}{c} \sum_{i}v^{i}.A^{i}\right) \mathrm dt$$
So 
$$\tag{4} L=-mc^2\sqrt {1-v^2/c^2}- \left(qA^{0}- \frac {q}{c} \sum_{i}v^{i}.A^{i}\right)$$
canonical momentum corresponding to $x^i$ can now be obtained as partial derivative of $L$ wrt $v^i=\mathrm dx^i/\mathrm dt$ and is given as :-
$$\tag{5} \pi_i=mv^i/\sqrt {1-v^2/c^2}+\frac{q}{c}A^{i}$$
Thus, as Qmechanic has answered, canonical momentum corresponding to $i$ th coordinate is physical momentum along that coordinate plus a contribution from gauge potential. 
Even without choosing a particular inertial frame we could find the canonical momentum $\pi_{\mu}$ corresponding to $x^\mu$ by taking the derivative of $L$ in its covariant form wrt $\dot x^\mu$. This would give -
$$\tag{6}\pi_{\mu}=-mc\eta_{\mu\nu}\dot x^{\nu}/\sqrt {\eta_{\alpha\beta}\dot x^\alpha\dot x^\beta}-\frac{q}{c}A_{\mu}$$
Now in classical mechanics above equation is nothing but a map from velocity space to phase space. Its only when we move to QM that we represent canonical momentums as derivatives wrt spatial coordinates. Again for convenience lets work in a particular inertial frame. Here momentum conjugate to $x^i$ is $\pi_i$ as given by the equation [5]. So as usual in QM we quantize by requiring the corresponding operators to satisfy
$$\tag{7}[ X^i,\Pi_j]=i\delta^i_{j}\hbar$$
We can represent this algebra on Hilbert space of functions on space $R^{n-1}$ (note that  spacetime is $R^n$) by defining $X^i$ to be the operator which acts as multiplication by $x^i$, and $\Pi_i$ to be the operator which acts as the derivative $-i\hbar\partial/\partial x^i$. From equation [5] we see that
$$\tag{8}\textrm{mechanical momentum operator}=-i\hbar\partial/\partial x^i-\frac{q}{c}A^{i}$$
Or taking $-i\hbar$ common we get 
$$\tag{9}\textrm{mechanical momentum operator}=-i\hbar\left(\partial/\partial x^i-i\frac{q}{c\hbar}A^{i}\right)=-i\hbar D_i$$
Where $D_i=\partial/\partial x^i-i\displaystyle\frac{q}{c\hbar}A^{i}$
