help geodist-------------------------------------------------------------------------------

Title

geodist-- Computes geodetic distances.

SyntaxIf one or more lat/lon coordinates are variables

geodistlat1 lon1 lat2 lon2[if] [in],generate(newvar)[options]

If all lat/lon coordinates are scalars or strings of numerical coordinates

geodistlat1 lon1 lat2 lon2,[options]

optionsDescription ------------------------------------------------------------------------- Mainmilesreport distances in milesEllipsoid

ellipsoid(#1,#2)alternative ellipsoid parameters(a,f)maxiter(#)maximum iterations#Sphere

spherecalculate great-circle distancesradius(#)custom radius#(in km) -------------------------------------------------------------------------

Description

geodistcomputes geodetic distances, i.e. the length of the shortest curve between two points along the surface of a mathematical model of the earth. By default, the input coordinates are assumed to be based on the WGS 1984 datum (the same used by Google Earth/Map and GPS devices) andgeodistcalculates ellipsoidal distances using Vincenty's (1975) equations.If the

sphereoption is selected,geodistuses a sphere with a radius of 6371km to approximate the shape of the earth and computes great-circle distances. This is significantly faster than computing ellipsoidal distances.Coordinates may be specified in any combination of variables, scalars, or strings of numerical coordinates. If one or more lat/lon coordinates are variables, then

generate(newvar)must be specified andnewvarwill be created to store the computed distances.Coordinates must be in signed decimal degrees, positive for north and east, and negative for south and west. Latitudes range from -90 to 90 and longitudes from -180 to 180.

+------+ ----+ Main +-------------------------------------------------------------Optionsmilesindicates that distances are to be reported in miles; if omitted, distances are in kilometers.+-----------+ ----+ Ellipsoid +--------------------------------------------------------

ellipsoid(#1,#2)is used to specify a reference ellipsoid other than the one used by default (WGS 1984). A reference ellipsoid is defined by two constants, (a,f). These are the length of the semi-major axis in meters (i.e. radius to the equator) and the flattening ratio. For example, the Airy 1830 reference ellipsoid can be specified withellipsoid(6377563.396,299.3249646).

maxiter(#)indicates the maximum number of iterations to perform when computing ellipsoidal distances. The default is 25. Note that virtually all cases require less than 10 iterations.+--------+ ----+ Sphere +-----------------------------------------------------------

sphererequests that great-circle distances be computed instead of ellipsoidal distances. By default, a radius of 6371km is used; this is the earth's mean radius (see http://en.wikipedia.org/wiki/Earth_radius#Mean_radii).

radius(#)specifies that great-circle distances are to be computed on a sphere with a radius of#km.

Accuracy and limitations

geodistuses the Vincenty (1975) equations when computing distances on a reference ellipsoid. According to Thomas and Featherstone (2005), these are thought to maintain submillimeter accuracy between all locations. However, the Vincenty algorithm cannot accurately calculate distances for near-antipodal points. In such cases,geodistwill replace the distance with the extended missing value.a. See http://geographiclib.sourceforge.net/cgi-bin/Geod for a calculator that does not have this limitation.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%.

Most longitude and latitude datasets are based on the WGS 1984 datum, which is what GPS devices and Google Earth/Map use. There are other datums and it is important to understand that distances will not be accurate if points are not based on the same datum. For example, here is a photo of my GPS lying on the Prime Meridian line at the Royal Observatory in Greenwich, England. It is part of the OSGB 1936 datum and is about 100 meters to the west of the Reference Meridian of WGS 1984 ( geodist 51.47803 -0.00145 51.47803 0).

When a sphere is used to approximate the shape of the earth,

geodistcomputes great-circle distances using two formulas that are equivalent and mathematically exact. However, each is susceptible to computational errors. As Sinnott (1984) pointed out, the cos function is incapable of distinguishing small angles. This makes the spherical law of cosines formula unsuitable for distances between points that are near each other. On the other hand, the haversine formula is incapable of distinguishing small differences at angles around sin(pi/2), which makes it unsuitable for near-antipodal points.geodistmakes the best of both worlds by using the haversine formula for almost all calculations and the law of cosines formula for near-antipodal points. This ensures submillimeter accuracy for all spherical distances computed bygeodist.

ExamplesCompute the distance between Ann Arbor, MI and Montreal, QC. The third line computes the distance on an ellipsoid with all flattening removed. Obviously, these are significantly different from driving distances, in this case 981km according to Google Maps.

.geodist 42.270873 -83.726329 45.545448 -73.639075.geodist 42.270873 -83.726329 45.545448 -73.639075, sphere.geodist 42.270873 -83.726329 45.545448 -73.639075, ellipsoid(6371000, > 1e20) With near anti-podal points, Vincenty's solution for distances along an ellipsoid will fail. Note that increasing the number of iterations may not makegeodistconverge:

.geodist 0 0 0 179.3.geodist 0 0 0 179.4.geodist 0 0 0 179.4, maxiter(1000) You can use the Karney (2011) calculator ifgeodistfails to return a distance.When computing spherical distances, the average earth radius is just that, and may not be the best to approximate ellipsoidal distances depending on the location of the points on earth. For example, at the poles, the radius of a sphere that matches the curvature of the earth is 6399.6km.

.geodist 88 0 88 180.geodist 88 0 88 180, radius(6399.6).geodist 88 0 88 180, sphereLet's illustrate the same point for locations distributed evenly over the State of Colorado. The northeast corner is at 41,-102 and the southwest corner is at 37,-109.

.clear.set obs 100.set seed 1234.gen double lat1 = 37 + (41 - 37) * uniform().gen double lon1 = -109 + (109 - 102) * uniform().gen double lat2 = 37 + (41 - 37) * uniform().gen double lon2 =-109 + (109 - 102) * uniform().sum lat1 lon1 lat2 lon2

.geodist lat1 lon1 lat2 lon2, gen(v).geodist lat1 lon1 lat2 lon2, gen(r6371) sphere.geodist lat1 lon1 lat2 lon2, gen(r6377) radius(6377).gen diff1 = abs(r6371-v).gen pcerr1 = 100 * diff1 / v.gen diff2 = abs(r6377-v).gen pcerr2 = 100 * diff2 / v.sumAny combination of variables and single lat/lon will work

.geodist lat1 lon1 41 -102, gen(v2).geodist lat1 lon1 lat2 -102, gen(v3).geodist lat1 -109 lat2 -102, gen(v4)If you want to compute the nearest neighbors, you can use cross to create all pairwise combinations of points and then sort.

.clear.set obs 100.gen id2 = 2000 + _n.gen double lat2 = 37 + (41 - 37) * uniform().gen double lon2 =-109 + (109 - 102) * uniform().tempfile f.save "`f'".clear.set obs 100.gen id1 = 1000 + _n.gen double lat1 = 37 + (41 - 37) * uniform().gen double lon1 = -109 + (109 - 102) * uniform().cross using "`f'".geodist lat1 lon1 lat2 lon2, gen(d).sort id1 d id2.by id1: keep if _n == 1 Alternatively, you can use geonear, available from SSC, to find the nearest neighbors:

.geonear id1 lat1 lon1 using "`f'", n(id2 lat2 lon2) ell

Saved results

geodistsaves the following inr()when no variable is specified:Scalars

r(iterations)number of iterations for ellipsoidal distancesr(distance)distance

References and acknowledgementsMany thanks to Chris Veness for the best web pages on the subject:

http://www.movable-type.co.uk/scripts/latlong-vincenty.html http://www.movable-type.co.uk/scripts/latlong.html The definition of the World Geodetic System 1984 is available from

http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf See Appendix A.1 for a list of reference ellipsoids.

C. F. F. Karney, Geodesics on an ellipsoid of revolution, Feb. 2011; preprint arxiv:1102.1215.

R. W. Sinnott, "Virtues of the Haversine", Sky and Telescope 68 (2), 159 (1984). Thanks to the University of Michigan's Shapiro Science Library.

C. M. Thomas and W. E. Featherstone, Validation of Vincenty's Formulas for the Geodesic Using a New Fourth-Order Extension of Kivioja's Formula, J. Surv. Engrg. Volume 131, Issue 1, pp. 20-26 (February 2005). Available from:

http://www.cage.curtin.edu.au/~will/thomas-featherstone.pdf

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

AuthorRobert Picard picard@netbox.com