Python Haversine – Compute Distance Between Coordinates and Print CSV

distancegpslatitude longitudepython

I have a .csv file with 3 columns. The first column is titled Distance and it is empty. The second and third columns are titled Longitude and Latitude which contain coordinates, for example:

Distance,Longitude,Latitude

,-77.60483,40.31117

,-77.58167,40.3045

,-77.52883,40.24917

,-77.508,40.14917

,-77.49617,40.13117

I have wrote a script that computes the distance between the two points using the haversine formula (see below), but I have to manually enter the coordinates, which is not desirable:

import math

lat1 = 40.31117
long1 = -77.60483
lat2 = 40.3045
long2 = -77.58167

def distance_on_unit_sphere(lat1, long1, lat2, long2):

# Converts lat & long to spherical coordinates in radians.
degrees_to_radians = math.pi/180.0

# phi = 90 - latitude
phi1 = (90.0 - lat1)*degrees_to_radians
phi2 = (90.0 - lat2)*degrees_to_radians

# theta = longitude
theta1 = long1*degrees_to_radians
theta2 = long2*degrees_to_radians

# Compute the spherical distance from spherical coordinates.
# For two locations in spherical coordinates:
# (1, theta, phi) and (1, theta', phi')cosine( arc length ) =
# sin phi sin phi' cos(theta-theta') + cos phi cos phi' distance = rho * arc    length

cos = (math.sin(phi1)*math.sin(phi2)*math.cos(theta1 - theta2) +
   math.cos(phi1)*math.cos(phi2))
arc = math.acos(cos)*6371 #radius of the earth in km

return arc

print distance_on_unit_sphere(lat1,long1,lat2,long2)

My goal is to populate the empty distance column using Python by computing the distance between the 1st location (-77.60483,40.31117) and the 2nd location (-77.58167,40.3045) and then print that value into the Distance column and then compute the distance between the 2nd (-77.58167,40.3045) and 3rd point (-77.52883,40.24917) and then print that value below the previously computed value in Distance column and so on for the subsequent coordinates. I need some help or an example provided to input the .csv file, run the haversine formula, and then write the result into the appropriate location in an output .csv.

Best Answer

If I've understood correctly, this code should do what you're after. Using your input

Distance,Longitude,Latitude
,-77.60483,40.31117
,-77.58167,40.3045
,-77.52883,40.24917
,-77.508,40.14917
,-77.49617,40.13117

I'm producing this output. I've assumed that you want the distances to go forward - i.e. that line 2 has its distance set to the distance to the next point (line 3).The last line is set to 0 as there is no following line to measure to

Distance,Longitude,Latitude
2.0992,-77.6048,40.3112
7.6122,-77.5817,40.3045
11.2593,-77.5288,40.2492
2.2399,-77.5080,40.1492
0.0000,-77.4962,40.1312

Keep your existing code, and append this snippet. The csv library makes it nice and easy to work with CSV files.

I'm also using a deque. This is a FIFO (first in, first out) data structure. I use this to keep a queue of points. Once this has two entries, I compute the distance between the two points, write it out, then discard the oldest point.

import csv
import collections
filein = open("test.csv","rU") # universal newline mode
fileout = open("output.csv","w")
fileout.write("Distance,Longitude,Latitude\n")
csvrdr = csv.reader(filein, delimiter=",")
linecount = 0
# use a deque. once this hits 2 entries, compute the distance
# and then remove the oldest entry
todo = collections.deque()
for line in csvrdr:
    if linecount > 0: # skip header
        _, longitude, latitude = line # _ means discard
        longitude = float(longitude)
        latitude = float(latitude)
        todo.append((longitude, latitude))
    if len(todo)==2:
        # got two points, compute distance and remove the oldest
        oldlon, oldlat = todo.popleft()
        distance = distance_on_unit_sphere(latitude, longitude, oldlat, oldlon)
        fileout.write("%2.4f,%2.4f,%2.4f\n" % (distance, oldlon, oldlat))
    linecount += 1
# got one entry in deque, spit that out with 0
assert(len(todo)==1)
oldlon, oldlat = todo.popleft()
fileout.write("%2.4f,%2.4f,%2.4f\n" % (0.0, oldlon, oldlat))
fileout.close()
filein.close()

Related Solutions

[GIS] Finding the minimum and maximum distance between an observer and a lat/lon bounding box

The "bounding box" can have several interpretations. The two most natural are (a) the sides are meridians and the top and bottom are circles of constant latitude and (b) the sides are still meridians, but now the top and bottom are geodesics (which therefore do not necessarily have constant latitude). Fortunately, the solution method is the same in both interpretations, because in all cases the sides are portions of circles. This reduces the problem to optimizing the distance between a fixed point x and an arc of a circle (in three dimensions). Do this for each of the four sides and retain the best results.

A simple way to represent the arc of a circle is to choose coordinates for x in which the circle's center is at the origin, the circle is determined by two vectors e and f of equal length (say, r) and at right angles to each other, and the arc is given by all points of the form cos(q) * e + sin(q) * f for angles q ranging from q0 through q1. Critical points for the optima occur either at the endpoints of the arc or where the tangent vector is perpendicular to the direction vector. The latter tells us the inner product of x - (cos(q) * e + sin(q) * f) with (-sin(q) * e + cos(q) * f) is zero. This easily reduces to the equation

tan(q) = -xe / xf

where xe is the inner product of x with e and xf is the inner product of x with f. There will usually be two solutions (differing by 180 degrees). Use ATan (or its equivalent) to find them. Any of these solutions that lies in the interval [q0, q1] is a critical point.

In this fashion you collect (in most cases) up to two critical points per side plus the four vertices of the bounding box: that gives between 4 and 12 distance calculations to make in order to find the largest and smallest distance. These are 3D (Pythagorean) distances: compute them using the Pythagorean formula, but don't bother to take square roots when comparing; just take square roots for the largest and smallest distances.

^{Portions of the four circles determined by the sides of this bounding box are shown in gray. The arcs of latitude are in red and the arcs of the meridians are in blue.}

There is a notable exception to this nice situation: x might lie on, or extremely close to, the axis of the circle. In this case, every point along the circle (and therefore every point along the arc) has the same distance to x. Make sure to check for this situation. It's easy to detect: both xf and xe are zero (or extremely small compared to x itself).

It remains, then, to describe how to find e, f, q0, and q1. That's easy:

For a meridian at longitude lambda, the circle's center is already at (0,0,0), so don't change the coordinates of x (which presumably is referenced to a geocentric system). Take e = (Rcos(lambda), Rsin(lambda), 0) and f = (0,0,R) where R is the earth's radius (in a spherical model). Here, q0 will be the minimum latitude among the bounding box coordinates and q1 will be the maximum latitude.
For a circle of latitude phi, take the center to be (0, 0, Rsin(phi)). Subtract these coordinates from x. In these new coordinates, we can let e = (Rcos(phi), 0, 0) and f = (0, R*cos(phi), 0). Now q0 is the westernmost longitude and q1 is the easternmost longitude in the bounding box.
Equations for e, f, q0, and q1 are more complicated for arbitrary geodesics (case (b) in the first paragraph). Because it is not apparent these are needed, I won't derive them here.

In summary,

Only simple calculations are needed to establish the 3D coordinate representation of x, e, f and to find q0 and q1. (This is all done four times, once for each side of the bounding box.)
Four arctangent calculations will determine between 0 and 8 critical points. At this time, you can inspect xe and xf (four times) to check whether any of the sides of the bounding box are equidistant from x. (You will have to decide how such a situation is to be handled).
Squared pythagorean distances (sums of squares) have to be determined between x and the four corners of the bounding box plus these other critical points. Take the square roots of the min and max of these distances.

GPS Satellite Coverage – Finding Ring of Coverage on WGS-84 Ellipsoid

The solution for an ellipsoid is pretty messy--it is an irregular shape, not a circle--and is best computed numerically rather than with a formula.

On a world map the difference between the WGS84 solution and a purely spherical solution will only barely be noticeable (it's about one pixel on a screen). The same difference would be created by changing the mask angle by about 0.2 degrees or in using a polygonal approximation. If these errors are acceptably small, then you can exploit the symmetry of the sphere to obtain a simple formula.

This map (using an Equirectangular projection) shows the coverage for a satellite at 22,164 kilometers (from the earth's center) with a mask angle of m = 15 degrees on the WGS84 spheroid. Recomputing the coverage for a sphere does not visibly change this map.

On the sphere, the coverage will truly be a circle centered at the satellite's location, so we only need to figure out its radius, which is an angle. Call this t. In cross-section there is a triangle OSP formed by the earth's center (O), the satellite (S), and any point (P) on the circle:

The side OP is the radius of the earth, R.
The side OS is the height of the satellite (above the earth's center). Call this h.
The angle OPS is 90 + m.
The angle SOP is t, which we want to find.
Because the three angles of a triangle sum to 180 degrees, the third angle OSP must equal 90 - (m + t).

The solution is now a matter of elementary trigonometry. The (planar) law of sines asserts that

sin(90 - (m+t)) / r = sin(90 + m) / h.

The solution is

t = ArcCos(cos(m) / (h/r)) - m.

As a check, consider some extreme cases:

When m = 0, t = ArcCos(r/h), which can be verified with elementary Euclidean geometry.
When h = r (the satellite hasn't launched), t = ArcCos(cos(m) / 1) - m = m - m = 0.
When m = 90 degrees, t = ArcCos(0) - 90 = 90 - 90 = 0, as it should be.

This reduces the problem to drawing a circle on the sphere, which can be solved in many ways. For instance, you can buffer the satellite's location by t * R * pi/180 using an equidistant projection centered at the satellite. Techniques for working with circles on the sphere directly are illustrated at https://gis.stackexchange.com/a/53323/664.

Edit

FWIW, for GPS satellites and small mask angles (less than 20 degrees or so), this non-trigonometric approximation is accurate (to a few tenths of a degree and less than a few hundredths of a degree when the mask angle is under 10 degrees):

t (degrees) = -0.0000152198628163333 * (-5.93410042925107*10^6 + 
              3.88800000000000*10^6 r/h + 65703.6145507725 m + 
              9.86960440108936 m^2 - 631.654681669719 r/h m^2)

For example, with a mask angle of m = 10 degrees and a satellite at 26,559.7 km above the earth's center (which is the nominal distance of a GPS satellite), this approximation gives 66.32159..., whereas the value (correct for the sphere) is 66.32023... .

(The approximation is based on a Taylor series expansion around m = 0, r/h = 1/4.)

Best Answer

Related Solutions

[GIS] Finding the minimum and maximum distance between an observer and a lat/lon bounding box

GPS Satellite Coverage – Finding Ring of Coverage on WGS-84 Ellipsoid

Edit

Related Question