On Fri, Aug 8, 2008 at 10:35 AM, Charles R Harris
<charlesr.har...@gmail.com> wrote:
>
>
> On Mon, Aug 4, 2008 at 11:48 AM, Zane Selvans <z...@ideotrope.org> 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
------------------------------------------------------------------------------
Stay on top of everything new and different, both inside and
around Java (TM) technology - register by April 22, and save
$200 on the JavaOne (SM) conference, June 2-5, 2009, San Francisco.
300 plus technical and hands-on sessions. Register today.
Use priority code J9JMT32. http://p.sf.net/sfu/p
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
