Find the error in theta for noisy points in $R^2$
Let $(x,y)$ be Cartesian coordinates and $(r,\theta)$ be polar coordinates:
$$
\begin{cases} x = r\cos(\theta) \\ y = r\sin(\theta) \end{cases}
$$
Then we have the Lemma:
$$
dx^2 + dy^2 = dr^2 + (r \, d\theta)^2
$$
A very elementary proof of this will be given in the Appendix below.
In order to accomplish things, let's replace the infinitesimals by finite differences:
$$
(\Delta x)^2 + (\Delta y)^2 = (\Delta r)^2 + (r \cdot \Delta \theta)^2
$$
Due to homogeneity of space, we can write:
$$
\Delta x = \Delta y = N \quad \mbox{and} \quad \Delta r = r \cdot \Delta \theta
$$
Herewith we easily derive:
$$
2 N^2 = 2 (r \cdot \Delta \theta)^2 \quad \Longrightarrow \\
\Delta \theta =\frac{N}{r} = \frac{N}{\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}}
$$
Appendix
$$
\begin{cases}
dx = \partial x / \partial r\,dr + \partial x / \partial \theta\,d\theta \\
dy = \partial y / \partial r\,dr + \partial y / \partial \theta\,d\theta
\end{cases} \quad \Longrightarrow \\ \begin{cases}
dx = \cos(\theta) dr - r\sin(\theta) d\theta
= \cos(\theta) dr - \sin(\theta) \, r d\theta \\
dy = \sin(\theta) dr + r\cos(\theta) d\theta
= \sin(\theta) dr + \cos(\theta)\, r d\theta \end{cases}
\quad \Longrightarrow \\ \begin{cases}
(dx)^2 = \cos^2(\theta) (dr)^2 - 2 dr\cos(\theta)\sin(\theta)\,r d\theta + \sin^2(\theta)(r d\theta)^2 \\
(dy)^2 = \sin^2(\theta) (dr)^2 + 2 dr\cos(\theta)\sin(\theta)\,r d\theta + \cos^2(\theta)(r d\theta)^2 \end{cases}
$$
Adding the latter equations together proves the Lemma.
Upon request. The purpose of this second answer is to show that the OP's method
leads to the same result as in my first answer (approximately).
It is noticed that $\,arctan2\,$ is rather a programming expression, not a mathematical one.
And the $2$ is redundant because the angle $\theta$ is small. Whatever.
The derivation with $\,\arctan\,$ can be done as follows:
$$
d\arctan(y/x) = \frac{1}{1+(y/x)^2}d(y/x) =
\frac{1}{1+(y/x)^2}\left(\frac{1}{x}dy-\frac{y}{x^2}dx\right) \quad \Longrightarrow \\
d\theta = \frac{x\,dy-y\,dx}{x^2+y^2}
$$
In order to accomplish things, let's replace the infinitesimals by finite differences again:
$$
\Delta \theta = \frac{(x_2-x_1)\, \Delta y - (y_2-y_1)\, \Delta x}{(x_2-x_1)^2+(y_2-y_1)^2}
$$
The numerator of this fraction represents twice the area of a triangle,
which is determinant of a parallelogram spanned by the vectors $(\Delta x,\Delta y)$ and $(x_2-x_1,y_2-y_1)$.
While the denominator represents, of course, the square length of the vector $(x_2-x_1,y_2-y_1)$. For small angles,
the numerator can be written with the error $N$ as : $\;\left|(x_2-x_1,y_2-y_1)\right|\,N$ .
One could make a picture to see this; think of "twice" and the fact that the area of a circle is $\pi\,r^2$
while the perimeter is $2\pi\,r$ . Or make the parallelogram a rectangle. So the result is (approximately)
the same as in the first answer.