R/geodesics.R
In guenardg/codep: Multiscale Codependence Analysis
Defines functions
geodesics
Documented in
geodesics
## **************************************************************************
##
## (c) 2018-2023 Guillaume Guénard
## Department de sciences biologiques,
## Université de Montréal
## Montreal, QC, Canada
##
## **Geodesic Distances Calculation**
##
## This file is part of codep
##
## codep is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## codep is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with codep. If not, see <https://www.gnu.org/licenses/>.
##
## R source code file
##
## **************************************************************************
##
#' Calculation of Geodesic Distances
#'
#' Function \code{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.
#'
#' @param x A set of geographic coordinates in the form of a two-column
#' \code{\link{matrix}} or \code{\link{data.frame}}.
#' @param y An other two-column \code{\link{matrix}} or \code{\link{data.frame}}
#' containing an optional second set of coordinates.
#' @param method The calculation method used to obtain the distances (default:
#' haversine method; see details).
#' @param radius Radius of the planetary body (when assuming a sphere; default:
#' 6371000 m).
#' @param sma Length of the semi-major axis of the planetary body (when assuming
#' a revolution ellipsoid; default: 6378137 m).
#' @param flat Flattening of the ellipsoid (default: 1/298.257223563).
#' @param maxiter Maximum number of iterations, whenever iterative calculation
#' is involved (default: 1024).
#' @param tol Tolerance used when iterative calculation is involved (default:
#' \code{.Machine$double.eps^0.75}; a machine dependent value).
#'
#' @return A `\link{dist}-class` object or, whenever argument \code{y} is
#' provided, a \code{\link{matrix}} with as many rows as the number of rows in
#' argument \code{x} and as many columns as the number of rows in argument
#' \code{y}.
#'
#' @details When only one set of coordinates is given to the function (i.e.,
#' when argument \code{y} is omitted), the function returns the pairwise
#' distances in the form of a `\link{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 \code{x} and
#' each coordinate of argument \code{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 \code{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 \code{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 \code{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
#' \code{maxiter} and \code{tol}, are ignored when using the haversine formula,
#' whereas the value of argument \code{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 \code{niter} attribute that gives the number of iterations that
#' were necessary to achieve convergence. Numbers greater than argument
#' \code{maxiter} are indicative of failed convergence; a warning is issued in
#' such a circumstance.
#'
#' Geodesic distance matrices are non metric.
#'
#' @author \packageAuthor{codep}
#' Maintainer: \packageMaintainer{codep}
#'
#' @seealso The \code{\link{dist}-class} and associated methods.
#'
#' @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
#'
#' @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
#' ##
#'
#' @useDynLib codep, .registration = TRUE
#'
#' @export
geodesics <- function(x, y, method = c("haversine", "Vincenty"),
radius = 6.371e6, sma = 6378137.0, flat = 1/298.257223563,
maxiter = 1024L, tol = .Machine$double.eps^0.75) {
if (NCOL(x)!=2L)
stop("'x' must be Lat Lon geodesic coordinates!")
if (!is.matrix(x))
x <- as.matrix(x)
storage.mode(x) <- "double"
method <- match.arg(method)
if (!missing(y)) {
if(NCOL(y)!=2L)
stop("'y' must be Lat Lon geodesic coordinates!")
if(!is.matrix(y))
y <- as.matrix(y)
storage.mode(y) <- "double"
N <- c(NROW(x),NROW(y))
out <- matrix(0, N[2L], N[1L])
rownames(out) <- dimnames(y)[[1L]]
colnames(out) <- dimnames(x)[[1L]]
} else {
N <- c(NROW(x),NROW(x))
out <- double(N[1L]*(N[1L]-1L)/2)
}
if (method == "haversine")
res <- .C(
"dist_geo_hvs",
x,
if(missing(y)) x else y,
N,
missing(y),
out,
radius
)
if (method == "Vincenty") {
niter <- integer(length(out))
res <- .C(
"dist_geo_vif",
x,
if(missing(y)) x else y,
N,
missing(y),
out,
niter,
sma,
flat,
maxiter,
tol
)
}
out[] <- res[[5L]]
attr(out, "method") <- sprintf("geodesic:%s", method)
attr(out, "call") <- match.call()
if (missing(y)) {
attr(out, "Size") <- N[1L]
attr(out, "Labels") <- dimnames(x)[[1L]]
attr(out, "Diag") <- FALSE
attr(out, "Upper") <- FALSE
if (method == "Vincenty")
attr(out, "niter") <- res[[6L]]
attr(out, "class") <- "dist"
}
if (method == "Vincenty" && any(res[[6L]] > maxiter))
warning(
"Convergence failed for ",
sum(res[[6L]]>maxiter),
" distance(s)!"
)
return(out)
}
#'
guenardg/codep documentation built on April 16, 2024, 9:01 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 geodesics
Documented in geodesics
## **************************************************************************
##
## (c) 2018-2023 Guillaume Guénard
## Department de sciences biologiques,
## Université de Montréal
## Montreal, QC, Canada
##
## **Geodesic Distances Calculation**
##
## This file is part of codep
##
## codep is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## codep is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with codep. If not, see <https://www.gnu.org/licenses/>.
##
## R source code file
##
## **************************************************************************
##
#' Calculation of Geodesic Distances
#'
#' Function \code{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.
#'
#' @param x A set of geographic coordinates in the form of a two-column
#' \code{\link{matrix}} or \code{\link{data.frame}}.
#' @param y An other two-column \code{\link{matrix}} or \code{\link{data.frame}}
#' containing an optional second set of coordinates.
#' @param method The calculation method used to obtain the distances (default:
#' haversine method; see details).
#' @param radius Radius of the planetary body (when assuming a sphere; default:
#' 6371000 m).
#' @param sma Length of the semi-major axis of the planetary body (when assuming
#' a revolution ellipsoid; default: 6378137 m).
#' @param flat Flattening of the ellipsoid (default: 1/298.257223563).
#' @param maxiter Maximum number of iterations, whenever iterative calculation
#' is involved (default: 1024).
#' @param tol Tolerance used when iterative calculation is involved (default:
#' \code{.Machine$double.eps^0.75}; a machine dependent value).
#'
#' @return A `\link{dist}-class` object or, whenever argument \code{y} is
#' provided, a \code{\link{matrix}} with as many rows as the number of rows in
#' argument \code{x} and as many columns as the number of rows in argument
#' \code{y}.
#'
#' @details When only one set of coordinates is given to the function (i.e.,
#' when argument \code{y} is omitted), the function returns the pairwise
#' distances in the form of a `\link{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 \code{x} and
#' each coordinate of argument \code{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 \code{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 \code{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 \code{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
#' \code{maxiter} and \code{tol}, are ignored when using the haversine formula,
#' whereas the value of argument \code{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 \code{niter} attribute that gives the number of iterations that
#' were necessary to achieve convergence. Numbers greater than argument
#' \code{maxiter} are indicative of failed convergence; a warning is issued in
#' such a circumstance.
#'
#' Geodesic distance matrices are non metric.
#'
#' @author \packageAuthor{codep}
#' Maintainer: \packageMaintainer{codep}
#'
#' @seealso The \code{\link{dist}-class} and associated methods.
#'
#' @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
#'
#' @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
#' ##
#'
#' @useDynLib codep, .registration = TRUE
#'
#' @export
geodesics <- function(x, y, method = c("haversine", "Vincenty"),
radius = 6.371e6, sma = 6378137.0, flat = 1/298.257223563,
maxiter = 1024L, tol = .Machine$double.eps^0.75) {
if (NCOL(x)!=2L)
stop("'x' must be Lat Lon geodesic coordinates!")
if (!is.matrix(x))
x <- as.matrix(x)
storage.mode(x) <- "double"
method <- match.arg(method)
if (!missing(y)) {
if(NCOL(y)!=2L)
stop("'y' must be Lat Lon geodesic coordinates!")
if(!is.matrix(y))
y <- as.matrix(y)
storage.mode(y) <- "double"
N <- c(NROW(x),NROW(y))
out <- matrix(0, N[2L], N[1L])
rownames(out) <- dimnames(y)[[1L]]
colnames(out) <- dimnames(x)[[1L]]
} else {
N <- c(NROW(x),NROW(x))
out <- double(N[1L]*(N[1L]-1L)/2)
}
if (method == "haversine")
res <- .C(
"dist_geo_hvs",
x,
if(missing(y)) x else y,
N,
missing(y),
out,
radius
)
if (method == "Vincenty") {
niter <- integer(length(out))
res <- .C(
"dist_geo_vif",
x,
if(missing(y)) x else y,
N,
missing(y),
out,
niter,
sma,
flat,
maxiter,
tol
)
}
out[] <- res[[5L]]
attr(out, "method") <- sprintf("geodesic:%s", method)
attr(out, "call") <- match.call()
if (missing(y)) {
attr(out, "Size") <- N[1L]
attr(out, "Labels") <- dimnames(x)[[1L]]
attr(out, "Diag") <- FALSE
attr(out, "Upper") <- FALSE
if (method == "Vincenty")
attr(out, "niter") <- res[[6L]]
attr(out, "class") <- "dist"
}
if (method == "Vincenty" && any(res[[6L]] > maxiter))
warning(
"Convergence failed for ",
sum(res[[6L]]>maxiter),
" distance(s)!"
)
return(out)
}
#'
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.