diff --git a/client/renderer/interface.yaml b/client/renderer/interface.yaml index b2917ca..97c7675 100644 --- a/client/renderer/interface.yaml +++ b/client/renderer/interface.yaml @@ -103,6 +103,26 @@ - { name: range, type: "integer[]" } - { name: options, type: "text[]", default: 'ARRAY[]::text[]' } +- name: cdb_mapzen_isochrone + return_type: SETOF cdb_dataservices_client.isoline + multi_row: true + multi_field: true + params: + - { name: source, type: "geometry(Geometry, 4326)" } + - { name: mode, type: text } + - { name: range, type: "integer[]" } + - { name: options, type: "text[]", default: 'ARRAY[]::text[]' } + +- name: cdb_mapzen_isodistance + return_type: SETOF cdb_dataservices_client.isoline + multi_row: true + multi_field: true + params: + - { name: source, type: "geometry(Geometry, 4326)" } + - { name: mode, type: text } + - { name: range, type: "integer[]" } + - { name: options, type: "text[]", default: 'ARRAY[]::text[]' } + - name: cdb_route_point_to_point return_type: cdb_dataservices_client.simple_route multi_field: true diff --git a/server/extension/sql/15_config_helper.sql b/server/extension/sql/15_config_helper.sql index 3d6fd72..118bbec 100644 --- a/server/extension/sql/15_config_helper.sql +++ b/server/extension/sql/15_config_helper.sql @@ -26,6 +26,20 @@ RETURNS boolean AS $$ return True $$ LANGUAGE plpythonu SECURITY DEFINER; +CREATE OR REPLACE FUNCTION cdb_dataservices_server._get_mapzen_isolines_config(username text, orgname text) +RETURNS boolean AS $$ + cache_key = "user_mapzen_isolines_routing_config_{0}".format(username) + if cache_key in GD: + return False + else: + from cartodb_services.metrics import MapzenIsolinesRoutingConfig + plpy.execute("SELECT cdb_dataservices_server._connect_to_redis('{0}')".format(username)) + redis_conn = GD["redis_connection_{0}".format(username)]['redis_metadata_connection'] + mapzen_isolines_config = MapzenIsolinesRoutingConfig(redis_conn, plpy, username, orgname) + GD[cache_key] = mapzen_isolines_config + return True +$$ LANGUAGE plpythonu SECURITY DEFINER; + CREATE OR REPLACE FUNCTION cdb_dataservices_server._get_internal_geocoder_config(username text, orgname text) RETURNS boolean AS $$ cache_key = "user_internal_geocoder_config_{0}".format(username) diff --git a/server/extension/sql/80_isolines_helper.sql b/server/extension/sql/80_isolines_helper.sql index b598cf6..f04e62c 100644 --- a/server/extension/sql/80_isolines_helper.sql +++ b/server/extension/sql/80_isolines_helper.sql @@ -53,3 +53,74 @@ RETURNS SETOF cdb_dataservices_server.isoline AS $$ finally: quota_service.increment_total_service_use() $$ LANGUAGE plpythonu SECURITY DEFINER; + + +CREATE OR REPLACE FUNCTION cdb_dataservices_server._cdb_mapzen_isolines( + username TEXT, + orgname TEXT, + isotype TEXT, + source geometry(Geometry, 4326), + mode TEXT, + data_range integer[], + options text[]) +RETURNS SETOF cdb_dataservices_server.isoline AS $$ + import json + from cartodb_services.mapzen import MatrixClient + from cartodb_services.mapzen import MapzenIsolines + from cartodb_services.metrics import QuotaService + + redis_conn = GD["redis_connection_{0}".format(username)]['redis_metrics_connection'] + user_mapzen_isolines_routing_config = GD["user_mapzen_isolines_routing_config_{0}".format(username)] + + # -- Check the quota + quota_service = QuotaService(user_mapzen_isolines_routing_config, redis_conn) + if not quota_service.check_user_quota(): + plpy.error('You have reached the limit of your quota') + + try: + client = MatrixClient(user_mapzen_isolines_routing_config.mapzen_matrix_api_key) + mapzen_isolines = MapzenIsolines(client) + + if source: + lat = plpy.execute("SELECT ST_Y('%s') AS lat" % source)[0]['lat'] + lon = plpy.execute("SELECT ST_X('%s') AS lon" % source)[0]['lon'] + origin = {'lat': lat, 'lon': lon} + else: + raise Exception('source is NULL') + + # -- TODO Support options properly + isolines = {} + if isotype == 'isodistance': + for r in data_range: + isoline = mapzen_isolines.calculate_isodistance(origin, mode, r) + isolines[r] = (isoline) + elif isotype == 'isochrone': + for r in data_range: + isoline = mapzen_isolines.calculate_isochrone(origin, mode, r) + isolines[r] = (isoline) + + result = [] + for r in data_range: + + # -- TODO encapsulate this block into a func/method + locations = isolines[r] + [ isolines[r][0] ] # close the polygon repeating the first point + wkt_coordinates = ','.join(["%f %f" % (l['lon'], l['lat']) for l in locations]) + sql = "SELECT ST_MPolyFromText('MULTIPOLYGON((({0})))', 4326) as geom".format(wkt_coordinates) + multipolygon = plpy.execute(sql, 1)[0]['geom'] + + result.append([source, r, multipolygon]) + + quota_service.increment_success_service_use() + quota_service.increment_isolines_service_use(len(isolines)) + return result + except BaseException as e: + import sys, traceback + type_, value_, traceback_ = sys.exc_info() + quota_service.increment_failed_service_use() + error_msg = 'There was an error trying to obtain isolines using mapzen: {0}'.format(e) + plpy.debug(traceback.format_tb(traceback_)) + raise e + #plpy.error(error_msg) + finally: + quota_service.increment_total_service_use() +$$ LANGUAGE plpythonu SECURITY DEFINER; diff --git a/server/extension/sql/85_isodistance.sql b/server/extension/sql/85_isodistance.sql index 3c3ac0e..6349015 100644 --- a/server/extension/sql/85_isodistance.sql +++ b/server/extension/sql/85_isodistance.sql @@ -19,3 +19,23 @@ RETURNS SETOF cdb_dataservices_server.isoline AS $$ return isolines $$ LANGUAGE plpythonu; + +-- mapzen isodistance +CREATE OR REPLACE FUNCTION cdb_dataservices_server.cdb_mapzen_isodistance(username TEXT, orgname TEXT, source geometry(Geometry, 4326), mode TEXT, range integer[], options text[] DEFAULT array[]::text[]) +RETURNS SETOF cdb_dataservices_server.isoline AS $$ + plpy.execute("SELECT cdb_dataservices_server._connect_to_redis('{0}')".format(username)) + redis_conn = GD["redis_connection_{0}".format(username)]['redis_metrics_connection'] + plpy.execute("SELECT cdb_dataservices_server._get_mapzen_isolines_config({0}, {1})".format(plpy.quote_nullable(username), plpy.quote_nullable(orgname))) + user_isolines_config = GD["user_mapzen_isolines_routing_config_{0}".format(username)] + type = 'isodistance' + + mapzen_plan = plpy.prepare("SELECT cdb_dataservices_server._cdb_mapzen_isolines($1, $2, $3, $4, $5, $6, $7) as isoline; ", ["text", "text", "text", "geometry(Geometry, 4326)", "text", "integer[]", "text[]"]) + result = plpy.execute(mapzen_plan, [username, orgname, type, source, mode, range, options]) + isolines = [] + for element in result: + isoline = element['isoline'] + isoline = isoline.translate(None, "()").split(',') + isolines.append(isoline) + + return isolines +$$ LANGUAGE plpythonu; diff --git a/server/extension/sql/90_isochrone.sql b/server/extension/sql/90_isochrone.sql index 6237bd7..11f11cb 100644 --- a/server/extension/sql/90_isochrone.sql +++ b/server/extension/sql/90_isochrone.sql @@ -19,3 +19,24 @@ RETURNS SETOF cdb_dataservices_server.isoline AS $$ return isolines $$ LANGUAGE plpythonu; + + +-- mapzen isochrones +CREATE OR REPLACE FUNCTION cdb_dataservices_server.cdb_mapzen_isochrone(username TEXT, orgname TEXT, source geometry(Geometry, 4326), mode TEXT, range integer[], options text[] DEFAULT array[]::text[]) +RETURNS SETOF cdb_dataservices_server.isoline AS $$ + plpy.execute("SELECT cdb_dataservices_server._connect_to_redis('{0}')".format(username)) + redis_conn = GD["redis_connection_{0}".format(username)]['redis_metrics_connection'] + plpy.execute("SELECT cdb_dataservices_server._get_mapzen_isolines_config({0}, {1})".format(plpy.quote_nullable(username), plpy.quote_nullable(orgname))) + user_isolines_config = GD["user_mapzen_isolines_routing_config_{0}".format(username)] + type = 'isochrone' + + mapzen_plan = plpy.prepare("SELECT cdb_dataservices_server._cdb_mapzen_isolines($1, $2, $3, $4, $5, $6, $7) as isoline; ", ["text", "text", "text", "geometry(Geometry, 4326)", "text", "integer[]", "text[]"]) + result = plpy.execute(mapzen_plan, [username, orgname, type, source, mode, range, options]) + isolines = [] + for element in result: + isoline = element['isoline'] + isoline = isoline.translate(None, "()").split(',') #--TODO what is this for? + isolines.append(isoline) + + return isolines +$$ LANGUAGE plpythonu; diff --git a/server/lib/python/cartodb_services/cartodb_services/mapzen/__init__.py b/server/lib/python/cartodb_services/cartodb_services/mapzen/__init__.py index ee7552c..e3dca35 100644 --- a/server/lib/python/cartodb_services/cartodb_services/mapzen/__init__.py +++ b/server/lib/python/cartodb_services/cartodb_services/mapzen/__init__.py @@ -1,2 +1,4 @@ from routing import MapzenRouting, MapzenRoutingResponse +from isolines import MapzenIsolines from geocoder import MapzenGeocoder +from matrix_client import MatrixClient diff --git a/server/lib/python/cartodb_services/cartodb_services/mapzen/isolines.py b/server/lib/python/cartodb_services/cartodb_services/mapzen/isolines.py new file mode 100644 index 0000000..de75459 --- /dev/null +++ b/server/lib/python/cartodb_services/cartodb_services/mapzen/isolines.py @@ -0,0 +1,160 @@ +from math import cos, sin, tan, sqrt, pi, radians, degrees, asin, atan2 +import logging + +class MapzenIsolines: + + NUMBER_OF_ANGLES = 24 + MAX_ITERS = 5 + TOLERANCE = 0.1 + + EARTH_RADIUS_METERS = 6367444 + + def __init__(self, matrix_client): + self._matrix_client = matrix_client + + """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 + #logging.basicConfig(level=logging.DEBUG, filename='/tmp/isolines.log') + #logging.basicConfig(level=logging.DEBUG) + logging.debug('origin = %s' % origin) + logging.debug('costing_model = %s' % costing_model) + logging.debug('isorange = %d' % 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 + 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 + + logging.debug('i = %d, costs = %s' % (i, costs)) + + 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) + } diff --git a/server/lib/python/cartodb_services/cartodb_services/mapzen/matrix_client.py b/server/lib/python/cartodb_services/cartodb_services/mapzen/matrix_client.py new file mode 100644 index 0000000..9b4df53 --- /dev/null +++ b/server/lib/python/cartodb_services/cartodb_services/mapzen/matrix_client.py @@ -0,0 +1,43 @@ +import requests +import json + +class MatrixClient: + + """ + A minimal client for Mapzen Time-Distance Matrix Service + + Example: + + client = MatrixClient('your_api_key') + locations = [{"lat":40.744014,"lon":-73.990508},{"lat":40.739735,"lon":-73.979713},{"lat":40.752522,"lon":-73.985015},{"lat":40.750117,"lon":-73.983704},{"lat":40.750552,"lon":-73.993519}] + costing = 'pedestrian' + + print client.one_to_many(locations, costing) + """ + + ONE_TO_MANY_URL = 'https://matrix.mapzen.com/one_to_many' + + def __init__(self, matrix_key): + self._matrix_key = matrix_key + + """Get distances and times to a set of locations. + See https://mapzen.com/documentation/matrix/api-reference/ + + Args: + locations Array of {lat: y, lon: x} + costing Costing model to use + + Returns: + A dict with one_to_many, units and locations + """ + def one_to_many(self, locations, costing): + request_params = { + 'json': json.dumps({'locations': locations}), + 'costing': costing, + 'api_key': self._matrix_key + } + response = requests.get(self.ONE_TO_MANY_URL, params=request_params) + + response.raise_for_status() # raise exception if not 200 OK + + return response.json() diff --git a/server/lib/python/cartodb_services/cartodb_services/metrics/__init__.py b/server/lib/python/cartodb_services/cartodb_services/metrics/__init__.py index 9bba57c..78d09e7 100644 --- a/server/lib/python/cartodb_services/cartodb_services/metrics/__init__.py +++ b/server/lib/python/cartodb_services/cartodb_services/metrics/__init__.py @@ -1,3 +1,3 @@ -from config import GeocoderConfig, MapzenGeocoderConfig, IsolinesRoutingConfig, InternalGeocoderConfig, RoutingConfig, ConfigException, ObservatorySnapshotConfig, ObservatoryConfig +from config import GeocoderConfig, MapzenGeocoderConfig, IsolinesRoutingConfig, MapzenIsolinesRoutingConfig, InternalGeocoderConfig, RoutingConfig, ConfigException, ObservatorySnapshotConfig, ObservatoryConfig from quota import QuotaService from user import UserMetricsService diff --git a/server/lib/python/cartodb_services/cartodb_services/metrics/config.py b/server/lib/python/cartodb_services/cartodb_services/metrics/config.py index d1b4a35..7ac7dba 100644 --- a/server/lib/python/cartodb_services/cartodb_services/metrics/config.py +++ b/server/lib/python/cartodb_services/cartodb_services/metrics/config.py @@ -241,6 +241,40 @@ class IsolinesRoutingConfig(ServiceConfig): return self._geocoder_type == self.GOOGLE_GEOCODER +class MapzenIsolinesRoutingConfig(ServiceConfig): + + PERIOD_END_DATE = 'period_end_date' + + def __init__(self, redis_connection, db_conn, username, orgname=None): + super(MapzenIsolinesRoutingConfig, self).__init__(redis_connection, db_conn, + username, orgname) + try: + self._mapzen_matrix_api_key = self._db_config.mapzen_matrix_api_key + self._isolines_quota = self._db_config.mapzen_matrix_monthly_quota + self._period_end_date = date_parse(self._redis_config[self.PERIOD_END_DATE]) + except Exception as e: + raise ConfigException("Malformed config for Mapzen isolines: {0}".format(e)) + + @property + def service_type(self): + return 'mapzen_isolines' + + @property + def isolines_quota(self): + return self._isolines_quota + + @property + def soft_isolines_limit(self): + return False + + @property + def period_end_date(self): + return self._period_end_date + + @property + def mapzen_matrix_api_key(self): + return self._mapzen_matrix_api_key + class InternalGeocoderConfig(ServiceConfig): def __init__(self, redis_connection, db_conn, username, orgname=None): @@ -438,6 +472,8 @@ class ServicesDBConfig: raise ConfigException('Mapzen configuration missing') else: mapzen_conf = json.loads(mapzen_conf_json) + self._mapzen_matrix_api_key = mapzen_conf['matrix']['api_key'] + self._mapzen_matrix_quota = mapzen_conf['matrix']['monthly_quota'] self._mapzen_routing_api_key = mapzen_conf['routing']['api_key'] self._mapzen_routing_quota = mapzen_conf['routing']['monthly_quota'] self._mapzen_geocoder_api_key = mapzen_conf['geocoder']['api_key'] @@ -492,6 +528,14 @@ class ServicesDBConfig: def heremaps_geocoder_cost_per_hit(self): return self._heremaps_geocoder_cost_per_hit + @property + def mapzen_matrix_api_key(self): + return self._mapzen_matrix_api_key + + @property + def mapzen_matrix_monthly_quota(self): + return self._mapzen_matrix_quota + @property def mapzen_routing_api_key(self): return self._mapzen_routing_api_key diff --git a/server/lib/python/cartodb_services/cartodb_services/metrics/quota.py b/server/lib/python/cartodb_services/cartodb_services/metrics/quota.py index c627769..4496e3c 100644 --- a/server/lib/python/cartodb_services/cartodb_services/metrics/quota.py +++ b/server/lib/python/cartodb_services/cartodb_services/metrics/quota.py @@ -70,6 +70,9 @@ class QuotaChecker: elif re.match('here_isolines', self._user_service_config.service_type) is not None: return self.__check_isolines_quota() + elif re.match('mapzen_isolines', + self._user_service_config.service_type) is not None: + return self.__check_isolines_quota() elif re.match('routing_mapzen', self._user_service_config.service_type) is not None: return self.__check_routing_quota() diff --git a/server/lib/python/cartodb_services/test/__init__.py b/server/lib/python/cartodb_services/test/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/server/lib/python/cartodb_services/test/test_helper.py b/server/lib/python/cartodb_services/test/test_helper.py index b0ffaf0..fe6fd8f 100644 --- a/server/lib/python/cartodb_services/test/test_helper.py +++ b/server/lib/python/cartodb_services/test/test_helper.py @@ -65,7 +65,7 @@ def _plpy_execute_side_effect(*args, **kwargs): if args[0] == "SELECT cartodb.CDB_Conf_GetConf('heremaps_conf') as conf": return [{'conf': '{"geocoder": {"app_id": "app_id", "app_code": "code", "geocoder_cost_per_hit": 1}, "isolines": {"app_id": "app_id", "app_code": "code"}}'}] elif args[0] == "SELECT cartodb.CDB_Conf_GetConf('mapzen_conf') as conf": - return [{'conf': '{"routing": {"api_key": "api_key_rou", "monthly_quota": 1500000}, "geocoder": {"api_key": "api_key_geo", "monthly_quota": 1500000}}'}] + return [{'conf': '{"routing": {"api_key": "api_key_rou", "monthly_quota": 1500000}, "geocoder": {"api_key": "api_key_geo", "monthly_quota": 1500000}, "matrix": {"api_key": "api_key_mat", "monthly_quota": 1500000}}'}] elif args[0] == "SELECT cartodb.CDB_Conf_GetConf('logger_conf') as conf": return [{'conf': '{"geocoder_log_path": "/dev/null"}'}] elif args[0] == "SELECT cartodb.CDB_Conf_GetConf('data_observatory_conf') as conf": diff --git a/server/lib/python/cartodb_services/test/test_mapzengeocoder.py b/server/lib/python/cartodb_services/test/test_mapzengeocoder.py index 0c7be54..7eee6c0 100644 --- a/server/lib/python/cartodb_services/test/test_mapzengeocoder.py +++ b/server/lib/python/cartodb_services/test/test_mapzengeocoder.py @@ -11,7 +11,7 @@ requests_mock.Mocker.TEST_PREFIX = 'test_' @requests_mock.Mocker() -class GoogleGeocoderTestCase(unittest.TestCase): +class MapzenGeocoderTestCase(unittest.TestCase): MAPZEN_GEOCODER_URL = 'https://search.mapzen.com/v1/search' EMPTY_RESPONSE = """{ diff --git a/server/lib/python/cartodb_services/test/test_mapzenisolines.py b/server/lib/python/cartodb_services/test/test_mapzenisolines.py new file mode 100644 index 0000000..4fb4e96 --- /dev/null +++ b/server/lib/python/cartodb_services/test/test_mapzenisolines.py @@ -0,0 +1,76 @@ +import unittest +from cartodb_services.mapzen import MapzenIsolines +from math import radians, cos, sin, asin, sqrt + +""" +This file is basically a sanity test on the algorithm. + + +It uses a mocked client, which returns the cost based on a very simple model: +just proportional to the distance from origin to the target point. +""" + +class MatrixClientMock(): + + def __init__(self, speed): + """ + Sets up the mock with a speed in km/h + """ + self._speed = speed + + def one_to_many(self, locations, costing): + origin = locations[0] + distances = [self._distance(origin, l) for l in locations] + response = { + 'one_to_many': [ + [ + { + 'distance': distances[i] * self._speed, + 'time': distances[i] / self._speed * 3600, + 'to_index': i, + 'from_index': 0 + } + for i in xrange(0, len(distances)) + ] + ], + 'units': 'km', + 'locations': [ + locations + ] + } + return response + + def _distance(self, a, b): + """ + Calculate the great circle distance between two points + on the earth (specified in decimal degrees) + http://stackoverflow.com/questions/4913349/haversine-formula-in-python-bearing-and-distance-between-two-gps-points + + Returns: + distance in meters + """ + + # convert decimal degrees to radians + lon1, lat1, lon2, lat2 = map(radians, [a['lon'], a['lat'], b['lon'], b['lat']]) + + # haversine formula + dlon = lon2 - lon1 + dlat = lat2 - lat1 + a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 + c = 2 * asin(sqrt(a)) + r = 6371 # Radius of earth in kilometers. Use 3956 for miles + return c * r + + +class MapzenIsolinesTestCase(unittest.TestCase): + + def setUp(self): + speed = 4 # in km/h + matrix_client = MatrixClientMock(speed) + self.mapzen_isolines = MapzenIsolines(matrix_client) + + def test_calculate_isochrone(self): + origin = {"lat":40.744014,"lon":-73.990508} + transport_mode = 'walk' + isorange = 10 * 60 # 10 minutes + solution = self.mapzen_isolines.calculate_isochrone(origin, transport_mode, isorange)