geodesics: Calculation of Geodesic Distances

In guenardg/codep: Multiscale Codependence Analysis

Calculation of Geodesic Distances

Description

Function geodesics carries out the calculation of pairwise geodesic distances within a set of coordinates or between two sets thereof, using one of two calculation approaches.

Usage

geodesics(
  x,
  y,
  method = c("haversine", "Vincenty"),
  radius = 6371000,
  sma = 6378137,
  flat = 1/298.257223563,
  maxiter = 1024L,
  tol = .Machine$double.eps^0.75
)

Arguments

x	A set of geographic coordinates in the form of a two-column matrix or data.frame.
y	An other two-column matrix or data.frame containing an optional second set of coordinates.
method	The calculation method used to obtain the distances (default: haversine method; see details).
radius	Radius of the planetary body (when assuming a sphere; default: 6371000 m).
sma	Length of the semi-major axis of the planetary body (when assuming a revolution ellipsoid; default: 6378137 m).
flat	Flattening of the ellipsoid (default: 1/298.257223563).
maxiter	Maximum number of iterations, whenever iterative calculation is involved (default: 1024).
tol	Tolerance used when iterative calculation is involved (default: .Machine$double.eps^0.75; a machine dependent value).

Details

When only one set of coordinates is given to the function (i.e., when argument y is omitted), the function returns the pairwise distances in the form of a 'dist-class' object representing a lower-triangle matrix. When the second coordinate set is given, the function calculates the distances between each coordinate of argument x and each coordinate of argument y.

Two calculation methods are implemented. The first is the haversine formula, which assume the planetary body to be a sphere. The radius of that sphere is given to the function as its argument radius, with the default value being the mean radius of planet earth. Of the two methods implemented, the haversine formula is fastest but its precision depend on how well the planetary body match the sphericity assumption. The second method implemented is Vincenty's inverse formula, which assumes the the planetary body is a revolution ellipsoid, which is expected for rotating semi-fluid such as planet earth. Argument sma, the length of the semi-major axis, corresponds to the radius of the circle obtained when the revolution ellipsoid at the equator, whereas argument flat correspond to the compression of the sphere, along the diameter joining the poles, to form the ellipsoid of revolution. Their default values corresponds to parameters for planet Earth according to WGS84. These values, along with arguments maxiter and tol, are ignored when using the haversine formula, whereas the value of argument radius is ignored when using Vincenty's inverse formula.

Vincenty's inverse formula is more precise on planet Earth (on the order of 0.5mm) than the haversine formula, but it involves more computation time and may sometimes fail to converge. This is more likely for pairs of locations that are nearly antipodal or both (numerically) very close to the equator. The results returned by the function when using Vincenty's inverse formula are given a niter attribute that gives the number of iterations that were necessary to achieve convergence. Numbers greater than argument maxiter are indicative of failed convergence; a warning is issued in such a circumstance.

Geodesic distance matrices are non metric.

Value

A 'dist-class' object or, whenever argument y is provided, a matrix with as many rows as the number of rows in argument x and as many columns as the number of rows in argument y.

Author(s)

Guillaume Guenard and Pierre Legendre, Bertrand Pages Maintainer: Guillaume Guenard <[e-mail]>

References

Vincenty, T. 1975. Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations. Survey Review XXIII (176): 88-93 doi:10.1179/sre.1975.23.176.88

Inman, J. 1835. Navigation and Nautical Astronomy: For the Use of British Seamen (3 ed.). London, UK: W. Woodward, C. & J. Rivington

See Also

The dist-class and associated methods.

Examples

##
### First example: locations spread throughout the world
##
coords <- cbind(c(43,22,9,12,-40,72,-86,-22),
                c(-135,22,0,1,-45,12,27,-139))
res_hav <- geodesics(coords)  ## Default: the haversine formula
res_hav
res_vif <- geodesics(coords, method = "Vincenty")
res_vif
attr(res_vif,"niter") ## The numbers of iterations
res_vif-res_hav       ## Absolute difference
200*(res_vif-res_hav)/(res_vif+res_hav) ## Large relative difference
##
### Second example: locations nearer from one another
##
coords <- cbind(c(45.01,44.82,45.23,44.74),
                c(72.03,72.34,71.89,72.45))
res_hav <- geodesics(coords)
res_vif <- geodesics(coords, method = "Vincenty")
res_vif-res_hav       ## Absolute difference
200*(res_vif-res_hav)/(res_vif+res_hav) ## Relative difference are smaller
##

guenardg/codep documentation built on April 16, 2024, 9:01 p.m.