160 lines
6.4 KiB
Python
160 lines
6.4 KiB
Python
from math import cos, sin, tan, sqrt, pi, radians, degrees, asin, atan2
|
|
|
|
|
|
class MapzenIsolines:
|
|
|
|
NUMBER_OF_ANGLES = 24
|
|
MAX_ITERS = 5
|
|
TOLERANCE = 0.1
|
|
|
|
EARTH_RADIUS_METERS = 6367444
|
|
|
|
def __init__(self, matrix_client, logger):
|
|
self._matrix_client = matrix_client
|
|
self._logger = logger
|
|
|
|
"""Get an isochrone using mapzen API.
|
|
|
|
The implementation tries to sick close to the SQL API:
|
|
cdb_isochrone(source geometry, mode text, range integer[], [options text[]]) -> SETOF isoline
|
|
|
|
But this calculates just one isoline.
|
|
|
|
Args:
|
|
origin dict containing {lat: y, lon: x}
|
|
transport_mode string, for the moment just "car" or "walk"
|
|
isorange int range of the isoline in seconds
|
|
|
|
Returns:
|
|
Array of {lon: x, lat: y} as a representation of the isoline
|
|
"""
|
|
def calculate_isochrone(self, origin, transport_mode, time_range):
|
|
if transport_mode == 'walk':
|
|
max_speed = 3.3333333 # In m/s, assuming 12km/h walking speed
|
|
costing_model = 'pedestrian'
|
|
elif transport_mode == 'car':
|
|
max_speed = 41.67 # In m/s, assuming 140km/h max speed
|
|
costing_model = 'auto'
|
|
else:
|
|
raise NotImplementedError('car and walk are the only supported modes for the moment')
|
|
|
|
upper_rmax = max_speed * time_range # an upper bound for the radius
|
|
|
|
return self.calculate_isoline(origin, costing_model, time_range, upper_rmax, 'time')
|
|
|
|
"""Get an isodistance using mapzen API.
|
|
|
|
Args:
|
|
origin dict containing {lat: y, lon: x}
|
|
transport_mode string, for the moment just "car" or "walk"
|
|
isorange int range of the isoline in seconds
|
|
|
|
Returns:
|
|
Array of {lon: x, lat: y} as a representation of the isoline
|
|
"""
|
|
def calculate_isodistance(self, origin, transport_mode, distance_range):
|
|
if transport_mode == 'walk':
|
|
costing_model = 'pedestrian'
|
|
elif transport_mode == 'car':
|
|
costing_model = 'auto'
|
|
else:
|
|
raise NotImplementedError('car and walk are the only supported modes for the moment')
|
|
|
|
upper_rmax = distance_range # an upper bound for the radius, going in a straight line
|
|
|
|
return self.calculate_isoline(origin, costing_model, distance_range, upper_rmax, 'distance', 1000.0)
|
|
|
|
|
|
"""Get an isoline using mapzen API.
|
|
|
|
The implementation tries to sick close to the SQL API:
|
|
cdb_isochrone(source geometry, mode text, range integer[], [options text[]]) -> SETOF isoline
|
|
|
|
But this calculates just one isoline.
|
|
|
|
Args:
|
|
origin dict containing {lat: y, lon: x}
|
|
costing_model string "auto" or "pedestrian"
|
|
isorange int Range of the isoline in seconds
|
|
upper_rmax float An upper bound for the binary search
|
|
cost_variable string Variable to optimize "time" or "distance"
|
|
unit_factor float A factor to adapt units of isorange (meters) and units of distance (km)
|
|
|
|
Returns:
|
|
Array of {lon: x, lat: y} as a representation of the isoline
|
|
"""
|
|
def calculate_isoline(self, origin, costing_model, isorange, upper_rmax, cost_variable, unit_factor=1.0):
|
|
|
|
# NOTE: not for production
|
|
# self._logger.debug('Calculate isoline', data={"origin": origin, "costing_model": costing_model, "isorange": isorange})
|
|
|
|
# Formally, a solution is an array of {angle, radius, lat, lon, cost} with cardinality NUMBER_OF_ANGLES
|
|
# we're looking for a solution in which abs(cost - isorange) / isorange <= TOLERANCE
|
|
|
|
# Initial setup
|
|
angles = self._get_angles(self.NUMBER_OF_ANGLES) # array of angles
|
|
rmax = [upper_rmax] * self.NUMBER_OF_ANGLES
|
|
rmin = [0.0] * self.NUMBER_OF_ANGLES
|
|
location_estimates = [self._calculate_dest_location(origin, a, upper_rmax / 2.0) for a in angles]
|
|
|
|
# Iterate to refine the first solution
|
|
for i in xrange(0, self.MAX_ITERS):
|
|
# Calculate the "actual" cost for each location estimate.
|
|
# NOTE: sometimes it cannot calculate the cost and returns None.
|
|
# Just assume isorange and stop the calculations there
|
|
|
|
response = self._matrix_client.one_to_many([origin] + location_estimates, costing_model)
|
|
costs = [None] * self.NUMBER_OF_ANGLES
|
|
if not response:
|
|
# In case the matrix client doesn't return any data
|
|
break
|
|
|
|
for idx, c in enumerate(response['one_to_many'][0][1:]):
|
|
if c[cost_variable]:
|
|
costs[idx] = c[cost_variable]*unit_factor
|
|
else:
|
|
costs[idx] = isorange
|
|
|
|
errors = [(cost - isorange) / float(isorange) for cost in costs]
|
|
max_abs_error = max([abs(e) for e in errors])
|
|
if max_abs_error <= self.TOLERANCE:
|
|
# good enough, stop there
|
|
break
|
|
|
|
# let's refine the solution, binary search
|
|
for j in xrange(0, self.NUMBER_OF_ANGLES):
|
|
|
|
if abs(errors[j]) > self.TOLERANCE:
|
|
if errors[j] > 0:
|
|
rmax[j] = (rmax[j] + rmin[j]) / 2.0
|
|
else:
|
|
rmin[j] = (rmax[j] + rmin[j]) / 2.0
|
|
|
|
location_estimates[j] = self._calculate_dest_location(origin, angles[j], (rmax[j]+rmin[j])/2.0)
|
|
|
|
# delete points that got None
|
|
location_estimates_filtered = []
|
|
for i, c in enumerate(costs):
|
|
if c != isorange:
|
|
location_estimates_filtered.append(location_estimates[i])
|
|
|
|
return location_estimates_filtered
|
|
|
|
|
|
|
|
# NOTE: all angles in calculations are in radians
|
|
def _get_angles(self, number_of_angles):
|
|
step = (2.0 * pi) / number_of_angles
|
|
return [(x * step) for x in xrange(0, number_of_angles)]
|
|
|
|
def _calculate_dest_location(self, origin, angle, radius):
|
|
origin_lat_radians = radians(origin['lat'])
|
|
origin_long_radians = radians(origin['lon'])
|
|
dest_lat_radians = asin(sin(origin_lat_radians) * cos(radius / self.EARTH_RADIUS_METERS) + cos(origin_lat_radians) * sin(radius / self.EARTH_RADIUS_METERS) * cos(angle))
|
|
dest_lng_radians = origin_long_radians + atan2(sin(angle) * sin(radius / self.EARTH_RADIUS_METERS) * cos(origin_lat_radians), cos(radius / self.EARTH_RADIUS_METERS) - sin(origin_lat_radians) * sin(dest_lat_radians))
|
|
|
|
return {
|
|
'lon': degrees(dest_lng_radians),
|
|
'lat': degrees(dest_lat_radians)
|
|
}
|