Nothing
R/sgo_distance.R
In sgo: Simple Geographical Operations (with OSGB36)
Defines functions
.direct.vicenty.ellipsoid
.inverse.vicenty.ellipsoid
.great.circle.harvesine
.local.scale.factor
.bng.distance
sgo_distance.sgo_points
#' @encoding UTF-8
#' @title 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.
#'
#' @name sgo_distance
#' @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)
#' @param x A \code{sgo_points} object describing a set of points in a geodetic
#' coordinate system.
#' @param y A \code{sgo_points} object, defaults to \code{x}.
#' @param by.element Logical variable. If \code{TRUE}, return a vector with
#' distance between the first elements of \code{x} and \code{y}, the second,
#' etc. If \code{FALSE}, return the dense matrix with all pairwise distances.
#' @param which Character vector. For geodetic coordinates one of
#' \code{Harvesine} or \code{Vicenty}. It defaults to \code{BNG} for points in
#' 'OS British National Grid' coordinates.
#' @param grid.true.distance Logical variable. Currently only used for BNG
#' coordinates. If \code{TRUE} it returns the true (geodesic) distance.
#' @param 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 \code{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, \code{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
#' \code{grid.true.distance} is \code{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
#' \code{FALSE} the Euclidean distance in the plane is calculated.
#'
#' \strong{Note}: Considering F as the scale factor, we have that
#' S (True distance) = s (grid distance) / F\cr
#' 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:\cr
#' F = 1/6(F1 + 4Fm + F2)
#'
#' @return
#' If \code{by.element} is \code{FALSE} \code{sgo_distance} returns a dense
#' numeric matrix of dimension length(x) by length(y). Otherwise it returns a
#' numeric vector of length \code{x} or \code{y}, the shorter one being
#' recycled. Distances involving empty geometries are \code{NA}.
#' All distances are returned in metres.
#' @references
#' Thaddeus Vincenty, 1975. \emph{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))
#' @export
sgo_distance <- function (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)
UseMethod("sgo_distance")
#' @export
sgo_distance.sgo_points <- function(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) {
if (missing(y))
y <- x
if (x$epsg != y$epsg)
stop("All points must have the same EPSG code")
if (isTRUE(x$epsg %in% c(4936, 3035, 4978, 3857)))
stop("This function doesn't support the input's EPSG")
default.simpson <- 20 #20 km
coords <- .sgo_points.2d.coords
if(isTRUE(x$epsg==27700 || x$epsg==7405)) {
p1 <- matrix(unlist(x[coords], use.names = FALSE), ncol = 2, byrow = FALSE)
p2 <- matrix(unlist(y[coords], use.names = FALSE), ncol = 2, byrow = FALSE)
if (by.element) {
.bng.distance(p1, p2, grid.true.distance, default.simpson)
} else {
rows.p1 <- nrow(p1)
rows.p2 <- nrow(p2)
# Kronecker product of our matrices with vectors of 1's:
# rep(1, nTimes) %x% mt would be 'similar' to rep(mt, times = nTimes)
# mt %x% rep(1, nTimes) would be 'similar' to rep(mt, each = nTimes)
m1 <- rep(1, rows.p2) %x% p1
m2 <- p2 %x% rep(1, rows.p1)
matrix(.bng.distance(m1, m2, grid.true.distance, default.simpson),
rows.p1, rows.p2)
}
} else {
p1 <- matrix(unlist(x[coords], use.names = FALSE),
ncol = 2, byrow = FALSE) / RAD.TO.DEG
p2 <- matrix(unlist(y[coords], use.names = FALSE),
ncol = 2, byrow = FALSE) / RAD.TO.DEG
if (by.element) {
if (which == "Harvesine") {
.great.circle.harvesine(p1, p2)
} else if (which == "Vicenty") {
.inverse.vicenty.ellipsoid(p1, p2, x$datum, iterations)$distance
}
} else {
rows.p1 <- nrow(p1)
rows.p2 <- nrow(p2)
m1 <- rep(1, rows.p2) %x% p1
m2 <- p2 %x% rep(1, rows.p1)
if (which == "Harvesine") {
matrix(.great.circle.harvesine(m1, m2), rows.p1, rows.p2)
} else if (which == "Vicenty") {
matrix(.inverse.vicenty.ellipsoid(m1, m2, x$datum, iterations)$distance,
rows.p1, rows.p2)
}
}
}
}
#' @noRd
#' @param p1 A matrix of coordinates
#' @param p2 A matrix of coordinates
#' @param grid.true.distance Logical value
#' @param dist.simpson: distance (in km) greater than this will apply Simpson's
#' Rule when calculating true (geodesic) distances
.bng.distance <- function(p1, p2, grid.true.distance = TRUE,
dist.simpson = 20) {
E1 <- p1[, 1]; E2 <- p2[, 1]
N1 <- p1[, 2]; N2 <- p2[, 2]
dE <- E2 - E1
dN <- N2 - N1
s <- sqrt(dE * dE + dN * dN)
if (!grid.true.distance)
#return(round(s, 3)) # round to mm
return(s)
Em <- E1 + (E2 - E1) / 2 # E at midpoint
Nm <- N1 + (N2 - N1) / 2 # N at midpoint
if (any(s >= (dist.simpson * 1000))) {
Ft <- .local.scale.factor(c(E1, E2, Em), c(N1, N2, Nm))
len.E1 <- length(E1) # E1, E2, Em have the same length
F1 <- Ft[(1:len.E1)]
Fm <- Ft[((len.E1 + 1):(len.E1 * 2))]
F2 <- Ft[length(Ft) - ((len.E1-1):0)]
F <- (F1 + 4 * Fm + F2) / 6
} else {
#F for mid point only
F <- .local.scale.factor(Em, Nm)
}
#round(s / F, 3) # S (true distance)
s / F
}
#' @noRd
#' @param E Vector of Easting coordinates
#' @param N Vector of Northing coordinates
.local.scale.factor <- function(E, N) {
# ellipsoid parameters
params <- lonlat.ellipsoid[lonlat.ellipsoid$ellipsoid=="Airy1830",
c("a","b","e2")]
a <- params$a
b <- params$b
e2 <- params$e2
F0 <- 0.9996012717 #Central meridian scale factor
aF0 <- a * F0
bF0 <- b * F0
n <- (a-b) / (a+b)
N0 <- -100000; E0 <- 400000 # True origin
phi0 <- 49 / RAD.TO.DEG
# Initial latitude φ'
dN <- N - N0
phi <- phi0 + dN/aF0
M <- NA
phi.plus <- NA
phi.minus <- NA
repeat {
phi.minus <- phi - phi0
phi.plus <- phi + phi0
M <- bF0 * (
(1 + n * (1 + 5/4 * n * (1 + n))) * phi.minus
- 3 * n * (1 + n * (1 + 7 / 8 * n)) * sin(phi.minus) * cos(phi.plus)
+ (15 / 8 * n * (n * (1 + n))) * sin(2 * phi.minus) * cos(2 * phi.plus)
- 35 / 24 * n^3 * sin(3 * phi.minus) * cos(3 * phi.plus)
) # meridional arc
if ( max(abs(dN - M)) < 0.00001 ) { break } # ie until < 0.01mm
phi <- phi + (dN - M) / aF0
}
#cos.phi <- cos(phi)
sin.phi <- sin(phi)
splat <- 1 - e2 * sin.phi * sin.phi
sqrtsplat <- sqrt(splat)
# radius of curvature at latitude φ perpendicular to a meridian
nu <- aF0 / sqrtsplat
# radius of curvature of a meridian at latitude φ
rho <- aF0 * (1 - e2) / (splat * sqrtsplat)
# East - west component of the deviation of the vertical, squared
eta2 <- nu / rho - 1
XXI <- 1 / (2 * rho * nu)
XXII <- (1 + 4 * eta2 * eta2) / (24 * rho * rho * nu * nu)
dE <- E - E0
dE2 <- dE * dE
F0 * (1 + dE2 * XXI + dE2 * dE2 * XXII) # return F
}
#' @noRd
#' @param p1 A matrix of coordinates in radians
#' @param p2 A matrix of coordinates in radians
#' @param R Default defined by the International Union of Geodesy and Geophysics
.great.circle.harvesine <- function(p1, p2, R=6371008) {
hav.dlat <- sin((p2[, 2] - p1[, 2]) / 2)
hav.dlon <- sin((p2[, 1] - p1[, 1]) / 2)
h <- hav.dlat * hav.dlat + cos(p1[, 2]) * cos(p2[, 2]) * hav.dlon * hav.dlon
# https://en.wikipedia.org/wiki/Haversine_formula:
# When using these formulae, one must ensure that h does not exceed 1 due to
# a floating point error (d is only real for 0 ≤ h ≤ 1). h only approaches 1
# for antipodal points (on opposite sides of the sphere)—in this region,
# relatively large numerical errors tend to arise in the formula when finite
# precision is used.
h <- pmin(h, 1)
d <- 2 * R * asin(sqrt(h))
#round(d, 3) #round to mm (problaby shouldn't expect accuracy greater than m)
d
}
# Vicenty (inverse) iterative method to compute the geographical distance
# between two given points.
#' @noRd
#' @param p1 A matrix of coordinates in radians
#' @param p2 A matrix of coordinates in radians
#' @param datum A string containing "OSGB36", "WGS84" or "ETRS89"
#' @param iterations Scalar value. Number of iterations to reach convergence
.inverse.vicenty.ellipsoid <- function(p1, p2, datum, iterations = 100L) {
# ellipsoid parameters
ellipsoid <- lonlat.datum[lonlat.datum$datum==datum, "ellipsoid"]
params <- lonlat.ellipsoid[lonlat.ellipsoid$ellipsoid==ellipsoid,
c("a","b","f")]
a2 <- params$a * params$a
b <- params$b
b2 <- b * b
f <- 1 / params$f
one.f <- 1 - f
dlon <- (p2[, 1] - p1[, 1]) # difference of longitudes
tan.U1 <- one.f * tan(p1[, 2]) # tan(reduced latitude)
tan.U2 <- one.f * tan(p2[, 2])
cos.U1 <- 1 / sqrt(1 + tan.U1 * tan.U1)
cos.U2 <- 1 / sqrt(1 + tan.U2 * tan.U2)
sin.U1 <- tan.U1 * cos.U1
sin.U2 <- tan.U2 * cos.U2
sin.U1.U2 <- sin.U1 * sin.U2
cos.U1.U2 <- cos.U1 * cos.U2
lambda <- dlon
no.convergence <- FALSE
repeat {
sin.lambda <- sin(lambda)
cos.lambda <- cos(lambda)
sin.sigma.1 <- cos.U2 * sin.lambda
sin.sigma.2 <- cos.U1 * sin.U2 - sin.U1 * cos.U2 * cos.lambda
sin2.sigma <- sin.sigma.1 * sin.sigma.1 + sin.sigma.2 * sin.sigma.2
sin.sigma <- sqrt(sin2.sigma)
cos.sigma <- sin.U1.U2 + cos.U1.U2 * cos.lambda
sigma <- atan2(sin.sigma, cos.sigma)
sin.alpha <- cos.U1.U2 * sin.lambda / sin.sigma
coincident.points <- is.nan(sin.alpha) # or antipodal (from sin.sigma=0)
sin.alpha[coincident.points] <- 0
cos2.alpha <- 1 - sin.alpha * sin.alpha
cos.2sigma.m <- cos.sigma - 2 * sin.U1.U2 / cos2.alpha
cos.2sigma.m[is.nan(cos.2sigma.m)] <- 0 # cos2.alpha = 0 (equatorial line)
cos2.2sigma.m <- cos.2sigma.m * cos.2sigma.m
C <- f / 16 * cos2.alpha * (4 + f * (4 - 3 * cos2.alpha))
lambda.tmp <- lambda
lambda <- dlon + (1 - C) * f * sin.alpha *
(sigma + C * sin.sigma * (cos.2sigma.m + C * cos.sigma *
(-1 + 2 * cos2.2sigma.m)))
iterations <- iterations - 1L
if ( max(abs(lambda - lambda.tmp)) <= 1e-12 || iterations == 0L) {
break # ie until <= 0.06mm
}
}
if (iterations == 0L) {
no.convergence <- abs(lambda - lambda.tmp) > 1e-12
}
u2 <- cos2.alpha * (a2 - b2) / b2
A <- 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
B <- u2 / 1024 * (256 + u2 * (-128 + u2 *(74 - 47 * u2)))
delta.sigma <- B * sin.sigma * (cos.2sigma.m +
B/4 * (cos.sigma * (-1 + 2 * cos2.2sigma.m) -
B/6 * cos.2sigma.m * (-3 + 4 * sin.sigma * sin.sigma) *
(-3 + 4 * cos2.2sigma.m)))
# alpha1, alpha2: azimuths of the geodesic, clockwise from north
# alpha2 in the direction p1 p2 produced
alpha1 <- atan2(cos.U2 * sin.lambda,
cos.U1 * sin.U2 - sin.U1 * cos.U2 * cos.lambda)
alpha2 <- atan2(cos.U1 * sin.lambda,
-sin.U1 * cos.U2 + cos.U1 * sin.U2 * cos.lambda)
# From Chris Veness (https://www.movable-type.co.uk):
# Special handling of exactly antipodal points where sin2.sigma = 0
# (due to discontinuity atan2(0, 0) = 0 but atan2(NaN, 0) = PI/2 / 90°) - in
# which case bearing is always meridional, due north (or due south!)
alpha1[coincident.points] <- 0
alpha2[coincident.points] <- PI
#s <- round(unname(b * A * (sigma - delta.sigma)), 3) #round to mm
s <- unname(b * A * (sigma - delta.sigma))
if (any(no.convergence)) {
s[no.convergence] <- NA_real_
alpha1[no.convergence] <- NA_real_
alpha2[no.convergence] <- NA_real_
warning("Vicenty formula failed to converge. Try to increase iterations.")
}
list(distance=s, initial.bearing=alpha1, final.bearing=alpha2)
}
# Vicenty (direct) iterative method to compute a destination point from a given
# point and an initial bearing (both in radians). 's' (distance) in m.
#' @noRd
#' @param p1 A matrix of coordinates in radians
#' @param s A numeric vector. Distances (to compute the destination point)
#' @param alpha1 A numeric vector. Initial bearing in radians
#' @param datum A string containing "OSGB36", "WGS84" or "ETRS89"
#' @param iterations Scalar value. Number of iterations to reach convergence
.direct.vicenty.ellipsoid <- function(p1, s, alpha1, datum, iterations = 100L) {
# ellipsoid parameters
ellipsoid <- lonlat.datum[lonlat.datum$datum==datum, "ellipsoid"]
params <- lonlat.ellipsoid[lonlat.ellipsoid$ellipsoid==ellipsoid,
c("a","b","f")]
a2 <- params$a * params$a
b <- params$b
b2 <- b * b
f <- 1 / params$f
one.f <- 1 - f
tan.U1 <- one.f * tan(p1[, 2])
cos.U1 <- 1 / sqrt(1 + tan.U1 * tan.U1)
sin.U1 <- tan.U1 * cos.U1
sin.alpha1 <- sin(alpha1)
cos.alpha1 <- cos(alpha1)
sigma1 <- atan2(tan.U1, cos.alpha1)
sin.alpha <- cos.U1 * sin.alpha1
cos2.alpha <- 1 - sin.alpha * sin.alpha
u2 <- cos2.alpha * (a2 - b2) / b2
A <- 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
B <- u2 / 1024 * (256 + u2 * (-128 + u2 *(74 - 47 * u2)))
sigma <- s / (b * A)
no.convergence <- FALSE
repeat {
cos.2sigma.m <- cos(2 * sigma1 + sigma)
cos2.2sigma.m <- cos.2sigma.m * cos.2sigma.m
sin.sigma <- sin(sigma)
cos.sigma <- cos(sigma)
delta.sigma <- B * sin.sigma * (cos.2sigma.m + B/4 *
(cos.sigma * (-1 + 2 * cos2.2sigma.m) -
B/6 * cos.2sigma.m * (-3 + 4 * sin.sigma * sin.sigma) *
(-3 + 4 * cos2.2sigma.m)))
sigma.tmp <- sigma
sigma <- s / (b * A) + delta.sigma
iterations <- iterations - 1L
if ( max(abs(sigma - sigma.tmp)) <= 1e-12 || iterations == 0L) {
break
}
}
if (iterations == 0L) {
no.convergence <- abs(sigma - sigma.tmp) > 1e-12
}
x <- sin.U1 * sin.sigma - cos.U1 * cos.sigma * cos.alpha1
phi2 <- atan2(sin.U1 * cos.sigma + cos.U1 * sin.sigma * cos.alpha1,
one.f * sqrt(sin.alpha * sin.alpha + x * x))
lambda <- atan2(sin.sigma * sin.alpha1,
cos.U1 * cos.sigma - sin.U1 * sin.sigma * cos.alpha1)
C <- f / 16 * cos2.alpha * (4 + f * (4 - 3 * cos2.alpha))
L <- lambda - (1 - C) * f * sin.alpha *
(sigma + C * sin.sigma * (cos.2sigma.m + C * cos.sigma *
(-1 + 2 * cos2.2sigma.m)))
lambda2 <- p1[, 1] + L
alpha2 <- atan2(sin.alpha, -x)
#p2 <- cbind(x=lambda2, y=phi2)
if (any(no.convergence)) {
#p2[no.convergence, ]<- NA_real_
lambda2[no.convergence]<- NA_real_
phi2[no.convergence]<- NA_real_
alpha2[no.convergence] <- NA_real_
warning("Vicenty formula failed to converge. Try to increase iterations.")
}
# coincident points. Set alpha2 = alpha1
p2.is.p1 <- (s == 0)
if (any(p2.is.p1)) {
#p2[p2.is.p1, ] <- p1[p2.is.p1, ]
lambda2[p2.is.p1] <- p1[p2.is.p1, 1]
phi2[p2.is.p1] <- p1[p2.is.p1, 2]
alpha2[p2.is.p1] <- alpha1[p2.is.p1]
}
list(lon=lambda2, lat=phi2, final.bearing=alpha2)
}
Try the sgo package in your browser
Any scripts or data that you put into this service are public.
sgo documentation built on Sept. 23, 2022, 5:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .direct.vicenty.ellipsoid .inverse.vicenty.ellipsoid .great.circle.harvesine .local.scale.factor .bng.distance sgo_distance.sgo_points
#' @encoding UTF-8
#' @title 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.
#'
#' @name sgo_distance
#' @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)
#' @param x A \code{sgo_points} object describing a set of points in a geodetic
#' coordinate system.
#' @param y A \code{sgo_points} object, defaults to \code{x}.
#' @param by.element Logical variable. If \code{TRUE}, return a vector with
#' distance between the first elements of \code{x} and \code{y}, the second,
#' etc. If \code{FALSE}, return the dense matrix with all pairwise distances.
#' @param which Character vector. For geodetic coordinates one of
#' \code{Harvesine} or \code{Vicenty}. It defaults to \code{BNG} for points in
#' 'OS British National Grid' coordinates.
#' @param grid.true.distance Logical variable. Currently only used for BNG
#' coordinates. If \code{TRUE} it returns the true (geodesic) distance.
#' @param 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 \code{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, \code{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
#' \code{grid.true.distance} is \code{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
#' \code{FALSE} the Euclidean distance in the plane is calculated.
#'
#' \strong{Note}: Considering F as the scale factor, we have that
#' S (True distance) = s (grid distance) / F\cr
#' 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:\cr
#' F = 1/6(F1 + 4Fm + F2)
#'
#' @return
#' If \code{by.element} is \code{FALSE} \code{sgo_distance} returns a dense
#' numeric matrix of dimension length(x) by length(y). Otherwise it returns a
#' numeric vector of length \code{x} or \code{y}, the shorter one being
#' recycled. Distances involving empty geometries are \code{NA}.
#' All distances are returned in metres.
#' @references
#' Thaddeus Vincenty, 1975. \emph{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))
#' @export
sgo_distance <- function (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)
UseMethod("sgo_distance")
#' @export
sgo_distance.sgo_points <- function(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) {
if (missing(y))
y <- x
if (x$epsg != y$epsg)
stop("All points must have the same EPSG code")
if (isTRUE(x$epsg %in% c(4936, 3035, 4978, 3857)))
stop("This function doesn't support the input's EPSG")
default.simpson <- 20 #20 km
coords <- .sgo_points.2d.coords
if(isTRUE(x$epsg==27700 || x$epsg==7405)) {
p1 <- matrix(unlist(x[coords], use.names = FALSE), ncol = 2, byrow = FALSE)
p2 <- matrix(unlist(y[coords], use.names = FALSE), ncol = 2, byrow = FALSE)
if (by.element) {
.bng.distance(p1, p2, grid.true.distance, default.simpson)
} else {
rows.p1 <- nrow(p1)
rows.p2 <- nrow(p2)
# Kronecker product of our matrices with vectors of 1's:
# rep(1, nTimes) %x% mt would be 'similar' to rep(mt, times = nTimes)
# mt %x% rep(1, nTimes) would be 'similar' to rep(mt, each = nTimes)
m1 <- rep(1, rows.p2) %x% p1
m2 <- p2 %x% rep(1, rows.p1)
matrix(.bng.distance(m1, m2, grid.true.distance, default.simpson),
rows.p1, rows.p2)
}
} else {
p1 <- matrix(unlist(x[coords], use.names = FALSE),
ncol = 2, byrow = FALSE) / RAD.TO.DEG
p2 <- matrix(unlist(y[coords], use.names = FALSE),
ncol = 2, byrow = FALSE) / RAD.TO.DEG
if (by.element) {
if (which == "Harvesine") {
.great.circle.harvesine(p1, p2)
} else if (which == "Vicenty") {
.inverse.vicenty.ellipsoid(p1, p2, x$datum, iterations)$distance
}
} else {
rows.p1 <- nrow(p1)
rows.p2 <- nrow(p2)
m1 <- rep(1, rows.p2) %x% p1
m2 <- p2 %x% rep(1, rows.p1)
if (which == "Harvesine") {
matrix(.great.circle.harvesine(m1, m2), rows.p1, rows.p2)
} else if (which == "Vicenty") {
matrix(.inverse.vicenty.ellipsoid(m1, m2, x$datum, iterations)$distance,
rows.p1, rows.p2)
}
}
}
}
#' @noRd
#' @param p1 A matrix of coordinates
#' @param p2 A matrix of coordinates
#' @param grid.true.distance Logical value
#' @param dist.simpson: distance (in km) greater than this will apply Simpson's
#' Rule when calculating true (geodesic) distances
.bng.distance <- function(p1, p2, grid.true.distance = TRUE,
dist.simpson = 20) {
E1 <- p1[, 1]; E2 <- p2[, 1]
N1 <- p1[, 2]; N2 <- p2[, 2]
dE <- E2 - E1
dN <- N2 - N1
s <- sqrt(dE * dE + dN * dN)
if (!grid.true.distance)
#return(round(s, 3)) # round to mm
return(s)
Em <- E1 + (E2 - E1) / 2 # E at midpoint
Nm <- N1 + (N2 - N1) / 2 # N at midpoint
if (any(s >= (dist.simpson * 1000))) {
Ft <- .local.scale.factor(c(E1, E2, Em), c(N1, N2, Nm))
len.E1 <- length(E1) # E1, E2, Em have the same length
F1 <- Ft[(1:len.E1)]
Fm <- Ft[((len.E1 + 1):(len.E1 * 2))]
F2 <- Ft[length(Ft) - ((len.E1-1):0)]
F <- (F1 + 4 * Fm + F2) / 6
} else {
#F for mid point only
F <- .local.scale.factor(Em, Nm)
}
#round(s / F, 3) # S (true distance)
s / F
}
#' @noRd
#' @param E Vector of Easting coordinates
#' @param N Vector of Northing coordinates
.local.scale.factor <- function(E, N) {
# ellipsoid parameters
params <- lonlat.ellipsoid[lonlat.ellipsoid$ellipsoid=="Airy1830",
c("a","b","e2")]
a <- params$a
b <- params$b
e2 <- params$e2
F0 <- 0.9996012717 #Central meridian scale factor
aF0 <- a * F0
bF0 <- b * F0
n <- (a-b) / (a+b)
N0 <- -100000; E0 <- 400000 # True origin
phi0 <- 49 / RAD.TO.DEG
# Initial latitude φ'
dN <- N - N0
phi <- phi0 + dN/aF0
M <- NA
phi.plus <- NA
phi.minus <- NA
repeat {
phi.minus <- phi - phi0
phi.plus <- phi + phi0
M <- bF0 * (
(1 + n * (1 + 5/4 * n * (1 + n))) * phi.minus
- 3 * n * (1 + n * (1 + 7 / 8 * n)) * sin(phi.minus) * cos(phi.plus)
+ (15 / 8 * n * (n * (1 + n))) * sin(2 * phi.minus) * cos(2 * phi.plus)
- 35 / 24 * n^3 * sin(3 * phi.minus) * cos(3 * phi.plus)
) # meridional arc
if ( max(abs(dN - M)) < 0.00001 ) { break } # ie until < 0.01mm
phi <- phi + (dN - M) / aF0
}
#cos.phi <- cos(phi)
sin.phi <- sin(phi)
splat <- 1 - e2 * sin.phi * sin.phi
sqrtsplat <- sqrt(splat)
# radius of curvature at latitude φ perpendicular to a meridian
nu <- aF0 / sqrtsplat
# radius of curvature of a meridian at latitude φ
rho <- aF0 * (1 - e2) / (splat * sqrtsplat)
# East - west component of the deviation of the vertical, squared
eta2 <- nu / rho - 1
XXI <- 1 / (2 * rho * nu)
XXII <- (1 + 4 * eta2 * eta2) / (24 * rho * rho * nu * nu)
dE <- E - E0
dE2 <- dE * dE
F0 * (1 + dE2 * XXI + dE2 * dE2 * XXII) # return F
}
#' @noRd
#' @param p1 A matrix of coordinates in radians
#' @param p2 A matrix of coordinates in radians
#' @param R Default defined by the International Union of Geodesy and Geophysics
.great.circle.harvesine <- function(p1, p2, R=6371008) {
hav.dlat <- sin((p2[, 2] - p1[, 2]) / 2)
hav.dlon <- sin((p2[, 1] - p1[, 1]) / 2)
h <- hav.dlat * hav.dlat + cos(p1[, 2]) * cos(p2[, 2]) * hav.dlon * hav.dlon
# https://en.wikipedia.org/wiki/Haversine_formula:
# When using these formulae, one must ensure that h does not exceed 1 due to
# a floating point error (d is only real for 0 ≤ h ≤ 1). h only approaches 1
# for antipodal points (on opposite sides of the sphere)—in this region,
# relatively large numerical errors tend to arise in the formula when finite
# precision is used.
h <- pmin(h, 1)
d <- 2 * R * asin(sqrt(h))
#round(d, 3) #round to mm (problaby shouldn't expect accuracy greater than m)
d
}
# Vicenty (inverse) iterative method to compute the geographical distance
# between two given points.
#' @noRd
#' @param p1 A matrix of coordinates in radians
#' @param p2 A matrix of coordinates in radians
#' @param datum A string containing "OSGB36", "WGS84" or "ETRS89"
#' @param iterations Scalar value. Number of iterations to reach convergence
.inverse.vicenty.ellipsoid <- function(p1, p2, datum, iterations = 100L) {
# ellipsoid parameters
ellipsoid <- lonlat.datum[lonlat.datum$datum==datum, "ellipsoid"]
params <- lonlat.ellipsoid[lonlat.ellipsoid$ellipsoid==ellipsoid,
c("a","b","f")]
a2 <- params$a * params$a
b <- params$b
b2 <- b * b
f <- 1 / params$f
one.f <- 1 - f
dlon <- (p2[, 1] - p1[, 1]) # difference of longitudes
tan.U1 <- one.f * tan(p1[, 2]) # tan(reduced latitude)
tan.U2 <- one.f * tan(p2[, 2])
cos.U1 <- 1 / sqrt(1 + tan.U1 * tan.U1)
cos.U2 <- 1 / sqrt(1 + tan.U2 * tan.U2)
sin.U1 <- tan.U1 * cos.U1
sin.U2 <- tan.U2 * cos.U2
sin.U1.U2 <- sin.U1 * sin.U2
cos.U1.U2 <- cos.U1 * cos.U2
lambda <- dlon
no.convergence <- FALSE
repeat {
sin.lambda <- sin(lambda)
cos.lambda <- cos(lambda)
sin.sigma.1 <- cos.U2 * sin.lambda
sin.sigma.2 <- cos.U1 * sin.U2 - sin.U1 * cos.U2 * cos.lambda
sin2.sigma <- sin.sigma.1 * sin.sigma.1 + sin.sigma.2 * sin.sigma.2
sin.sigma <- sqrt(sin2.sigma)
cos.sigma <- sin.U1.U2 + cos.U1.U2 * cos.lambda
sigma <- atan2(sin.sigma, cos.sigma)
sin.alpha <- cos.U1.U2 * sin.lambda / sin.sigma
coincident.points <- is.nan(sin.alpha) # or antipodal (from sin.sigma=0)
sin.alpha[coincident.points] <- 0
cos2.alpha <- 1 - sin.alpha * sin.alpha
cos.2sigma.m <- cos.sigma - 2 * sin.U1.U2 / cos2.alpha
cos.2sigma.m[is.nan(cos.2sigma.m)] <- 0 # cos2.alpha = 0 (equatorial line)
cos2.2sigma.m <- cos.2sigma.m * cos.2sigma.m
C <- f / 16 * cos2.alpha * (4 + f * (4 - 3 * cos2.alpha))
lambda.tmp <- lambda
lambda <- dlon + (1 - C) * f * sin.alpha *
(sigma + C * sin.sigma * (cos.2sigma.m + C * cos.sigma *
(-1 + 2 * cos2.2sigma.m)))
iterations <- iterations - 1L
if ( max(abs(lambda - lambda.tmp)) <= 1e-12 || iterations == 0L) {
break # ie until <= 0.06mm
}
}
if (iterations == 0L) {
no.convergence <- abs(lambda - lambda.tmp) > 1e-12
}
u2 <- cos2.alpha * (a2 - b2) / b2
A <- 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
B <- u2 / 1024 * (256 + u2 * (-128 + u2 *(74 - 47 * u2)))
delta.sigma <- B * sin.sigma * (cos.2sigma.m +
B/4 * (cos.sigma * (-1 + 2 * cos2.2sigma.m) -
B/6 * cos.2sigma.m * (-3 + 4 * sin.sigma * sin.sigma) *
(-3 + 4 * cos2.2sigma.m)))
# alpha1, alpha2: azimuths of the geodesic, clockwise from north
# alpha2 in the direction p1 p2 produced
alpha1 <- atan2(cos.U2 * sin.lambda,
cos.U1 * sin.U2 - sin.U1 * cos.U2 * cos.lambda)
alpha2 <- atan2(cos.U1 * sin.lambda,
-sin.U1 * cos.U2 + cos.U1 * sin.U2 * cos.lambda)
# From Chris Veness (https://www.movable-type.co.uk):
# Special handling of exactly antipodal points where sin2.sigma = 0
# (due to discontinuity atan2(0, 0) = 0 but atan2(NaN, 0) = PI/2 / 90°) - in
# which case bearing is always meridional, due north (or due south!)
alpha1[coincident.points] <- 0
alpha2[coincident.points] <- PI
#s <- round(unname(b * A * (sigma - delta.sigma)), 3) #round to mm
s <- unname(b * A * (sigma - delta.sigma))
if (any(no.convergence)) {
s[no.convergence] <- NA_real_
alpha1[no.convergence] <- NA_real_
alpha2[no.convergence] <- NA_real_
warning("Vicenty formula failed to converge. Try to increase iterations.")
}
list(distance=s, initial.bearing=alpha1, final.bearing=alpha2)
}
# Vicenty (direct) iterative method to compute a destination point from a given
# point and an initial bearing (both in radians). 's' (distance) in m.
#' @noRd
#' @param p1 A matrix of coordinates in radians
#' @param s A numeric vector. Distances (to compute the destination point)
#' @param alpha1 A numeric vector. Initial bearing in radians
#' @param datum A string containing "OSGB36", "WGS84" or "ETRS89"
#' @param iterations Scalar value. Number of iterations to reach convergence
.direct.vicenty.ellipsoid <- function(p1, s, alpha1, datum, iterations = 100L) {
# ellipsoid parameters
ellipsoid <- lonlat.datum[lonlat.datum$datum==datum, "ellipsoid"]
params <- lonlat.ellipsoid[lonlat.ellipsoid$ellipsoid==ellipsoid,
c("a","b","f")]
a2 <- params$a * params$a
b <- params$b
b2 <- b * b
f <- 1 / params$f
one.f <- 1 - f
tan.U1 <- one.f * tan(p1[, 2])
cos.U1 <- 1 / sqrt(1 + tan.U1 * tan.U1)
sin.U1 <- tan.U1 * cos.U1
sin.alpha1 <- sin(alpha1)
cos.alpha1 <- cos(alpha1)
sigma1 <- atan2(tan.U1, cos.alpha1)
sin.alpha <- cos.U1 * sin.alpha1
cos2.alpha <- 1 - sin.alpha * sin.alpha
u2 <- cos2.alpha * (a2 - b2) / b2
A <- 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
B <- u2 / 1024 * (256 + u2 * (-128 + u2 *(74 - 47 * u2)))
sigma <- s / (b * A)
no.convergence <- FALSE
repeat {
cos.2sigma.m <- cos(2 * sigma1 + sigma)
cos2.2sigma.m <- cos.2sigma.m * cos.2sigma.m
sin.sigma <- sin(sigma)
cos.sigma <- cos(sigma)
delta.sigma <- B * sin.sigma * (cos.2sigma.m + B/4 *
(cos.sigma * (-1 + 2 * cos2.2sigma.m) -
B/6 * cos.2sigma.m * (-3 + 4 * sin.sigma * sin.sigma) *
(-3 + 4 * cos2.2sigma.m)))
sigma.tmp <- sigma
sigma <- s / (b * A) + delta.sigma
iterations <- iterations - 1L
if ( max(abs(sigma - sigma.tmp)) <= 1e-12 || iterations == 0L) {
break
}
}
if (iterations == 0L) {
no.convergence <- abs(sigma - sigma.tmp) > 1e-12
}
x <- sin.U1 * sin.sigma - cos.U1 * cos.sigma * cos.alpha1
phi2 <- atan2(sin.U1 * cos.sigma + cos.U1 * sin.sigma * cos.alpha1,
one.f * sqrt(sin.alpha * sin.alpha + x * x))
lambda <- atan2(sin.sigma * sin.alpha1,
cos.U1 * cos.sigma - sin.U1 * sin.sigma * cos.alpha1)
C <- f / 16 * cos2.alpha * (4 + f * (4 - 3 * cos2.alpha))
L <- lambda - (1 - C) * f * sin.alpha *
(sigma + C * sin.sigma * (cos.2sigma.m + C * cos.sigma *
(-1 + 2 * cos2.2sigma.m)))
lambda2 <- p1[, 1] + L
alpha2 <- atan2(sin.alpha, -x)
#p2 <- cbind(x=lambda2, y=phi2)
if (any(no.convergence)) {
#p2[no.convergence, ]<- NA_real_
lambda2[no.convergence]<- NA_real_
phi2[no.convergence]<- NA_real_
alpha2[no.convergence] <- NA_real_
warning("Vicenty formula failed to converge. Try to increase iterations.")
}
# coincident points. Set alpha2 = alpha1
p2.is.p1 <- (s == 0)
if (any(p2.is.p1)) {
#p2[p2.is.p1, ] <- p1[p2.is.p1, ]
lambda2[p2.is.p1] <- p1[p2.is.p1, 1]
phi2[p2.is.p1] <- p1[p2.is.p1, 2]
alpha2[p2.is.p1] <- alpha1[p2.is.p1]
}
list(lon=lambda2, lat=phi2, final.bearing=alpha2)
}
Try the sgo package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.