geodesics: Calculation of Geodesic Distances In guenardg/codep: Multiscale Codependence Analysis

 geodesics R Documentation

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"),
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 <guillaume.guenard@gmail.com>

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

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 15, 2023, 6:47 a.m.