Merge pull request #216 from CartoDB/mapzen-isolines-reloaded

Add mapzen isolines

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

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)
+        }

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()

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

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

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()

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":

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 = """{

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)