sgo_distance {sgo}

Calculate distance(s) between points

Description

Calculates the distance between OS National Grid Reference points. Points with angular coordinates will use the Harvesine or Vicenty formulae.

Usage

sgo_distance(x, y, by.element = FALSE,
  which = ifelse(isTRUE(x$epsg==27700 || x$epsg==7405), "BNG", "Vicenty"),
  grid.true.distance = ifelse(isTRUE(x$epsg==27700 || x$epsg==7405),
  TRUE, FALSE), iterations = 100L)

Arguments

x

A sgo_points object describing a set of points in a geodetic coordinate system.

y

A sgo_points object, defaults to x.

by.element

Logical variable. If TRUE, return a vector with distance between the first elements of x and y, the second, etc. If FALSE, return the dense matrix with all pairwise distances.

which

Character vector. For geodetic coordinates one of Harvesine or Vicenty. It defaults to BNG for points in 'OS British National Grid' coordinates.

grid.true.distance

Logical variable. Currently only used for BNG coordinates. If TRUE it returns the true (geodesic) distance.

iterations

Numeric variable. Maximum number of iterations used in the Vicenty method.

Details

This function can use two different methods when working with geodetic coordinates: When which = "Vicenty" the Vincenty's formula is used to calculate the geodesics (distance) on an ellipsoid to an accuracy of up to a millimetre. If such accuracy is not needed, which can also accept the string "Harvesine" which calculates the great-circle distance between two points on a sphere. Harvesines are faster to compute than the Vicenty distances but can result in an error of up to 0.5%.

When working with (BNG) planar coordinates the Local Scale Factor is the scale distortion inherent in the map projection at a point. When grid.true.distance is TRUE the function computes a line scale factor using Simpson's Rule to achieve greater accuracy and approximate the distance to the true geodesic distance. When it is FALSE the Euclidean distance in the plane is calculated.

Note: Considering F as the scale factor, we have that S (True distance) = s (grid distance) / F
For most purposes the scale factor can be taken as constant for distances up to 20km (errors not exceeding 1 or 2 parts er million) and equal to the mid point of the line. For longer lines, this routine computes a scale factor for both ends (F1 and F2) and the mid point (Fm) and then uses the Simpson's Rule:
F = 1/6(F1 + 4Fm + F2)

Value

If by.element is FALSE sgo_distance returns a dense numeric matrix of dimension length(x) by length(y). Otherwise it returns a numeric vector of length x or y, the shorter one being recycled. Distances involving empty geometries are NA. All distances are returned in metres.

References

Thaddeus Vincenty, 1975. Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations. Survey Review, 23:176, 88-93, DOI: 10.1179/sre.1975.23.176.88

Examples

p1 <- sgo_points(list(-3.9369, 56.1165), epsg=4326)
lon <- c(-4.25181,-3.18827)
lat <- c(55.86424, 55.95325)
pts <- sgo_points(list(longitude=lon, latitude=lat), epsg=4326)
p1.to.pts <- sgo_distance(p1, pts, by.element = TRUE)

## Perimeter of a polygon defined as a series of ordered points:
lon <- c(-6.43698696, -6.43166843, -6.42706831, -6.42102546,
-6.42248238, -6.42639092, -6.42998435, -6.43321409)
lat <- c(58.21740316, 58.21930597, 58.22014035, 58.22034112,
58.21849188, 58.21853606, 58.21824033, 58.21748949)
pol <- sgo_points(list(lon, lat), epsg=4326)
# Create a copy of the polygon with its coordinates shifted one position
# so that we can calculate easily the distance between vertices
coords <- sgo_coordinates(pol)
pol.shift.one <- sgo_points(rbind(coords[-1, ], coords[1, ]), epsg=pol$epsg)
perimeter <- sum(sgo_distance(pol, pol.shift.one, by.element=TRUE))