Re: Lat/Long Distance Calculation Syntax

[Richard Ristow]

At 11:56 AM 2/25/2010, Adam B. Troy wrote:

I'm trying to refine a syntax to calculate the distance between a pair of latitude/longitudes.  The syntax is from Ray Levesque's site, but as the distance between points increases, the error of the calculation greatly increases,

The logic of the syntax you posted seems to be the same as the logic I posted about a year ago(1), and that seemed to work against at least a small test data set. Here's a comparison run. I've changed the logic you posted in only one way: dividing the calculated distance by .999670, corresponding to an Earth radius of 3,959.873 statute miles. In the test data, Providence to Chicago distance is from Jon Peck's posting "Re: Function for arc cosine", Thu, 7 Jun 2007 09:53:06 -0500;  others test points are compared with distance calculated at site http://www.movable-type.co.uk/scripts/latlong.html.

On the face of it, it's working. What lat-lon pairs give you values that seem wrong? How do those values compare with the results form the above Web site, or other independent calculation?

*  ...............   Initialize constants      ................. .

*  The following code requires that angles in the base system    .
*  (SPSS) be in radians, so that the trigonometric distance is   .
*  in radians, and can be multiplied directly by the Earth's     .
*  radius.                                                       .

DO IF $CASENUM EQ 1.
*  These initializations >MUST< be performed:                    .
*  #EarthRad is the Earth's radius in whatever units you please; .
*  the calculated distance will be in those units:               .
*              6,372.7976    km,                                 .
*              3,959.873     statute miles,                      .
*              3,441.035     nautical miles.                     .
.  COMPUTE
   #EarthRad =  3959.873  /* statute miles */.

*  #AngleCvt is the number of your angle units (degrees, here)   .
*  in one of SPSS's angle units (radians). It uses that          .
*  ARCTAN(1)is PI/4 radians or (in any angle measure) 1/8 circle .
.
.  COMPUTE
   #AngleCvt =  360  /* Number of input units in a full circle */
              /(8*ARTAN(1)).

END IF.

*  ...............   Compute distance          ................  .
*  Compute distance between points with coordinates              .
*  (lat1,lon1) and (lat2,lon2)                                   .

compute distance = #EarthRad*
    (2*artan(1)-arsin(      sin(lat1/#AngleCvt)   /* (sin(lat1)  */
                           *sin(lat2/#AngleCvt)   /* .sin(lat2)  */
                         +  cos(lat1/#AngleCvt)   /* +cos(lat1)  */
                           *cos(lat2/#AngleCvt)   /* .cos(lat2)  */
                           *cos(lon2/#AngleCvt    /* .cos(long2  */
                               -lon1/#AngleCvt)   /*     -long1))*/
                      )).

FORMAT    lat1 lat2 lon1 lon2 (F5.2).
FORMAT    Distance GivenDist  (F6.2).



*  From Levesque website: ... .

compute theta = lon2 - lon1.
compute d1 = (sin((artan(1)/45)*(lat2))
             *sin((artan(1)/45)*(lat1)))
           + (cos((artan(1)/45)*(lat2))
            * cos((artan(1)/45)*(lat1))
            * cos((artan(1)/45)*(theta))).

if range(d1,-1,1) d2 = 2*artan(1) - arsin(d1) .
compute d3 = (45/artan(1))*(d2).
compute distmile = d3 * 60 * 1.1515 / .999670  .

FORMAT theta    (F5.2)
       d1 d2 d3 (F4.3)
       distmile (F6.2).

LIST.

List
|-----------------------------|---------------------------|
|Output Created               |26-FEB-2010 00:05:36       |
|-----------------------------|---------------------------|
City             City             GivenD
1     lat1  lon1 2     lat2  lon2 ist    distance theta   d1   d2   d3 distmile

Pvd  41.90 87.65 Chi  41.73 71.43 836.27  834.33  -16.2 .978 .211 12.1  834.33
A1   42.00 80.00 A2   39.00 70.00 564.33  564.52  -10.0 .990 .143 8.17  564.52
B1   44.00 70.00 B2   49.00 85.00 791.00  791.03  15.00 .980 .200 11.4  791.03

Number of cases read:  3    Number of cases listed:  3
=============================
APPENDIX: Test data, and code
=============================
DATA LIST LIST /
   City1    lat1    lon1 City2    lat2    lon2   GivenDist
   (A4,     F6.2,   F6.2,A4,      F6.2,   F6.2,  F7.2).
BEGIN DATA
   Pvd     41.90   87.65 Chi     41.73   71.43   836.27
   A1      42.00   80.00 A2      39.00   70.00   564.33
   B1      44.00   70.00 B2      49.00   85.00   791.00
END DATA.

.  /*--  LIST /*-*/.


*  ...............   Initialize constants      ................. .

*  The following code requires that angles in the base system    .
*  (SPSS) be in radians, so that the trigonometric distance is   .
*  in radians, and can be multiplied directly by the Earth's     .
*  radius.                                                       .

DO IF $CASENUM EQ 1.
*  These initializations >MUST< be performed:                    .
*  #EarthRad is the Earth's radius in whatever units you please; .
*  the calculated distance will be in those units:               .
*              6,372.7976    km,                                 .
*              3,959.873     statute miles,                      .
*              3,441.035     nautical miles.                     .
.  COMPUTE
   #EarthRad =  3959.873  /* statute miles */.

*  #AngleCvt is the number of your angle units (degrees, here)   .
*  in one of SPSS's angle units (radians). It uses that          .
*  ARCTAN(1)is PI/4 radians or (in any angle measure) 1/8 circle .
.
.  COMPUTE
   #AngleCvt =  360  /* Number of input units in a full circle */
              /(8*ARTAN(1)).

END IF.

*  ...............   Compute distance          ................  .
*  Compute distance between points with coordinates              .
*  (lat1,lon1) and (lat2,lon2)                                   .

compute distance = #EarthRad*
    (2*artan(1)-arsin(      sin(lat1/#AngleCvt)   /* (sin(lat1)  */
                           *sin(lat2/#AngleCvt)   /* .sin(lat2)  */
                         +  cos(lat1/#AngleCvt)   /* +cos(lat1)  */
                           *cos(lat2/#AngleCvt)   /* .cos(lat2)  */
                           *cos(lon2/#AngleCvt    /* .cos(long2  */
                               -lon1/#AngleCvt)   /*     -long1))*/
                      )).

FORMAT    lat1 lat2 lon1 lon2 (F5.2).
FORMAT    Distance GivenDist  (F6.2).



*  From Levesque website: ... .

compute theta = lon2 - lon1.
compute d1 = (sin((artan(1)/45)*(lat2))
             *sin((artan(1)/45)*(lat1)))
           + (cos((artan(1)/45)*(lat2))
            * cos((artan(1)/45)*(lat1))
            * cos((artan(1)/45)*(theta))).

if range(d1,-1,1) d2 = 2*artan(1) - arsin(d1) .
compute d3 = (45/artan(1))*(d2).
compute distmile = d3 * 60 * 1.1515 / .999670  .

FORMAT theta    (F5.2)
       d1 d2 d3 (F4.3)
       distmile (F6.2).

LIST.

-----------------------------
(1)Date: Sun, 29 Mar 2009 00:57:59 -0400
From:    Richard Ristow <[hidden email]>
Subject: Re: Calculating distance between latitude and longitude
To:      [hidden email]

===================== To manage your subscription to SPSSX-L, send a message to [hidden email] (not to SPSSX-L), with no body text except the command. To leave the list, send the command SIGNOFF SPSSX-L For a list of commands to manage subscriptions, send the command INFO REFCARD