# R/pointdistance.R In raster: Geographic Data Analysis and Modeling

#### Documented in pointDistance

```# Author: Robert J. Hijmans  and Jacob van Etten
# Date :  June 2008
# Version 0.9
# Licence GPL v3

.pointsToMatrix <- function(p) {
if (inherits(p, 'SpatialPoints')) {
p <- coordinates(p)
} else if (is.data.frame(p)) {
p <- as.matrix(p)
} else if (is.vector(p)){
if (length(p) != 2) {
stop('Wrong length for a vector, should be 2')
} else {
p <- matrix(p, ncol=2)
}
}
if (is.matrix(p)) {
if (ncol(p) != 2) {
stop( 'A points matrix should have 2 columns')
}
cn <- colnames(p)
if (length(cn) == 2) {
if (toupper(cn[1]) == 'Y' | toupper(cn[2]) == 'X')  {
stop('Highly suspect column names (x and y reversed?)')
}
if (toupper(substr(cn[1],1,3) == 'LAT' | toupper(substr(cn[2],1,3)) == 'LON'))  {
stop('Highly suspect column names (longitude and latitude reversed?)')
}
}
} else {
stop('points should be vectors of length 2, matrices with 2 columns, or a SpatialPoints* object')
}

return(p)
}

.distm <- function (x, longlat) {
if (longlat) {
n <- nrow(x)
dm <- matrix(ncol = n, nrow = n)
dm[cbind(1:n, 1:n)] <- 0
if (n > 1) {
for (i in 2:n) {
j = 1:(i - 1)
dm[i, j] = .geodist(x[i, 1], x[i, 2], x[j, 1], x[j, 2])
}
}
return(dm)
} else {
return(.planedist2(x, x))
}
}

.distm2 <- function (x, y, longlat) {
if (longlat) {
n <- nrow(x)
m <- nrow(y)
dm <- matrix(ncol=m, nrow=n)
for (i in 1:n) {
dm[i,] <- .geodist(x[i, 1], x[i, 2], y[, 1], y[, 2])
}
return(dm)
} else {
return(.planedist2(x, y))
# fun <- .planedist
}
}

.distm2new <- function (x, y, longlat, a=6378137, f=1/298.257223563) {
if (longlat) {
n <- nrow(x)
m <- nrow(y)

xx <- cbind(rep(x[,1], m), rep(x[,2], m))
yy <- cbind(rep(y[,1], each=n), rep(y[,2], each=n))

g <- .Call("_raster_point_distance", xx, yy, TRUE, a, f, PACKAGE='raster')
return(matrix(g, n, m))
} else {
return(.planedist2(x, y))
# fun <- .planedist
}
}

pointDistance <- function (p1, p2, lonlat, allpairs=FALSE, ...) {

longlat <- list(...)\$longlat
if (!is.null(longlat)) {
lonlat <- longlat
}

if (missing(lonlat)) {
if (isLonLat(p1)) {
lonlat <- TRUE
} else if (! is.na(projection(p1)) ) {
lonlat <- FALSE
} else {
stop('you must provide a "lonlat" argument (TRUE/FALSE)')
}
}
stopifnot(is.logical(lonlat))

p1 <- .pointsToMatrix(p1)
if (missing(p2)) {
return(.distm(p1, lonlat))
}

p2 <- .pointsToMatrix(p2)

if (nrow(p1) != nrow(p2)) {
allpairs <- TRUE
}

if (allpairs) {
if(nrow(p1) > 1 & nrow(p2) > 1) {
return(.distm2(p1, p2, lonlat))
}
}

if (lonlat ) {
#		return( .haversine(p1[,1], p1[,2], p2[,1], p2[,2], r=6378137) )
return( .geodist(p1[,1], p1[,2], p2[,1], p2[,2]) )
} else {
return( .planedist(p1[,1], p1[,2], p2[,1], p2[,2]) )
}
}

.planedist <- function(x1, y1, x2, y2) {
sqrt(( x1 -  x2)^2 + (y1 - y2)^2)
}

.planedist2 <- function(p1, p2) {
# code by Bill Venables
# https://stat.ethz.ch/pipermail/r-help/2008-February/153841.html
z0 <- complex(, p1[,1], p1[,2])
z1 <- complex(, p2[,1], p2[,2])
outer(z0, z1, function(z0, z1) Mod(z0-z1))
}

.geodist <- function(x1, y1, x2, y2, a=6378137, f=1/298.257223563) {
# recycle
p <- cbind(x1, y1, x2, y2)
.Call("_raster_point_distance", p[,1:2, drop=FALSE], p[, 3:4,drop=FALSE], TRUE, a, f, PACKAGE='raster')
#	.Call("inversegeodesic", as.double(p[,1]), as.double(p[,2]), as.double(p[,3]), as.double(p[,4]), as.double(a), as.double(f), PACKAGE='raster')
}

.old_haversine <- function(x1, y1, x2, y2, r=6378137) {