R Programming
Just for fun, I am placing a short article on Vincenty’s Formula for calculating the ending latitude and longitude from a given lat/lon and a departure angle and distance (in meters) here using the “R” programming language.
“R” is mostly used for “big data” or statistics, but it works remarkably well in the scientific field as well, as this short program shows.
I will only present the R program itself without presenting the “pretty print” formulae. Those can be found in many places on the Internet, such as Wikipedia, among others. Also, as R already supplies a library for computing the “inverse” solutions, I saw no need to attempt to “reinvent the wheel,” as the saying goes.
Now, the “inverse functions” are contained in the “Geosphere” package, available at The Comprehensive R Archive Network, and are called “distVincentyEllipsoid” and “finalBearing.” They are more accurate than the “Haversine” formulae, which use a spheroid. They default to using the WGS84 ellipsoid, which more closely describes the actual shape of the Earth.
Anyway, on to the subject… I had a need to calculate distances between two sets of GPS datapoints, and determining those points if not known. This led me to the above referenced formulae. These calculations can be useful for Geocaching as well as surveying applications.
A little history. I started this project by attempting to determine more exact GPS coordinates than the commonly available C/A codes, as the P codes are not available except to the military. Even though the Selective Availability (SA) function was turned off (?) many years ago, the GPS packages still show a disturbingly large variance in readings. So, to do this I am averaging GPS readings at different times of the day, with various numbers of available satellites. I am creeping up on 25,000 readings from two disparately located GPS receivers, hence the need to determine the location differences between the two, to sub-millimeter accuracy.
I decided to use R as a platform simply because I have been investigating it, and to explore the mathematical capabilities inherent to it. I have used Python as well for the same calculations, but wanted to do it in R.
So without further ado…
We start with some values for the radius of the equator, flattening of the ellipsoid for WGS-84, and radius at the poles.
# WGS-84 constants.
a <- 6378137
b <- 6356752.314245
f <- 0.00335281066
print("Vincenty direct calculation.")
phi <- as.double(readline(prompt="Enter latitude: "))
L1 <- as.double(readline(prompt="Enter longitude: "))
a1 <- as.double(readline(prompt="Enter azimuth: "))
s <- as.double(readline(prompt="Enter distance: ")))
# Test values:
phi = 32; L1 = -106.5; a1 = -23; s = 60243
Answers: lat2 = 32.49982392, lon2 = -106.13365295
Depending on how your platform performs math to what significant number of digits, you may see different numbers, but they should not be too far off. In DMS notation, within fractions of a second.
The remainder of the program is the actual formulas, and consist of some initial calculations, followed by iterations of the core formulas until the result variation is insignificant, within whatever precision you wish. Repeating what I have read, 1e-12 gives precision in the sub-millimeter range. I have not verified that!
This section is to establish a initial value going into the loop.
# First pass to get an initial value of sigma
sigma_m <- 2 * sigma1 + sigma
delta_sigma <- B * sin(sigma) * (cos(sigma_m) + .25 * B * (cos(sigma) *
(-1 + 2 * cos(sigma_m)^2) - (B/6) * cos(sigma_m) * (-3 + 4 * sin(sigma)^2
* (-3 + 4 * cos(sigma_m)^2))))
sigma_prev <- sigma
sigma <- (s/(b * A)) + delta_sigma
I have not discovered a way in R to run through a loop before evaluation of the loop’s statements. Hence, the above section. Fortunately, on the Casio Prizm series of calculators, there is a “Do/LpWhile” construct which does that, so the above repeated formulas are not necessary.
# Loop until no significant change in sigma.
while ( abs(sigma_prev - sigma) > 1e-12) {
sigma_m <- 2 * sigma1 + sigma
delta_sigma <- B * sin(sigma) * (cos(sigma_m) + .25 * B *
(cos(sigma) * (-1 + 2 * cos(sigma_m)^2) - (B/6) *
cos(sigma_m) * (-3 + 4 * sin(sigma)^2 *
(-3 + 4 * cos(sigma_m)^2))))
sigma_prev <- sigma
sigma <- (s/(b * A)) + delta_sigma
}
phi2 <- atan2(sin(U1) * cos(sigma) + cos(U1) * sin(sigma) * cos(a1), (1-f) *
sqrt(sin(alpha)^2 + (sin(U1) * sin(sigma) - cos(U1) * cos(sigma) *
cos(a1))^2))
lambda <- atan2(sin(sigma) * sin(a1), cos(U1) * cos(sigma) - sin(U1) *
sin(sigma) * cos(a1))
C <- (f/16) * cos(alpha)^2 * (4 + f * (4 - 3 * cos(alpha)^2))
L <- lambda - (1-C) * f * sin(alpha) * ( sigma + C * sin(sigma) *
(cos(sigma_m) + C * cos(sigma) *
( -1 + 2 * cos(sigma_m)^2)))
L2 <- L + L1
a2 <- atan2(sin(alpha), - sin(U1) * sin(sigma) + cos(U1) * cos(sigma) *
cos(a1))
To make things more readable, we add the following section at the end. As we are doing this a couple of times, we use a function (DMS) and use it to create and return the output.
# Convert input value to deg, min, sec.
DMS <- function(az) {
degrees <- as.integer(az)
minutes <- as.integer(abs(az - degrees)*60)
seconds <- ((abs(az - degrees)*60)-minutes)*60
seconds1 <- format(seconds, digits=4, nsmall=2)
# Unicode \U00b0 is degree symbol
dms <- paste(" DMS =",degrees,'\U00b0',minutes,"'",seconds1,"'' ")
return(dms)
}
# Convert values back to degrees.
a2D <- rad2deg(a2)
phi2D <- rad2deg(phi2)
L2D <- rad2deg(L2)
# Show outputs.
print("")
print(paste("Lat2:", phi2D, "Lon2:", L2D))
print(paste("Azimuth:", a2D, "Distance:", s, "m"))
print(paste("Starting angle:", DMS(a1D)))
print(paste("Ending angle:", DMS(a2D)))
The entire program is available here as an “R” file. To use it within R type: source(“vincentyDirect.R”) or whatever you decide to name it. Also, as I am casually trying to learn more about R, I’m sure this program could be made much more compact and efficient, but I did it so I could make it work, not be good! Perhaps, it will be useful to someone, and I’m sure they will improve it.
God Bless and have a great day! Bye for now…