Transferring between field and single-particle versions of the Dirac equation Schwartz's quantum field theory text moves between spinor fields and individual particles slickly, and I'm wondering what the justification is. To review, the Lagrangian 
$$\mathcal{L} = \overline{\psi} (i \gamma^\mu \partial_\mu - m)\psi$$ 
gives the Dirac equation, 
$$(i \gamma^\mu \partial_\mu - m)\psi = 0$$
which may be rearranged to the form $i \partial_t \psi = H_d \psi$ where
$$H_d = (i\partial_i \gamma_0 \gamma_i + \gamma_0 m).$$
Since this looks like normal QM time evolution, we can apparently interpret $H_d$ as the Hamiltonian, substituting in the momentum operator $p_i = -i\partial_i$. 
So though we started with a relativistic field equation, we now apparently have a nonrelativistic particle equation. The field $\psi$ is now a wavefunction. Field equations don't work anymore; we can start with $\mathcal{L}$ to get the field momentum $\pi = i\psi^\dagger$ and then get the Hamiltonian (density) $\mathcal{H}$, but this is totally different from $H_d$.


*

*How is this procedure legit? What implicit steps were taken to convert $\psi$ from a field to a wavefunction? Do they always work?

*How are the field momentum $\pi$ and the momentum operators $p_i$ related? Same for the Hamiltonian density $\mathcal{H}$ and Hamiltonian operator $H_d$.

 A: I am going to provide some further technical details, also in support of my original comments.
The Hamiltonian density corresponding to Lagrangian density $\mathcal L = \bar\psi (i\gamma^\mu\partial_\mu - m)\psi$ is 
$$
\mathcal H(\vec{x}, x^0) = \pi\dot\psi - \mathcal L = i\psi^\dagger\partial_0\psi - \psi^\dagger\gamma^0(i\gamma^0\partial_0 + i\gamma^i\partial_i - m)\psi = \bar\psi(-i\gamma^i\partial_i + m)\psi
$$ 
which gives the total Hamiltonian
$$
H = \int{d^3x \; \mathcal H(\vec{x}, x^0)} = \int{ d^3x \; \bar\psi(-i\gamma^i\partial_i + m)\psi} \equiv \int{ d^3x \; \psi^\dagger(-i\gamma^0\gamma^i\partial_i + m\;\gamma^0)\psi}
$$
Here $\psi$, $\psi^\dagger$ (and $\bar\psi = \psi^\dagger\gamma^0$) are understood to be field operators satisfying 
$$
\{\psi_\alpha(\vec{x},x^0), \psi_\beta(\vec{y},x^0)\} = \{\psi^\dagger_\alpha(\vec{x},x^0), \psi^\dagger_\beta(\vec{y},x^0)\} = 0\\
\{\psi_\alpha(\vec{x},x^0), \psi^\dagger_\beta(\vec{y},x^0)\} = \delta_{\alpha\beta}\delta(\vec{x} - \vec{y})
$$
From this in the Heisenberg representation one easily obtains
$$
i\partial_0\psi = [\psi, H] = (-i\gamma^0\gamma^i\partial_i + m\;\gamma^0)\psi \equiv (\gamma^0\gamma^i p_i + m\;\gamma^0)\psi
$$
Alternatively, the same form of the EOM follows from the original one derived from the Lagrangian, $(i\gamma^\mu\partial_\mu - m)\psi = 0$, upon multiplication by $\gamma^0$:
$$
(i\gamma^0\gamma^\mu\partial_\mu - m\; \gamma^0)\psi = 0 \;\;\Leftrightarrow \;\; i\partial_0\psi - (-i\gamma^0\gamma^i\partial_i + m\;\gamma^0)\psi = 0
$$
Note: There seems to be a sign discrepancy in the OP's form for $H_d$, but judging from the rest of the question the issue raised is not about a possible typo, but about the interpretation of the resulting Hamiltonian equation, and so I assume the latter is the case.
The Dirac field $\psi$ is defined on the associated Fock space. In order to make the transition to the Dirac equation for single particle states, recall that in the Heisenberg representation $\psi^\dagger(x)|0\rangle = |\vec{x}\rangle_H$. Given an arbitrary single-particle spinor $|\Psi\rangle_H$ consider now the EOM for the amplitude $_H\langle\vec{x}|\Psi\rangle_H = \langle 0|\psi(x)|\Psi\rangle_H \equiv \Psi(\vec{x}, x^0)$:
$$
i\partial_0 \Psi(\vec{x}, x^0) \equiv \langle 0| i\partial_0 \psi(x) |\Psi\rangle_H = \langle 0| (-i\gamma^0\gamma^i\partial_i + m\;\gamma^0)\psi(x) |\Psi\rangle_H \equiv \\
\equiv (-i\gamma^0\gamma^i\partial_i + m\;\gamma^0) \langle 0 |\psi(x)|\Psi\rangle_H = (-i\gamma^0\gamma^i\partial_i + m\;\gamma^0) _H\langle x |\Psi\rangle_H = (-i\gamma^0\gamma^i\partial_i + m\;\gamma^0) \Psi(\vec{x}, x^0)
$$
Here use was made of 
$$
_H\langle x |\Psi\rangle_H = \langle 0 | e^{iHx^0}\psi(\vec{x},x^0=0)e^{-iHx^0}|\Psi(x^0=0)\rangle = \\
= \langle 0 |\psi(\vec{x},x^0=0)|\Psi(x^0)\rangle  = \langle \vec{x}, x^0=0 | \Psi(x^0) \rangle = \Psi(\vec{x}, x^0)
$$
A similar approach can be used with many-particle states, the latter being defined by Slater determinants.
As for the relation between the canonical momentum $\pi$ and the momentum operator $\vec{p} = -i \nabla$, as giulio already noted in his answer, they are two different things: one is canonically conjugate to $\psi$, the other is the generator of translations in 3d-space. However please note the difference between the single-particle $\vec{p}$ and the total momentum $\bf P$. The latter acts on the Fock space. 
A: The Dirac field may be re-interpreted as a "second quantization" of a single particle theory with Hamiltonian $H_d$. There are two ways to do so. Firstly, the observation you made above: the field respects a Schroedinger-like equation with $H_d$ as Hamiltonian and therefore you may consider the state $\psi (x) = \Psi(x) |vacuum \rangle $ (I used $\Psi$ to denote the field, $\psi$ to denote the wave-function) as a genuine wave-function. Secondly, as remarked below, you may write the Hamiltonian of the Dirac field in a way that suggests a second-quantization formulation of a many-body operator: $\mathcal{H} = \int \bar{\psi}  H_d \psi$.
However, this interpretation poses a problem: the Hamiltonian $H_d$ is not positive definite. In fact, for every state of positive energy you also have one with the opposite energy. This problem may be overcome by requesting that the vacuum state is the one with all the negative energy states filled up. Because we are dealing with fermions, those state are effectively ruled out, and we are left with positive energy states. (This was the original interpretation given by Dirac, the "Dirac sea").
Now, this is a little bit artificious, and not necessary. In fact, as you already saw, if you follow the "Canonical Formalism" approach, starting from the Lagrangian you get to a different Hamiltonian $\mathcal{H}$ which is manifestly positive since is written as the sum of operators $a^{\dagger} a$ with positive coefficients (which are the energies). 
In this interpretation, you see the Dirac theory as a second quantization of states with definite quadri-momenta and spin-half, together with their "anti-particles", and you give up on any interpretation of $\psi$ as a wave-function. This is customary in QFT: you don't have a position operator and a wave-function $\psi(x)$. Instead, the observables you consider are just quadri-momenta, spin and helicities. 
Now, to explicitly address your questions:
1) I would say that this procedure is not completely legit, because it gets you to the negative energy states, a problem solved in artificious way with the Dirac sea. Moreover, the canonical formalism is a (more or less) automatic machinery that when fed with a Lagrangian will give you canonical variables written as sums of single-particle creators/annihilators and a positive Hamiltonian build with these variables. It always works, and does not requires you to put extra efforts in interpreting its results (like the Dirac sea), so why should we bother with other interpretations?
2) I guess the two are not strictly related. $p_i$ is a single particle operator, you may use it to define a translation operator in the fock space and I guess you obtain something like $P = \sum_p p a_p^{\dagger} a_p$. $\pi$ is just the variable canonically conjugate to $\psi$, it has a "x" dependance and its expansion on $a,a^{\dagger}$ does not contain their products. Hence the $p_i$ (and $P$) is not related to $pi$. For $H_d$ and $\mathcal{H}$, the two are related by $\mathcal{H} = \int \pi H_d \psi$. (See Weinberg, vol 1, pag 323).
