I return to this issue because it is very similar to How do I find vector line bearing in QGIS or GRASS? and it can be solved with Python in the same way:
1) Haversine distance
One can find lots of scripts by searching Haversine distance with Python on the Internet and I choose one of them in Haversine Formula in Python (Bearing and Distance between two GPS points)
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
c = 2 * math.asin(math.sqrt(a))
km = 6367 * c
return km
We have a series of lines (points) in the file that must be treated in pairs (point1 - point2) to calculate the distance. For this we will use a simple iterator from Most pythonic way to get the previous element
def offset(iterable):
prev = None
for elem in iterable:
yield prev, elem
prev = elem
Now it is possible to read the file (example of Kerrie) in pairs of lines/points
import csv
with open('testhavers.csv', 'rb') as f:
reader = csv.DictReader(f)
for pair in offset(reader):
print pair
(None, {'LAT': '10.08527', 'LON': '124.59833', 'ID': '1', 'TIME': '21:24:37'})
({'LAT': '10.08527', 'LON': '124.59833', 'ID': '1', 'TIME': '21:24:37'},
{'LAT': '10.08523', 'LON': '124.59830', 'ID': '2', 'TIME': '21:25:07'})
({'LAT': '10.08523', 'LON': '124.59830', 'ID': '2', 'TIME': '21:25:07'},
{'LAT': '10.08526', 'LON': '124.59832', 'ID': '3', 'TIME': '21:25:37'})
({'LAT': '10.08526', 'LON': '124.59832', 'ID': '3', 'TIME': '21:25:37'},
{'LAT': '10.08526', 'LON': '124.59831', 'ID': '4', 'TIME': '21:26:07'})
Then create a shapefile containing the original fields of the csv file and a new field for the distance with the Python modules Shapely and Fiona of Sean Gillies:
import fiona
from shapely.geometry import Point, mapping
# creation of the schema of the shapefile (geometry and fields)
schema = { 'geometry': 'Point', 'properties':{'ID': 'int', 'LAT':'float', 'LON':'float', 'TIME':'str','distance' : 'float'}}
# creation of the shapefile:
with fiona.collection("result.shp", "w", "ESRI Shapefile", schema) as output:
# reading the csv file
with open('testhavers.csv', 'rb') as f:
reader = csv.DictReader(f)
# we need here to eliminate the first pair of point with None
for i, pair in enumerate(offset(reader)):
if i == 0: (pair with None)
# writing of the point geometry and the attributes
point = Point(float(pair[1]['LON']), float(pair[1]['LAT']))
dist = 0 # None
output.write({'properties': {'ID':int(pair[1]['ID']),'LAT':float(pair[1]['LAT']),'LON':float(pair[1]['LON']), 'TIME':pair[1]['TIME'],'distance': dist},'geometry': mapping(point)})
else:
# writing of the point geometry and the attributes
point = Point(float(pair[1]['LON']), float(pair[1]['LAT']))
# Haversine distance between pairs of points
dist = haversine(float(pair[0]['LON']), float(pair[0]['LAT']), float(pair[1]['LON']),float(pair[1]['LAT']))
output.write({'properties': {'ID':int(pair[1]['ID']),'LAT':float(pair[1]['LAT']),'LON':float(pair[1]['LON']), 'TIME':pair[1]['TIME'],'distance': dist},'geometry': mapping(point)})
and the result:
It is also possible to do it with PyQGIS but it is more complex than Fiona which uses simple dictionaries to create shapefiles.
You can use another function to calculate the Haversine distance (Why is law of cosines more preferable than haversine when calculating distance between two latitude-longitude points?) without any problem, only the distance calculation changes, not the process of creating the shapefile.
Best Answer
Do a spatial join!
First, set up your data frame in a projected coordinate system of your choice (whatever units you want your distances to show up in). So, say you're working in State Plane Feet, make sure all your layers are in State Plane Feet, so if they're not project them into it. From there, Right click on the points layer and click Joins & Relates > Joins.... Make sure to select "Join Layer from another layer based on spatial location." In there, select your layers that you're spatially joining together in the drop down list.
Then, after selected your layer, #2 asks "You are Joining: Points to _" (Points, Lines, Polygons, whatever).
It then gives you two options, the first to give a summary of numeric attributes, and the second option says all the points will be given the attributes of the layer being joined to it, plus it will calculate a distance field. Bingo!
From there, once you run the spatial join your new shapefile will have the distances in the attribute table. Keep in mind if you don't project your shapefiles the distances might show up in decimal degrees though. I hope this helps!