[Proj] Computing geodesics with Matlab; accuracy of spherical distances.

Charles Karney charles.karney at sri.com
Sun Dec 15 11:40:40 EST 2013


About a year ago, I released a native Matlab implementation of the
geodesic algorithms in GeographicLib.  The standalone package is
available here:

   http://www.mathworks.com/matlabcentral/fileexchange/39108

This works with Octave and Matlab.  The solution is vectorized; solving
1000000 inverse problems in parallel takes a few seconds on my home
machine.

This capability allow you to extend Tobler's 1964 paper

   http://dx.doi.org/10.1111/j.0033-0124.1964.009_q.x

comparing spherical and ellipsoidal distances.  For example, a rather
accurate way of estimating the ellipsoidal distance is to transfer the
problem to a sphere of radius (2*a+b)/3 with the spherical latitude set
to

     atan((1-f)^(3/2)*tan(phi))

(which is roughly the rectifying latitude).  The RMS relative error in
the distance is 0.05% and the max relative error is 0.11% (WGS84
ellipsoid).  The following short Matlab program illustrates
this result

     ell=defaultellipsoid;                   % WGS84
     m=1000000;                              % number of points
     % sample points uniformly
     lat1=asind(2*rand(m,1)-1); lon1=180*(2*rand(m,1)-1);
     lat2=asind(2*rand(m,1)-1); lon2=180*(2*rand(m,1)-1);
     % the ellipsoidal distance
     sell=geoddistance(lat1,lon1,lat2,lon2,ell);
     % the mean sphere
     a=ell(1); f=ecc2flat(ell(2)); b=a*(1-f); sph=[(2*a+b)/3,0];
     p=1.5;                                  % convert the latitude
     mu1=atand((1-f)^p*tand(lat1)); mu2=atand((1-f)^p*tand(lat2));
     % the spherical distance
     ssph=geoddistance(mu1,lon1,mu2,lon2,sph);
     err=(ssph-sell)./sell;                  % the relative error
     100*[sqrt(mean(err(:).^2)),max(abs(err(:)))]

If the geographic latitude is used (set p=0 above), the errors are RMS =
0.18%, max = 0.56%.  (Note that Tobler gives too optimistic an estimate
of the error because he only uses 27 pairs of points in his tests,
instead of 1000000 used here.)

-- 
Charles Karney <charles.karney at sri.com>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662


More information about the Proj mailing list