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

http://zaneselvans.org

+1 303 815 6866