Sunday, 14 January 2018

shapefile - Calculating distance between latitude and longitude points using python


I have a database of 2000 points, which is in CSV format containing latitude and longitude. I want to find the distance between all these points, and then show these points in a map by using pyQGIS; I have seen one example using Fiona but I want to do this in pyQGIS.


I have used following code:


def offset(iterable):

prev = None
for elem in iterable:
yield prev, elem
prev = elem

import csv
import math

def haversine(lon1, lat1, lon2, lat2):
"""

Calculate the great circle distance between two points
on the earth (specified in decimal degrees).
Source: http://gis.stackexchange.com/a/56589/15183
"""
# 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

with open('1.csv', 'rb') as f1:
reader1 = csv.reader(f1)
header1 = reader1.next()
with open('2.csv', 'rb') as f2:
reader2 = csv.reader(f2)
header2 = reader2.next()

for row1 in offset(header1):
for row2 in offset(header2):
floats1 = map(float, row1[1:])
floats2 = map(float, row2[1:])
print(floats1)
print(floats2)
print haversine(floats1[1],floats1[0],floats2[1],floats2[0])

Showing error:


Traceback (most recent call last):

File "", line 13, in
IndexError: list index out of range

My .csv file contains following data with about 2000 data pointt:


ID  LAT         LONG
1 12.953578 77.592453
2 12.953511 77.592445
3 12.953145 77.593147
4 12.951835 77.594612

Answer




pseudocode for your problem would be:


for point1 in csv:
for point2 in csv:
distance = haversine(point1, point2)

where haversine is defined as (from e.g. https://gis.stackexchange.com/a/56589/15183) :


def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees).

Source: https://gis.stackexchange.com/a/56589/15183
"""
# 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

...so the full code would be:


import csv
import math

def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees).

Source: https://gis.stackexchange.com/a/56589/15183
"""
# 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

with open('1.csv', 'rb') as f1:
reader1 = csv.reader(f1, delimiter=';', quoting=csv.QUOTE_NONE)
header1 = reader1.next()
with open('2.csv', 'rb') as f2:
reader2 = csv.reader(f2, delimiter=';', quoting=csv.QUOTE_NONE)
header2 = reader2.next()
for row1 in reader1:
for row2 in reader2:

floats1 = map(float, row1[1:])
floats2 = map(float, row2[1:])
print(floats1)
print(floats2)
print haversine(floats1[1],floats1[0],floats2[1],floats2[0])

with the two files 1.csv:


name;lat;lon
NYC;40.58;-74.03


and 2.csv:


name;lat;lon
Cape Town;-33.90;18.46

...result checked against http://www.movable-type.co.uk/scripts/latlong.html seems correct.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...