The problem you are describing (i.e. finding an approximation of some state given some number of identical copies of it and some set of measurements) is known as quantum state tomography or state tomography for short.
In practise, the most efficient schemes for state tomography will depend on a specific experiment's setup and limitations, for which different protocols exist (see the wikipedia page for an overview). Because the efficiency scaling of state tomography is notoriously bad in general, finding optimal tomography schemes for particular scenarios is an active area of research.
That being said, for the rest of this answer I will assume that efficiency scalings are not important for what you want to do and describe the general theory behind state tomography on a single qubit.
It can be shown that any density matrix can be written out as a linear combination of Pauli matrices, such that
$$
\rho = \frac{1}{2}(I + x\sigma_x + y\sigma_y + z\sigma_z) = \frac{1}{2} (I + \vec{r} \cdot \vec{\sigma})
$$
where $\vec{\sigma} = (\sigma_x, \sigma_y, \sigma_z)$ and $\vec{r} = (x, y, z)$ is the so-called Bloch vector representing the coordinates on the Bloch sphere you reference above. To find these, we observe that
$$
x = \textrm{Tr}(\rho \sigma_x) = \langle \sigma_x \rangle_\rho, \quad y = \textrm{Tr}(\rho \sigma_y) = \langle \sigma_y \rangle_\rho, \quad z = \textrm{Tr}(\rho \sigma_z) = \langle \sigma_z \rangle_\rho.
$$
and so each the amplitude of each component of the Bloch vector is given by the expectation value of it's associated Pauli operator on $\rho$. Recall that the expectation value $\langle A \rangle_\rho$ is given by the average over all eigenvalues returned by measurement of $A$ on $\rho$.
So, to find the state's Bloch vector, just perform each Pauli measurements many times to find the respective expectation value.
Finally, to retrieve the matrix representation of $\rho$, simply recall that
$$
I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}, \quad
\sigma_x = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \quad
\sigma_y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \quad
\sigma_z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix},
$$
and so
$$
\rho = \frac{1}{2} \begin{pmatrix} 1 + z & x - iy \\ x + iy & 1 - z \end{pmatrix}.
$$
Useful links and further reading:
Afterword:
By extending this protocol to two qubits we can flavour of why state tomography is not scalable in general. For two qubits $\rho_1$ and $\rho_2$ the combined state $\rho_{1,2} = \rho_1 \otimes \rho_2$ is given by
\begin{align}
\rho_{1,2} = \frac{1}{2} (I + \vec{r}_1 \cdot \vec{\sigma}_1) \otimes \frac{1}{2} (I + \vec{r}_2 \cdot \vec{\sigma}_2)
\end{align}
which is going to contain $4^2 = 16$ components, one for each $2$-fold Pauli operator. Clearly a strategy that requires $4^n$ different expectation values to be estimated to recover an $n$-qubit state is unfeasible for even small systems of qubits. So this is where the aforementioned research comes in.
