-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathdistancelib.py
160 lines (129 loc) · 4.86 KB
/
distancelib.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created on Thu Jul 5 14:34:54 2012
@author: mag
"""
import datetime
__author__ = 'Alexander Myasoedov'
__email__ = 'mag@rshu.ru'
__created__ = datetime.datetime(2012, 7, 5)
__modified__ = datetime.datetime(2012, 7, 5)
__version__ = "1.0"
__status__ = "Development"
# Adapted from code & formulas by David Z. Creemer and others
# http://www.zachary.com/blog/2005/01/12/python_zipcode_geo-programming
# http://williams.best.vwh.net/avform.htm
#
from math import sin,atan,acos,asin,atan2,sqrt,pi,radians,degrees, modf
from numpy import mod, cos
# At the equator / on another great circle???
nauticalMilePerLat = 60.00721
nauticalMilePerLongitude = 60.10793
rad = pi / 180.0
milesPerNauticalMile = 1.15078
kmsPerNauticalMile = 1.85200
degreeInMiles = milesPerNauticalMile * 60
degreeInKms = kmsPerNauticalMile * 60
# Semi-axes of WGS-84 geoidal reference
WGS84_a = 6378.1370 # Major semiaxis [km]
WGS84_b = 6356.7523 # Minor semiaxis [km]
# Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
def WGS84EarthRadius(lat):
# http://en.wikipedia.org/wiki/Earth_radius
lat = lat * rad
An = WGS84_a*WGS84_a * cos(lat)
Bn = WGS84_b*WGS84_b * sin(lat)
Ad = WGS84_a * cos(lat)
Bd = WGS84_b * sin(lat)
return sqrt( (An*An + Bn*Bn)/(Ad*Ad + Bd*Bd) )
def getDistanceLL(loc1, loc2):
"Haversine formula - give coordinates as (lat_decimal,lon_decimal) tuples"
lat1, lon1 = loc1
lat2, lon2 = loc2
# earth's mean radius = 6,371km
earthradius = WGS84EarthRadius((lat1+lat2)/2)
#if type(loc1[0]) == type(()):
# # convert from DMS to decimal
# lat1,lon1 = DMSToDecimal(loc1[0]),DMSToDecimal(loc1[1])
#if type(loc2[0]) == type(()):
# lat2,lon2 = DMSToDecimal(loc2[0]),DMSToDecimal(loc2[1])
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2.0))**2
c = 2.0 * atan2(sqrt(a), sqrt(1.0-a))
km = earthradius * c
return km
def bearing(loc1, loc2):
"Haversine formula - give coordinates as (lat_decimal,lon_decimal) tuples"
lat1, lon1 = loc1
lat2, lon2 = loc2
# earth's mean radius = 6,371km
earthradius = WGS84EarthRadius((lat1+lat2)/2)
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
bearing = atan2(sin(dlon) * cos(lat2), cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dlon))
bearing = mod(degrees(bearing) + 360, 360)
return bearing
def getPixelResolution(lat, lon, shape, units="km"):
"""
get the resolution of the pixel in Row/Col directions
"""
if str(lat.size) == '4':
col = ([lat[0], lon[0]], [lat[2], lon[2]])
row = ([lat[0], lon[0]], [lat[3], lon[3]])
else:
col = ([lat[0,0], lon[0,0]], [lat[0,-1], lon[0,-1]])
row = ([lat[0,0], lon[0,0]], [lat[-1,0], lon[-1,0]])
PxResRow = getDistanceLL(col[0], col[1])
PxResCol = getDistanceLL(row[0], row[1])
PxResRow = PxResRow/shape[1]
PxResCol = PxResCol/shape[0]
if units == "deg":
PxResRow = min(getCoordinateDiffForDistance(lat[0], lon[0], PxResRow, units="km"))
PxResCol = min(getCoordinateDiffForDistance(lat[0], lon[0], PxResCol, units="km"))
return PxResCol, PxResRow
def getDistancePx(loc1, loc2, lat, lon, shape):
"""Get the distance between two points given coordinates in loc[row,col]"""
PxResCol, PxResRow = getPixelResolution(lat, lon, shape)
km = sqrt( (PxResCol*(loc1[0]-loc2[0]))**2 + (PxResRow*(loc1[1]-loc2[1]))**2 )
return km
def DecimalToDMS(decimalvalue):
"convert a decimal value to degree,minute,second tuple"
m,s = divmod(dd*3600,60)
d,m = divmod(mnt,60)
# d = modf(decimalvalue)[0]
# m=0
# s=0
return (d,m,s)
def DMSToDecimal((degrees,minutes,seconds)):
"Convert a value from decimal (float) to degree,minute,second tuple"
d = abs(degrees) + (minutes/60.0) + (seconds/3600.0)
if degrees < 0:
return -d
else:
return d
def getCoordinateDiffForDistance(originlat, originlon, distance, units="km"):
"""return longitude & latitude values that, when added to & subtraced from
origin longitude & latitude, form a cross / 'plus sign' whose ends are
a given distance from the origin"""
degreelength = 0
if units == "km":
degreelength = degreeInKms
elif units == "miles":
degreelength = degreeInMiles
else:
raise Exception("Units must be either 'km' or 'miles'!")
lat = distance / degreelength
lon = min(distance / (cos(originlat * rad) * degreelength))
return (lat, lon)
def isWithinDistance(origin, loc, distance):
"boolean for checking whether a location is within a distance"
if getDistance(origin, loc) <= distance:
return True
else:
return False