I finally ended up needing to implement great circle fitting, and went
ahead and implemented the suggestion below (yes, more than 8 months
later...). I don't think it's quite right though. The quantity to
minimize is u^T*A*u, but doing so doesn't make u one of the
eigenvectors of A, does it?
I generated a perfect great circle segment (as a series of lon,lat
points using spherical reckoning) and found its corresponding pole
(good_u) by going pi/2 radians away from it on the sphere,
perpendicular to its path.
I also converted the series of lon,lat points into x,y,z vectors (on
the unit sphere), and constructed A as:
In [552]: A = dot(array([x,y,z]),array([x,y,z]).transpose())
and then found the eigenvalues/vectors:
In [553]: eigvals_A, eigvecs_A = eig(A)
and none of them corresponds to the good_u which I found above, and
they don't minimize the product:
In [554]: [ dot(dot(eigvecs_A[N],A),eigvecs_A[N].transpose()) for N in
range(3) ]
Out[554]: [42.058431684800112, 25.787798426739176, 22.153769888460737]
whereas good_u does.
In [555]: dot(dot(good_u,A),good_u.transpose())
Out[555]: -5.9120729242764932e-15
Have I misunderstood something here? It's not immediately obvious to
me why choosing u such that it minimizes u^T * A * u would make it an
eigenvector of A. It has been a long time since I took linear algebra
though...
Thanks for any insights,
Zane
···
On Fri, Aug 8, 2008 at 10:35 AM, Charles R Harris <charlesr.harris@...287...> wrote:
On Mon, Aug 4, 2008 at 11:48 AM, Zane Selvans <zane@...1923...> wrote:
Does anyone out there happen to know a simple algorithm for least
squares fitting a great circle to a given set of lat/lon points on a
sphere? Seems like it might not be a crazy thing to add to the library.Depends on whether you need distance along the sphere surface or not. But if
you can deal with saggital distances it reduces to an eigenvalue problem.
Represent the great circle by a unit vector u perpindicular to the plane
that determined by the great circle, then minimize the sum\sum_{i=1}^{n} |dot(u,x_i)|^2
which you can rewrite by setting
A = \sum_{i=1}^{n} (x_i * x^T_i)
so that you end up minimizing
u^T * A * u
subject to the constraint u^T * u = 1. The vector u is then the eigenvector
corresponding to the smallest eigenvalue of A.Chuck
--
Zane A. Selvans
Amateur Earthling
+1 303 815 6866