The density matrix of a single qubit can be obtained from the density matriz of the whole system through the partial trace operation:
\begin{align}
\rho_i = \text{tr}_{\neq i}(\rho).
\end{align}
In order to do this operation, you must first find a basis $\{\vert \alpha\rangle\}$ for the Hilbert space without the qubit. Then the partial trace can be written as
\begin{align}
\text{tr}_{\neq i}(\rho) = \sum_{\{\vert \alpha\rangle\}} \langle \alpha\vert \rho \vert \alpha \rangle
\end{align}
To calculate this, you should just notice that every ket or bra of the basis of the full Hilbert space can be written as a tensor product of a ket (bra) of the qubit and a ket (bra) of the rest of the system: $\vert n_1, n_2... n_i ... n_N \rangle = \vert n_1... n_{i-1}, n_{i+1},...n_N\rangle \otimes \vert n_i \rangle $.
After doing this you obtain $\rho_i$, which is a $2\times2$ matrix, as you said.
As for your final question, no, in general you cannot obtain $\vert\psi\rangle$ from the $\rho_i$s, because the set of these matrices contains less information than the full $\rho$. We can see this from the following example:
Consider the pure states
\begin{align}
\vert + \rangle &= \frac{1}{\sqrt{2}} \left(\vert \uparrow \downarrow \rangle + \vert \downarrow \uparrow \rangle \right) \\
\vert - \rangle &= \frac{1}{\sqrt{2}} \left( \vert \uparrow \downarrow \rangle - \vert \downarrow \uparrow \rangle \right).
\end{align}
The density matrices which describes these pure states are then
\begin{align}
\rho^{\pm} = \vert \pm \rangle \langle \pm \vert
\end{align}
The reduced density matrices are given by
\begin{align}
\rho^{\pm}_1
&= {}_2\langle \uparrow \vert \rho^{\pm} \vert \uparrow \rangle_{2} + {}_2\langle \downarrow \vert \rho^{\pm} \vert \downarrow \rangle_{2} \\
&= \frac{1}{2} \left( (\pm \vert \downarrow \rangle_{1}) (\pm {}_{1} \langle \downarrow \vert) + (\vert \uparrow \rangle_{1}) ({}_{1} \langle \uparrow \vert) \right) \\
&= \frac{1}{2} \left( \vert \downarrow \rangle \langle \downarrow \vert_{1} + \vert \uparrow \rangle \langle \uparrow \vert_{1} \right) \\
\end{align}
and
\begin{align}
\rho^{\pm}_2
&= {}_1\langle \uparrow \vert \rho^{\pm} \vert \uparrow \rangle_{1} + {}_1\langle \downarrow \vert \rho^{\pm} \vert \downarrow \rangle_{1} \\
&= \frac{1}{2} \left( (\vert \downarrow \rangle_{2}) ({}_{2} \langle \downarrow \vert) + ( \pm\vert \uparrow \rangle_{2}) ( \pm{}_{2} \langle \uparrow \vert) \right) \\
&= \frac{1}{2} \left( \vert \downarrow \rangle \langle \downarrow \vert_{2} + \vert \uparrow \rangle \langle \uparrow \vert_{2} \right) \\
\end{align}.
So you can see that although $\rho^{+}$ and $\rho^{-}$ represent distinct bipartite states, their reduced density matrices are the same. This means that with only the reduced density matrices we cannot tell which bigger state they came from, $\rho^{+}$ or $\rho^{-}$.
