R/get_city_dist.R

In Gabegit/gmdata: My some useful functions for data tools.

#' 计算城市间球面距离
#' @author 李国民
#' @param long1: longtitude of city 1
#' @param lat1: latitude of city 1
#' @param long2: longtitude of city 2
#' @param lat2: lattitude of city 2
#' @return a distance numeric
#' @export
#' @examples
#' get_city_dist(115.88333,28.68333,112.97087,28.19874)
#南昌与长沙距离

# Calculates the geodesic distance between two points specified by degrees (DD) latitude/longitude using
# Haversine formula (hf), Spherical Law of Cosines (slc) and Vincenty inverse formula for ellipsoids (vif)
get_city_dist <- function(long1, lat1, long2, lat2) {

	# Convert degrees to radians
	long1 <- deg2rad(long1)
	lat1 <- deg2rad(lat1)
	long2 <- deg2rad(long2)
	lat2 <- deg2rad(lat2)

	return(list(haversine = gcd.hf(long1, lat1, long2, lat2),
							sphere = gcd.slc(long1, lat1, long2, lat2)
						) )
}

# 计算城市间距离
# Calculates the geodesic distance between two points specified by radian latitude/longitude using the
# Spherical Law of Cosines (slc)
gcd.slc <- function(long1, lat1, long2, lat2) {
	R <- 6371 # Earth mean radius [km]
	d <- acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(long2-long1)) * R
	return(d) # Distance in km
}

# Calculates the geodesic distance between two points specified by radian latitude/longitude using the
# Haversine formula (hf)
gcd.hf <- function(long1, lat1, long2, lat2) {
	R <- 6371 # Earth mean radius [km]
	delta.long <- (long2 - long1)
	delta.lat <- (lat2 - lat1)
	a <- sin(delta.lat/2)^2 + cos(lat1) * cos(lat2) * sin(delta.long/2)^2
	c <- 2 * asin(min(1,sqrt(a)))
	d = R * c
	return(d) # Distance in km
}


# Convert degrees to radians
deg2rad <- function(deg) return(deg*pi/180)

# Calculates the geodesic distance between two points specified by radian latitude/longitude using
# Vincenty inverse formula for ellipsoids (vif)
gcd.vif <- function(long1, lat1, long2, lat2) {

	# WGS-84 ellipsoid parameters
	a <- 6378137         # length of major axis of the ellipsoid (radius at equator)
	b <- 6356752.314245  # ength of minor axis of the ellipsoid (radius at the poles)
	f <- 1/298.257223563 # flattening of the ellipsoid

	L <- long2-long1 # difference in longitude
	U1 <- atan((1-f) * tan(lat1)) # reduced latitude
	U2 <- atan((1-f) * tan(lat2)) # reduced latitude
	sinU1 <- sin(U1)
	cosU1 <- cos(U1)
	sinU2 <- sin(U2)
	cosU2 <- cos(U2)

	cosSqAlpha <- NULL
	sinSigma <- NULL
	cosSigma <- NULL
	cos2SigmaM <- NULL
	sigma <- NULL

	lambda <- L
	lambdaP <- 0
	iterLimit <- 100
	while (abs(lambda-lambdaP) > 1e-12 & iterLimit>0) {
		sinLambda <- sin(lambda)
		cosLambda <- cos(lambda)
		sinSigma <- sqrt( (cosU2*sinLambda) * (cosU2*sinLambda) +
												(cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda) )
		if (sinSigma==0) return(0)  # Co-incident points
		cosSigma <- sinU1*sinU2 + cosU1*cosU2*cosLambda
		sigma <- atan2(sinSigma, cosSigma)
		sinAlpha <- cosU1 * cosU2 * sinLambda / sinSigma
		cosSqAlpha <- 1 - sinAlpha*sinAlpha
		cos2SigmaM <- cosSigma - 2*sinU1*sinU2/cosSqAlpha
		if (is.na(cos2SigmaM)) cos2SigmaM <- 0  # Equatorial line: cosSqAlpha=0
		C <- f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))
		lambdaP <- lambda
		lambda <- L + (1-C) * f * sinAlpha *
			(sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))
		iterLimit <- iterLimit - 1
	}
	if (iterLimit==0) return(NA)  # formula failed to converge
	uSq <- cosSqAlpha * (a*a - b*b) / (b*b)
	A <- 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
	B <- uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)))
	deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM^2) -
																					 	B/6*cos2SigmaM*(-3+4*sinSigma^2)*(-3+4*cos2SigmaM^2)))
	s <- b*A*(sigma-deltaSigma) / 1000

	return(s) # Distance in km
}