5.1: What are some algorithms for calculating the distance between 2 points? Summaries : Use the following to calculate distance between 2 positions: float range(lat1,long1,lat2,long2) float lat1,long1,lat2,long2; /* * Given Arguments * lat1 (f) = latitude of position 1 (radians) * long1 (f) = longitude of position 1 (radians) * lat2 (f) = latitude of position 2 (radians) * long2 (f) = longitude of position 2 (radians) * Yielded Arguments * range (f) = range in radians */ { float a; float result; a = long1-long2; if (a < 0) a = -a; if (a > M_PI) a = 2. * M_PI - a; result = acos(sin(lat2)*sin(lat1)+cos(lat2)*cos(lat1)*cos(a)); return (result); } To convert the result from radians to nautical miles, multiply the result by 3437.74677. ---------------------------------------------------- /* ================================================================ /* Decimal Degree Distance Program /* ================================================================ /*N dd_distance.aml /* ---------------------------------------------------------------- /* Purpose: /* This program calculates great circle distance based on /* decimal degree coordinates supplied by the user. The output /* units are dependant upon the circumference units used. The /* formula used is extracted directly from ELEMENTS OF /* CARTOGRAPHY, 4e - Robinson, Sale, Morrison - John Wiley & /* Sons, Inc. - pp. 44-45. /*E /* ---------------------------------------------------------------- /* Parameters: /*A /* lat1,lon1 == First coordinate pair (in decimal degrees) /* lat2,lon2 == Second coordinate pair (in decimal degrees) /*E /* ---------------------------------------------------------------- /* Globals: /*X /* None. /*E /* ---------------------------------------------------------------- /* Functions Called: /*F /* None. /*E /* ---------------------------------------------------------------- /* Portability: /*O /* No system dependancies. /*E /* ---------------------------------------------------------------- /* History: /*H /* Lance Shipman October 5, 1990 Orignal AML coding /*E /* ================================================================*/ /* &args lat1, lon1, lat2, lon2 /* Set standard error routine... &severity &error &routine stderr /* /* Set radius, changing the units of radius will change the /* distance units... &s R 3958.754 /* statue miles /* Set PI &s PI 3.141592654 /* /* Convert degrees to radians &s a [angrad %lat1%] &s b [angrad %lat2%] &s p [angrad [abs [calc %lon1% - %lon2% ]]] /* /* Calculate greate circle distance in radians... &s dtheta [acos [calc ( [sin %a% ] * [ sin %b%] ) + ( [cos %a%] * ~ [cos %b%] * [cos %p%] ) ]] /* Convert the distance in radians to miles... &s .distance [calc ( [radang %dtheta%] * %PI% * %R% ) / 180 ] /* &ty %.distance% /* This line should be removed for program use... /* &return /* &routine stderr /* --------------------------------------------------- /* This is a standard error routine for general usage. /* Insert the severity line a the start of your AML, /* after the header and arge statement. /* /* Lance Shipman October 1990 @ESRI, Redlands /* --------------------------------------------------- /* &severity &error &ignore /* Ignore errors in this routine. &mess &on /* Turn on messages. &ty \\\\\ &ty Error in dd_distance.aml &ty &return;&return ----------------------------------------------------------------- Try Maling's book: Coordinate Systems and Map Projections. ----------------------------------------------------------------- Just so happens I ran across the following recently: > Angular distance between two places is given by > > cos (dist)=sin(lat1)sin(lat2)+cos(lat1)cos(lat2)cos(difference in long.) > > Then convert angle to distance by 60miles = 1 minute of angle. These are > nautical miles. > > Glen (gkm@cc.uow.edu.au ) --- From the longtitude and latitude of a point, you can easily calculate the components of a Cartesian vector originating at the centre of the sphere. For simplicity, assume for now that the sphere has radius 1: x1 = cos(longtitude1) * cos(latitude1) y1 = sin(longtitude1) * cos(latitude1) z1 = sin(latitude1) You have two points, so compute two such vectors, one for each point. Now take advantage of the fact that the scalar product of the two unit vectors is the cosine of the angle between them. [ Scalar product: x1*x2 + y1*y2 + z1*z2 ] On a sphere, the angle between those two vectors (if measured in radians) is conveniently also the shortest distance between the two points on the surface of the sphere (with radius 1). So take the inverse cosine and multiply by the radius, and you have the distance you were looking for. Anyway, here's the expression you want, simplified a little bit: cosdist = (cos(latitude1) * cos(latitude2) * cos(longtitude1 - longtitude2)) + (sin(latitude1) * sin(latitude2)) distance = radius * arccos(cosdist) Remember to use radians for all the longtitudes and latitudes, and make sure your arccos function always returns an angle (in radians) between 0 and pi. -------------------------------------------------------------------------- Computation of distances can be a complex subject, especially if you want precision and are thus dealing with an elliptical Earth. There are many algorithms out there and some of them make poor assumptions and thus do not produce accurate results for all cases. To compute geodesic ("Great Circle" or shortest) distances between two points you can use the "geod" program distributed with the proj system on charon.er.usgs.gov's (128.128.40.24) ftp anon (see pub/README for details on file needed). It will determine distances accurate to about a cm over 10,000km. There is a geod.1 file which explains usage.