## matplotlib-users

 [Matplotlib-users] least squares fitting of great circle to points on a sphere From: Zane Selvans - 2008-08-04 16:48:27 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. Thanks! Zane -- Zane Selvans Amateur Earthling http://zaneselvans.org zane@... 303/815-6866 PGP Key: 55E0815F 
 Re: [Matplotlib-users] least squares fitting of great circle to points on a sphere From: Zane Selvans - 2009-04-21 21:53:15 On Fri, Aug 8, 2008 at 10:35 AM, Charles R Harris wrote: > > > On Mon, Aug 4, 2008 at 11:48 AM, Zane Selvans 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 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 -- Zane A. Selvans Amateur Earthling http://zaneselvans.org +1 303 815 6866