Re: [R] Converting coordinates to actual distances

From: Jim Lemon <bitwrit_at_ozemail.com.au>
Date: Fri 16 Sep 2005 - 17:35:27 EST

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 radians conversion as well.

Jim

#define APOSTROPHE 39
#define EARTH_RADIUS 6367

typedef struct {
  char *basename;
  char *number[2];
  double lat[2];
  double lon[2];
  char delimiter[4];

  char *data02;
  char *data03;
  char *data07;
  char *data08;
  FILE *input;
  FILE *output;

} PSTS_INFO; /* DegreesToRadians
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
*/

double DegreesToRadians(char *dms) {
  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@solstice.jpl.nasa.gov) Calculates 'great circle' distances using the Haversine formula, based upon a spherical earth of EARTH_RADIUS radius. 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[1] - psts_info->lon[0];   dlat = psts_info->lat[1] - psts_info->lat[0];   a = pow(sin(dlat / 2),2) + cos(psts_info->lat[0]) *    cos(psts_info->lat[1]) * 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));*/
  return(d * EARTH_RADIUS);
}



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Sep 21 02:47:15 2005

This archive was generated by hypermail 2.1.8 : Sun 23 Oct 2005 - 17:24:35 EST