Experience with geodist package

[jviana]

Hi everyone,
I am wondering if anyone has experience with the geodist package (http://ideas.repec.org/c/boc/bocode/s457147.html) for calculating geodetic distances. I have recently used it to calculate distances of 2 miles or less, however I am not familiar enough with geospatial data to think critically about how the calculations are done and how valid the results are. I am going to try to get my hands on ArcGIS to run the same procedure and compare results between the two methods, but in the meantime perhaps someone here can speak to their experiences and if they have compared the results generated from geodist with another "industry standard" like ArcGIS.

Thanks,
Joe

[Joe Canner]

Joe,

I can't speak to your specific question about comparison with ArcGIS, but a quick Google search yields quite a bit of information on geodetic distances, including sites (such as http://www.ga.gov.au/earth-monitorin...lgorithms.html) that provide several different methods and illustrate the relative accuracy of each. The different methods have different levels of complexity, but most would be relatively easy to program. Since ArcGIS is hard to set up and learn to use, if it were me, I would see what algorithm geodist uses, see whether its accuracy is sufficient and, if not, find a more accurate method and write my own program.

Regards,
Joe

[Brendan Halpin]

By coincidence, I have spent the evening coding and playing with the Haversine distance as described at http://www.movable-type.co.uk/scripts/latlong.html

The program below calculates distances between successive observations (with optionally, a lag greater than 1). Do

distance lat lon, gen(d0)

to create a distance variable, d0 (in kilometres), from latitude and longitude variables lat and lon.

I find it gives results that are almost indistinguishable from the "spherical law of cosines" method also described at that link.

Undoubtedly the person who wrote the geodist package has a better understanding of cartography that I do, but here is my code for the curiosity value:

program define distance
syntax varlist(min=2 max=2) , GEn(name) [LAg(real 1)]
local lat : word 1 of `varlist'
local lon : word 2 of `varlist'
local earthradius 6371
tempvar latrad lonrad hda

qui {
  gen double `latrad' = `lat'*_pi/180
  gen double `lonrad' = `lon'*_pi/180

  gen double `hda' = 0
  replace `hda' = sin((`latrad' - `latrad'[_n-`lag'])/2)^2 + cos(`latrad'[_n-`lag'])*cos(`latrad')*sin((`lonrad' - `lonrad'[_n-`lag'])/2)^2 if _n>1
  gen double `gen' = `earthradius'*2*atan2(sqrt(`hda'), sqrt(1-`hda'))
}
end

[Robert Picard]

I'm the author of geodist, which is a standalone version of the distance routines I wrote for geonear, a program that finds nearest neighbors using geodetic distances. Both are available from SSC.

As explained in detail in the help file, geodist calculates the length of the shortest curve between two points along the surface of a mathematical model of the earth. By default, the shape of the earth is approximated using the WGS 1984 reference ellipsoid (what GPS devices and Google Earth/Map use) and distances are computed using Vincenty's inverse solution equations. The Haversine and Laws of Cosine equations (which are mathematically equivalent) are used if the earth is modeled as a sphere.

I'm no expert on ArcGIS, but as far as I can tell, ArcGIS has historically calculated Euclidean distances between features. That's because ArcGIS works with projected coordinate systems. See http://blogs.esri.com/esri/arcgis/20...ting_geodesic/.

With respect to the validity of the results produced by my programs, I'm always happy to provide all certification scripts that I use upon request. Here's a small subset:

* Choose to run this script under the executable version
    
    version `c(stata_version)'
    cscript "geodist - Stata version `c(stata_version)' " adofile geodist
    

* =============================================================================
* Duplicate all 5 examples Vincenty used to compare his results with those of
* Rainsford. From Vincenty, T. (1975) Direct and inverse solutions of geodesics
* on the ellipsoid with application of nested equations, Survey Review 22(176):
* 88-93. Available from:
* http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf

    set obs 5
    
    // Rainsford line (a). Bessel ellipsoid. 
    gen double lat1 = -(33 + 26/60 + 0/3600) in 1
    gen double lon1 = 108 + 13/60 + 0/3600 in 1
    gen double lat2 = 55 + 45/60 in 1
    
    // Rainsford line (b), International ellipsoid
    replace lat1 = 26 + 7/60 + 42.83946/3600 in 2 
    replace lon1 = 41 + 28/60 + 35.50729/3600 in 2
    replace lat2 = 37 + 19/60 + 54.95367/3600 in 2
    
    // Rainsford line (c), International ellipsoid
    replace lat1 = 67 + 22/60 + 14.77638/3600 in 3
    replace lon1 = 137 + 47/60 + 28.31435/3600 in 3
    replace lat2 = 35 + 16/60 + 11.24862/3600 in 3
    
    // Rainsford line (d), International ellipsoid
    replace lat1 = -(0 + 59/60 + 53.83076/3600) in 4
    replace lon1 = 179 + 17/60 + 48.02997/3600 in 4
    replace lat2 = 1 in 4

    // Rainsford line (e), International ellipsoid
    replace lat1 = 1 + 1/60 + 15.18952/3600 in 5
    replace lon1 = 179 + 46/60 + 17.84244/3600 in 5
    replace lat2 = 1 in 5

    // forepoints in azimuths of 0
    gen double lon2 = 0
    
    geodist lat1 lon1 lat2 lon2 in 1,  e(6377397.155, 299.1528128) gen(d1)
    geodist lat1 lon1 lat2 lon2 in 2/5,  e(6378388, 297) gen(d2to5)
    
* Match Vincenty's delta S, Table II, col D
    
    assert round(d1[1]*1e6 - 14110526170,.1) == round(-0.4,.1)
    assert round(d2to5[2]*1e6 - 4085966703,.1) == round(-0.4,.1)
    assert round(d2to5[3]*1e6 - 8084823839,.1) == round(-0.7,.1)
    assert round(d2to5[4]*1e6 - 19960000000,.1) == round(-0.2,.1)
    assert round(d2to5[5]*1e6 - 19780006558,.1) == round(0.8,.1)
    

* =============================================================================
* Repeat using immediate form

    geodist `=-(33+26/60+0/3600)' `=108+13/60+0/3600' `=55+45/60' 0, ///
        e(6377397.155, 299.1528128)
    assert round(r(distance)*1e6 - 14110526170,.1) == round(-0.4,.1)
    
    geodist `=26+7/60+42.83946/3600' `=41+28/60+35.50729/3600' `=37+19/60+54.95367/3600' 0, ///
        e(6378388, 297)
    assert round(r(distance)*1e6 - 4085966703,.1) == round(-0.4,.1)
    
    geodist `=67+22/60+14.77638/3600' `=137+47/60+28.31435/3600' `=35+16/60+11.24862/3600' 0, ///
        e(6378388, 297)
    assert round(r(distance)*1e6 - 8084823839,.1) == round(-0.7,.1)
    
    geodist `=-(0+59/60+53.83076/3600)' `=179+17/60+48.02997/3600' 1 0, ///
        e(6378388, 297)
    assert round(r(distance)*1e6 - 19960000000,.1) == round(-0.2,.1)
    
    geodist `=1+1/60+15.18952/3600' `=179+46/60+17.84244/3600' 1 0, ///
        e(6378388, 297)
    assert round(r(distance)*1e6 - 19780006558,.1) == round(0.8,.1)
    
    
* =============================================================================
* Use Vincenty's equations to check spherical distances. Remove all flattening
* of the ellipsoid by using a very large flattening ratio. The default radius
* used by -geodist- for distances on a shpere is 6371km.

    geodist lat1 lon1 lat2 lon2, gen(dWGS84)    // just to compare
    geodist lat1 lon1 lat2 lon2, gen(dv) e(6371000,1e99)
    geodist lat1 lon1 lat2 lon2, gen(ds) sphere
    assert abs(dv-ds) < 1e-9
    gen delta = dWGS84 - ds
    sum d*
    

* =============================================================================
* From http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp

    // FlindersPeak to Buninyong using GRS80 ellipsoid
    geodist `=-(37 + 57/60 + 03.72030/3600)' `=144 + 25/60 + 29.52440/3600' ///
        `=-(37 + 39/60 + 10.15610/3600)' `=143 + 55/60 + 35.38390/3600', ///
        e(6378137,298.257222101)
        
    // match results in meters
    assert round(r(distance)*1e3 ,.3) == round(54972.271,.3)

Both ellipsoid and spherical distances calculated by geodist are accurate to machine precision (to the millimeter). However, as is pointed out in the help file:

Obviously, most real-life data points will not be located on the reference ellipsoid. Since the reference ellipsoid is designed to approximate mean sea level (or more precisely the geoid), expect distances at sea level to be particularly accurate. However, when computing distances on Lake Titicaca, at 3.2km above the reference ellipsoid, expect an error in proportion to the increase in radius of the earth (6374.2/6371) or 0.05%.

I also suspect that in many cases, driving distances would be more useful than "as the crow flies" distances. According to Google Maps, the distance between Ann Arbor, MI and Montréal, QC is 981km. The distance according to geodist is 887km. Travel time might be even more relevant.

[Jasper Tjaden]

Dear Robert,

Thank you very much for geodist. I have a question:
I keep getting an error message after geodist var1 var 2 var 3 var4, gen(x) ("latitude 1 must be between -90 and 90"). All my coordinates are in the right range. I had several attempts with 'format' and with gen, floor() to change but geodist does still not work.

Do you have any idea?
Many thanks.

[Robert Picard]

I suggest that you use variable names that are more representative of their content. How did you check that the coordinates are within the right range. What does

summarize var1

report?.

Note that you will get this error message if var1 is missing for all observations.