#! /usr/bin/python -tt # vim: set sw=4 sts=4 et tw=80 fileencoding=utf-8: # # edist - quick calculations for distances between two points # Copyright (C) 2007 James Rowe; # All rights reserved. # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA. # __version__ = "0.1.0" __date__ = "2007-05-04" __author__ = "U{James Rowe }" __copyright__ = "Copyright (C) 2007 James Rowe" __license__ = "GNU General Public License Version 2" __credits__ = "" __history__ = "See Mercurial repository" __doc__ = """ edist - quick calculations for distances between two points @version: %s @author: %s @copyright: %s @status: complete @license: %s """ % (__version__, __author__, __copyright__, __license__) import math def to_dms(angle): """ Convert decimal angle to degrees, minutes and seconds notation >>> to_dms(52.015) (True, 52, 0, 54) >>> to_dms(-0.221) (False, 0, 13, 15) @type angle: C{float} or coercible to C{float} @param angle: Angle to convert @rtype: C{tuple} of C{bool} for sign and C{int}s for values @return: Angle converted to degrees, minutes and seconds """ if angle < 0: sign = False angle = abs(angle) else: sign = True degrees = math.floor(angle) minutes = math.floor((angle - degrees) * 60) seconds = math.floor((angle - degrees - minutes / 60) * 3600) return sign, int(degrees), int(minutes), int(seconds) class Point: """ Simple class for representing a location on a sphere @ivar latitude: Location's latitude @ivar longitude: Locations's longitude """ def __init__(self, latitude=0, longitude=0): """ Initialise a new Point object @type latitude: C{float} or coercible to C{float} @param latitude: Location's latitude @type longitude: C{float} or coercible to C{float} @param longitude: Location's longitude """ self.latitude = float(latitude) self.longitude = float(longitude) def __repr__(self): """ Self-documenting string representation >>> Point(52.015, -0.221) Point(52.015000, -0.221000) @rtype: C{str} @return: String to recreate Point object """ return "Point(%f, %f)" % (self.latitude, self.longitude) def __str__(self, mode="decimal"): """ Pretty printed location string @todo: Require Python v2.5+, and this could be a I{whole lot} cleaner using conditional expressions. >>> print Point(52.015, -0.221) N52.015°; W0.221° >>> print Point(52.015, -0.221).__str__(mode="dms") 052°00'54"N, 000°13'15"W @type mode: C{str} @param mode: Coordinate formatting system to use @rtype: C{str} @return: Human readable string representation of Point object """ if mode == "decimal": if self.latitude < 0: text = "S" else: text = "N" text += "%s°; " % abs(self.latitude) if self.longitude < 0: text += "W" else: text += "E" text += "%s°" % abs(self.longitude) elif mode == "dms": latitude_dms = to_dms(self.latitude) longitude_dms = to_dms(self.longitude) format = '''%03i°%02i'%02i"''' text = format % latitude_dms[1:] if latitude_dms[0]: text += "N, " else: text += "S, " text += format % longitude_dms[1:] if longitude_dms[0]: text += "W" else: text += "W" else: raise ValueError("Unknown mode type `%s'" % mode) return text def as_rad(self): """ Calculate latitude and longitude in radians >>> "%.5f %.5f" % Point(52.015, -0.221).as_rad() '0.90783 -0.00386' @rtype: C{tuple} of C{float}s @return: Latitude and longitude in radians """ return math.radians(self.latitude), math.radians(self.longitude) def distance(self, other, method="haversine"): """ Calculate the distance from self to other >>> "%i kM" % Point(52.015, -0.221).distance(Point(52.6333, -2.5)) '169 kM' >>> "%i kM" % Point(52.015, -0.221).distance(Point(52.6333, -2.5), ... method="sloc") '169 kM' As a smoke test this check uses the example from Wikipedia's U{Great-circle distance entry } of Nashville International Airport to Los Angeles International Airport, and is correct to within 2 kilometres of the calculation there. >>> "%i kM" % Point(36.1200, -86.6700).distance(Point(33.9400, ... -118.4000)) '2885 kM' >>> "%i kM" % Point(36.1200, -86.6700).distance(Point(33.9400, ... -118.4000), ... method="sloc") '2885 kM' @type other: C{Point} instance @param other: Location to calculate distance to @type method: C{str} @param method: Method used to calculate distance @rtype: C{float} @return: Distance in kilometres between self and other """ self_latitude, self_longitude = self.as_rad() other_latitude, other_longitude = other.as_rad() longitude_difference = other_longitude - self_longitude latitude_difference = other_latitude - self_latitude earth_radius = 6370 #kM if method == "haversine": temp = math.sin(latitude_difference / 2) ** 2 + \ math.cos(self_latitude) * \ math.cos(other_latitude) * \ math.sin(longitude_difference / 2) ** 2 return 2 * earth_radius * math.atan2(math.sqrt(temp), math.sqrt(1-temp)) elif method == "sloc": return math.acos(math.sin(self_latitude) * \ math.sin(other_latitude) + \ math.cos(self_latitude) * \ math.cos(other_latitude) * \ math.cos(longitude_difference)) * earth_radius else: raise ValueError("Unknown method type `%s'" % method) def bearing(self, other): """ Calculate the initial bearing from self to other >>> "%i" % Point(52.015, -0.221).bearing(Point(52.6333, -2.5)) '294' >>> "%i" % Point(52.6333, -2.5).bearing(Point(52.015, -0.221)) '113' >>> "%i" % Point(36.1200, -86.6700).bearing(Point(33.9400, ... -118.4000)) '274' >>> "%i" % Point(33.9400, -118.4000).bearing(Point(36.1200, ... -86.6700)) '76' @type other: C{Point} instance @param other: Location to calculate bearing to @rtype: C{float} @return: Initial bearing from self to other in degrees """ self_latitude, self_longitude = self.as_rad() other_latitude, other_longitude = other.as_rad() longitude_difference = other_longitude - self_longitude y = math.sin(longitude_difference) * math.cos(other_latitude) x = math.cos(self_latitude) * math.sin(other_latitude) - \ math.sin(self_latitude) * math.cos(other_latitude) * \ math.cos(longitude_difference) bearing = math.degrees(math.atan2(y,x)) # Always return positive North-aligned bearing return (bearing + 360) % 360 if __name__ == '__main__': import doctest import sys sys.exit(doctest.testmod()[0])