I know how to use least square for estimating a constant value given a bunch of measurements. It is

The linear least squares estimation formula you consider implies that the observation $z\in <!--\; \in \; -->{\mathbb{R}}^m$ and the parameter of interest $x\in <!--\; \in \; -->{\mathbb{R}}^n$, ($n\le <!--\; \le \; -->m$), are linearly related as
$z=Hx+n(1)$
where $H\in <!--\; \in \; -->{\mathbb{R}}^{m\times <!--\; \times \; -->n}$, and $n\in <!--\; \in \; -->{\mathbb{R}}^m$ is a noise vector with uncorellated elements of zero mean and equal variance.
For the sensor problem, since the sensors provide the cartesian coordinates of the point (if I understood correct), it is convenient to consider as parameter of interest the cartesian coordinates $x\in <!--\; \in \; -->{\mathbb{R}}^2$ of the point. Noting that the $i$-th sensor provides a measurement $z_i=x+n_i\in <!--\; \in \; -->{\mathbb{R}}^2$, stacking all, say, $N$ measurements to a single observation vector $z\triangleq <!--\; \triangleq \; -->[z_1^T,z_2^T,\dots <!--\; \dots \; -->,z_N^T{]}^T\in <!--\; \in \; -->{\mathbb{R}}^{2N}$ results in the linear model of (1). I will leave it to you to determine the form of $H\in <!--\; \in \; -->{\mathbb{R}}^{2N\times <!--\; \times \; -->2}$. You may then proceed in determining the LS estimate of $x$, from which you can then obtain an estimate of the polar coordinates by the well-known transformation formulas.
Note that the above approach estimates the polar coordinates implicitly (as a by-product of the cartesian coordinates estimate). The reason is that the polar coordinates are not linearly related to the observations (i.e., for $x=[r,\mathrm{\theta <!--\; \theta \; -->}]$, you cannot find a matrix $H$ so that (1) holds). If you wish to obtain a direct estimate of the polar coordinates from the observations you would have to perform a non-linear parameter estimation procedure.
EDIT: Since your question considers only $N=1$ sensor, the estimate of the cartesian coordinates is simply $z=(z_1,z_2)$ and the corresponding polar coordinates estimate is $(\stackrel{^<!--\; ^\; -->}{r},\stackrel{^<!--\; ^\; -->}{\mathrm{\theta <!--\; \theta \; -->}})=(\sqrt{z_1^2+z_2^2},{\tan }^{-<!--\; -\; -->1}<!--\; \; -->(z_2/z_1))$