[GRASSLIST:5809] Re: geodesic distance
H Bowman
hamish_nospam at yahoo.com
Wed Mar 19 20:41:18 EST 2003
> is there a module available in GRASS to measure geodesic distances in
> a projected location (eqc in my case)?
FYI, there's:
d.geodesic (lat/lon only)
d.measure
d.rhumbline
but I don't think any of them are what you are really looking for.
For an algorithm, see:
http://mathworld.wolfram.com/GreatCircle.html
this might get you started:
data = [ lon1 lat1; lon2 lat2 ];
%lon (western hemisphere is 360 - lon!)
lam1 = data(1,1) * (2*pi/360); % convert to radians
lam2 = data(2,1) * (2*pi/360);
%lat
del1 = data(1,2) * (2*pi/360);
del2 = data(2,2) * (2*pi/360);
LatDist = ((data(1,2) - data(2,2) ) *111.191);
% convert to km (111.191km/degree lat)
r = 6367.444; % mean radius of the earth (km)
R1dotR2 = cos(del1) * cos(del2) * cos(lam1-lam2) + sin(del1) * sin(del2);
d = r * acos( R1dotR2 );
% mercator heading!!!!! A trip from 0 lon, 80N to 180 lon, 80N will
% tell you to go East not North!!! Only good locally!
theta = acos(LatDist/d);
Hamish
More information about the grass-user
mailing list