How to calculate center point in geographic coordinates?

Question

How can I calculate center point of $N$ latitude/longitude pairs, if the distance is at the most $30m$? Assuming altitude is same for all points.

Answer 1

For completeness (I know this is pretty late; no need to change your accepted answer):

You have $n$ points on the globe, given as latitude $\phi_i$ and longitude $\lambda_i$ for $i=1,\ldots,n$ (adopting Wikipedia's notation). Consider a Cartesian coordinate system in which the Earth is a sphere centered at the origin, with $z$ pointing to the North pole and $x$ crossing the Equator at the $\lambda=0$ meridian. The 3D coordinates of the given points are $$\begin{align} x_i &= r\cos\phi_i\cos\lambda_i,\\ y_i &= r\cos\phi_i\sin\lambda_i,\\ z_i &= r\sin\phi_i, \end{align}$$ (compare spherical coordinates, which uses $\theta=90^\circ-\phi$ and $\varphi=\lambda$). The centroid of these points is of course $$(\bar x,\bar y,\bar z) = \frac1n \sum (x_i, y_i, z_i).$$ This will not in general lie on the unit sphere, but we don't need to actually project it to the unit sphere to determine the geographic coordinates its projection would have. We can simply observe that $$\begin{align} \sin\phi &= z/r, \\ \cos\phi &= \sqrt{x^2+y^2}/r, \\ \sin\lambda &= y/(r\cos\phi), \\ \cos\lambda &= x/(r\cos\phi), \end{align}$$ which implies, since $r$ and $r\cos\phi$ are nonnegative, that $$\begin{align} \bar\phi &= \operatorname{atan2}\left(\bar z, \sqrt{\bar x^2+\bar y^2}\right), \\ \bar\lambda &= \operatorname{atan2}(\bar y, \bar x). \end{align}$$

So yes, the code you linked to does appear to switch the latitude and longitude in the output. You should submit a patch to the author.

Answer 2

If your points are no more than 30m apart, then doing any trigonometric computations will likely introduce more errors than it avoids. So I'd say simply treat the coordinates as coordinates on a plane, and average them to get the centroid.

In order to avoid issues with the 180° meridian, you might want to pick one point as reference and ensure that all the others don't differ in longitude by more than 180°. If they do, then add or subtract 360° to remedy that. The end result might be outside the [-180°, 180°] range but can be adjusted by another 360° if desired.

Near the poles, the compted longitude will likely have a great deal of uncertainty due to input distributed over a wide range of values. But this only corresponds to the fact that at the poles, large differences in longiude correspond to small differences in distance, so nothing wrong there. If you are even closer to the pole, there might be situations where the geodesics between your data points would be badly approximated by the planar interpretation. Roughly speaking, the computation would connect them along a parallel while the most direct connection might be a lot closer to the pole. But I'd expect such effects to only matter once you are within 100m or so of the pole, so probably they are not worth the effort, as even the badly computed result isn't completely off, like it would be for the 180° meridian case.