# [R] Converting coordinates to actual distances

Jim Lemon bitwrit at ozemail.com.au
Fri Sep 16 09:35:27 CEST 2005

```Paul Brewin wrote:
> Hello,
>
> I've been searching for a method of converting Lat/Lon decimal
> coordinates into actual distances between points, and taking into
> account the curvature of the earth.  Is there such a package in R?  I've
> looked at the GeoR package, but this does not seem to contain what I am
> looking for.  Ideally the output would be a triangular matrix of
> distances.

Hi Paul,

Below is an implementation of the Haversine formula that I once had to
use for calculating distances between telephone exchanges. It's in C,
but could be translated to R fairly easily. I've included a degrees to

Jim

#define APOSTROPHE 39

typedef struct {
char *basename;
char *number;
double lat;
double lon;
char delimiter;
char *data02;
char *data03;
char *data07;
char *data08;
FILE *input;
FILE *output;
} PSTS_INFO;

Converts a string of the form nnndnn'nn"[NSEW] to a value in radians.
East (E) and South (S) are arbitrarily negative.  Note that this routine
will generate garbage if not passed a string of the correct form,
although it is unlikely to do too much damage.
Return value    value of the string in radians
-10     input not correct format
*/

int index = 0;
double degrees;
double minutes = 0;
double seconds = 0;
int mult;

degrees = atof(dms);
while(isdigit(dms[index])) index++;
if(dms[index++] == 'd') {
minutes = atof(&dms[index]);
minutes /= 60;
while(isdigit(dms[index])) index++;
if(dms[index++] == APOSTROPHE) {
seconds = atof(&dms[index]);
seconds /= 360;
while(isdigit(dms[index]) && dms[index]) index++;
if(dms[index++] == QUOTES) {
mult = (dms[index] == 'S' || dms[index] == 'E') ? -1 : 1;
degrees += minutes + seconds;
return(mult * degrees * M_PI/180);
}
}
}
return(-10);
}

/* GreatCircleDistance
(Formula and recommendations for calculation taken from:
http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
by Bob Chamberlain - rgc at solstice.jpl.nasa.gov)
Calculates 'great circle' distances using the Haversine formula,
This version extracts the two latitude/longitude pairs (as radians)
from the PSTS_INFO structure.
*/

double GreatCircleDist(PSTS_INFO *psts_info) {
double a,d,dlon,dlat;

dlon = psts_info->lon - psts_info->lon;
dlat = psts_info->lat - psts_info->lat;
a = pow(sin(dlat / 2),2) + cos(psts_info->lat) *
cos(psts_info->lat) * pow(sin(dlon / 2),2);
/* The penultimate result may be calculated either way - the first
is bulletproof, the second is a bit faster */
d = 2 * atan2(sqrt(a),sqrt(1 - a));
/* d = 2 * asin(sqrt(a));*/