# R/distHaversine.R In geosphere: Spherical Trigonometry

#### Documented in distHaversine

```# author of original JavaScript code: Chris Vennes
# (c) 2002-2009 Chris Veness
# http://www.movable-type.co.uk/scripts/latlong.html
# Licence: LGPL, without any warranty express or implied

# Port to R by Robert Hijmans
# October 2009
# version 0.1

distHaversine <- function(p1, p2, r=6378137) {
#* Haversine formula to calculate distance between two points specified by
#* from: Haversine formula - R.W. Sinnott, "Virtues of the Haversine",
#*  Sky and Telescope, vol 68, no 2, 1984
#*  http:#//www.census.gov/cgi-bin/geo/gisfaq?Q5.1

# source http://www.movable-type.co.uk/scripts/latlong.html
# (c) 2002-2009 Chris Veness
if (missing(p2)) {
p2 <- p1[-1, ,drop=FALSE]
p1 <- p1[-nrow(p1), ,drop=FALSE]
} else {
}

p = cbind(p1[,1], p1[,2], p2[,1], p2[,2], as.vector(r))

dLat <- p[,4]-p[,2]
dLon <- p[,3]-p[,1]
a <- sin(dLat/2) * sin(dLat/2) + cos(p[,2]) * cos(p[,4]) * sin(dLon/2) * sin(dLon/2)
# to avoid values of 'a' that are a sliver above 1
# which may occur at antipodes
# https://stackoverflow.com/questions/45889616/why-does-disthaversine-return-nan-for-some-pairs-of-coordinates#
a <- pmin(a, 1)
dist <- 2 * atan2(sqrt(a), sqrt(1-a)) * p[,5]
return( as.vector(dist))
}

#	lon1 <- p[,1]
#	lat1 <- p[,2]
#	lon2 <- p[,3]
#	lat2 <- p[,4]
#	r <- p[,5]
#	dLat <- (lat2-lat1)
#	dLon <- (lon2-lon1)
#	a <- sin(dLat/2) * sin(dLat/2) + cos(lat1) * cos(lat2) * sin(dLon/2) * sin(dLon/2)
#	dist <- 2 * atan2(sqrt(a), sqrt(1-a)) * r
```

## Try the geosphere package in your browser

Any scripts or data that you put into this service are public.

geosphere documentation built on May 2, 2019, 5:16 p.m.