Collection of useful distance calculation functions for the earth

geodist.py

# A collection of useful distance calculation functions for the earth
from __future__ import division
from math import sin, cos, acos, asin, radians, degrees, sqrt

EARTH_RADIUS = 6371009                      # meters


class Geolocation(object):
	def __init__(self, latitude, longitude):
		self.latitude = latitude            # Restrict to be within -90 and 90
		self.longitude = longitude          # Restrict to be within -180 and 180


def great_circle_distance_coarse(pointA, pointB):
	"""
	Computes the great circle distance in meters between two points using spherical law of cosines
	Requires: EARTH_RADIUS in meters, decimal degree point.latitude, point.longitude
	"""
	dlat = pointA.latitude - pointB.latitude
	dlon = pointA.longitude - pointB.longitude

	term1 = sin(radians(pointA.latitude)) * sin(radians(pointB.latitude))
	term2 = cos(radians(pointA.latitude)) * cos(radians(pointB.latitude)) * cos(radians(dlon))
	angle = acos(term1 + term2)             # Radians
	return EARTH_RADIUS * angle


def great_circle_distance_fine(pointA, pointB):
	"""
	Computes the great circle distance in meters between two points using the Haversine formula
	Requires: EARTH_RADIUS in meters, decimal degree point.latitude, point.longitude
	"""
	dlat = pointA.latitude - pointB.latitude
	dlon = pointA.longitude - pointB.longitude

	term1 = sin(radians(dlat)/2)**2
	term2 = cos(radians(pointA.latitude)) * cos(radians(pointB.latitude)) * sin(radians(dlon)/2)**2
	angle = 2* asin( sqrt(term1 + term2) )             # Radians
	return EARTH_RADIUS * angle


def great_circle_distance_very_fine(pointA, pointB):
	"""
	Computes the great circle distance in meters between two points using special case of Vincenty formula
	Requires: EARTH_RADIUS in meters, decimal degree point.latitude, point.longitude
	"""
	pass

# Gold Standard:  ehhh? Google Maps?


macg = Geolocation(42.355404, -71.099742)
dome = Geolocation(42.359781, -71.092316)
stata_stairs = Geolocation(42.360873, -71.09081)
capital = Geolocation(42.358148, -71.063715)

north_side_road = Geolocation(42.355761, -71.099834)
south_side_road = Geolocation(42.355715, -71.099798)

print great_circle_distance_coarse(macg, capital)
print great_circle_distance_fine(macg, capital)

print great_circle_distance_coarse(dome, stata_stairs)
print great_circle_distance_fine(dome, stata_stairs)

# real width of road is what, maybe 4 meters?
print great_circle_distance_coarse(north_side_road, south_side_road)
print great_circle_distance_fine(north_side_road, south_side_road)