The answer is... it is not so simple. Some quantum mechanics follow, but the TL;DR version is that while $m_l=0$ corresponds to $p_z$, the orbitals for $m_l=+1$ and $m_l=-1$ lie in the $xy$-plane, but not on the axes. The reason for this outcome is that the wavefunctions are usually formulated in spherical coordinates to make the maths easier, but graphs in the Cartesian coordinates make more intuitive sense for humans. The $p_x$ and $p_y$ orbitals are constructed via a linear combination approach from radial and angular wavefunctions and converted into $xyz$. Thus, it is not possible to directly correlate the values of $m_l=\pm1$ with specific orbitals. The notion that we can do so is sometimes presented in introductory courses to make a complex mathematical model just a little bit simpler and more intuitive.
From Physical Chemistry by Atkins and DePaula, the three wavefunctions for $n=2$ and $l=1$ are as follows.
$$\begin{align}
&\Psi_{2,1,0}&&=r\cos{\theta}f(r)\\
&\Psi_{2,1,+1}&&=-\dfrac{r}{\sqrt{2}}\sin{\theta}\mathrm{e}^{\mathrm{i}\phi}f(r)\\
&\Psi_{2,1,-1}&&=\dfrac{r}{\sqrt{2}}\sin{\theta}\mathrm{e}^{-\mathrm{i}\phi}f(r)\\
&\end{align}$$
The notation is $\Psi_{n,l,m_l}$, $r$ is the radius, $\theta$ is the angle with respect to the $z$-axis and $\phi$ is the angle with respect to the $xz$-plane.
$$f(r)=\sqrt{\dfrac{Z^5}{32\pi a_0^5}}\mathrm{e}^{-Zr/2a_0}$$
in which $Z$ is the atomic number (or probably better nuclear charge) and $a_0$ is the Bohr radius.
In switching from spherical to Cartesian coordinates, we make the substitution $z=r\cos{\theta}$, so:
$$\Psi_{2,1,0}=zf(r)$$
This is $\Psi_{2p_z}$ since the value of $\Psi$ is dependent on $z$: when $z=0;\ \Psi=0$, which is expected since $z=0$ describes the $xy$-plane.
The other two wavefunctions are unhelpfully degenerate in the $xy$-plane. An equivalent statement is that these two orbitals do not lie on the $x$- and $y$-axes, but rather bisect them. Thus it is typical to take linear combinations of them to make the equation look prettier. Linear combinations are allowed by the maths of quantum mechanics. If any set of wavefunctions is a solution to the Schrödinger equation, then any set of linear combinations of these wavefunctions must also be a solution. We can do this because orbitals and the wavefunctions that describe them are not real physical objects. They constitute a mathematical model.
In the equations below, we're going to make use of some trigonometry, notably Euler's formula:
$$\mathrm{e}^{\mathrm{i}\phi}=\cos{\phi}+\mathrm{i}\sin{\phi}$$
$$\sin{\phi} = \frac{\mathrm{e}^{\mathrm{i}\phi}-\mathrm{e}^{-\mathrm{i}\phi}}{2\mathrm{i}}$$
$$\cos{\phi} = \frac{\mathrm{e}^{\mathrm{i}\phi}+\mathrm{e}^{-\mathrm{i}\phi}}{2}$$
We're also going to use $x=\sin{\theta}\cos{\phi}$ and $y=\sin{\theta}\sin{\phi}$.
$$\begin{align}
\Psi_{2p_x}=\frac{1}{\sqrt{2}}\left(\Psi_{2,1,+1}-\Psi_{2,1,-1}\right)=\frac{1}{2}\left(\mathrm{e}^{\mathrm{i}\phi}+\mathrm{e}^{-\mathrm{i}\phi} \right)r\sin{\theta}f(r)=r\sin{\theta}\cos{\phi}f(r)=xf(r) \\
\Psi_{2p_y}=\frac{\mathrm{i}}{\sqrt{2}}\left(\Psi_{2,1,+1}+\Psi_{2,1,-1}\right)=\frac{1}{2\mathrm{i}}\left(\mathrm{e}^{\mathrm{i}\phi}-\mathrm{e}^{-\mathrm{i}\phi} \right)r\sin{\theta}f(r)=r\sin{\theta}\sin{\phi}f(r)=yf(r)\\
\end{align}$$
So, while $m_l=0$ corresponds to $\Psi_{p_z}$, $m_l=+1$ and $m_l=-1$ cannot be directly assigned to $\Psi_{p_x}$ and $\Psi_{p_y}$. Rather $m_l=\pm1$ corresponds to $\{\Psi_{p_x},\Psi_{p_y} \}$. Put another way, I suppose we could say that $m_l=+1$ might correspond to $\Psi_{p_{x+y}}$ and $m_l=-1$ might correspond to $\Psi_{p_{x-y}}$.