Match two sets of latitude/longitude points by distance

latlongmatch.py

#script that attempts to match compare lat,long points two csv files, match by distance
#justin grimes (@justgrimes) 04/16/12

import math
import csv
import sys
import re

# takes two lat/long points and returns distance, modified code from -> http://www.johndcook.com/python_longitude_latitude.html
def distance_on_unit_sphere(lat1, long1, lat2, long2):

    # convert to 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
        
    cos = (math.sin(phi1)*math.sin(phi2)*math.cos(theta1 - theta2) + 
           math.cos(phi1)*math.cos(phi2))
	# some reason cos can be greater than 1.0 in formula?? 
	# causes math error cos(>1) is NaaN, so if estimate higher than 1.0 assume distance 0
    if cos > 1.0: 
		cos = 1.0
    arc = math.acos(cos)*3960 #multiple by 3960 for miles, 6373 for kilometers
    return arc #returns distance

#open files
a = csv.reader(open('imls.csv','rb')) #open data file #1 ie imls.csv
z = csv.writer(open('output.csv','wb'))
z.writerow(['id1','lat1','long1','id2','lat2','long2','distance']) #header
for row in a: #iterate through file
	a_id = row[0] #unique ids in first row
	a_lat = float(row[1]) # extract lat
	a_long = float(row[2]) # extract long
	h = [] #create/wipe array, maybe better way to wipe/reset?
	b = csv.reader(open('ntia.csv','rb')) #open data file #2 ie ntia.csv
	for row in b:
		b_id = row[0] #unique ids in first row
		b_lat = float(row[1]) # extract lat
		b_long = float(row[2]) # extract long
		distance = distance_on_unit_sphere(a_lat,a_long,b_lat,b_long)
		h.append([distance,b_id,b_lat,b_long])
	h.sort() #sort to find shortest distance, assumes shorts distance is best match
	r = h.pop(0) #extract shortest info
	r = str(r).strip("[") #strip [ 
	r = str(r).strip("]") #strip ]
	dis, b_id, b_lat, b_long = re.split(',',r) #separate data back out
	z.writerow([a_id,a_lat,a_long,b_id,b_lat,b_long,dis]) #write data to file
	
# useful to quick check results http://www.movable-type.co.uk/scripts/latlong.html

# more precise numbers for urth's radius
# 6378137 meters
# 6378.137 km
# 3963.191 miles
# 3441.596 nautical miles