From 55675f90006a4df448fc16f4f782b819c943ffb8 Mon Sep 17 00:00:00 2001
From: Matt Davis <jiffyclub@gmail.com>
Date: Fri, 19 Dec 2014 13:05:24 -0800
Subject: [PATCH] New function for calculating great circle distance

Calculated using the Haversine formula.

I also removed some of the existing functions in pandana.utils
that depended on GeoPandas and other geo libraries.
---
 pandana/tests/test_great_circle_dist.py | 16 +++++
 pandana/utils.py                        | 89 +++++++++++--------------
 2 files changed, 54 insertions(+), 51 deletions(-)
 create mode 100644 pandana/tests/test_great_circle_dist.py

diff --git a/pandana/tests/test_great_circle_dist.py b/pandana/tests/test_great_circle_dist.py
new file mode 100644
index 00000000..c34b35b7
--- /dev/null
+++ b/pandana/tests/test_great_circle_dist.py
@@ -0,0 +1,16 @@
+import numpy.testing as npt
+
+from pandana.utils import great_circle_dist as gcd
+
+
+def test_gcd():
+    # tested against geopy
+    # https://geopy.readthedocs.org/en/latest/#module-geopy.distance
+    lat1 = 41.49008
+    lon1 = -71.312796
+    lat2 = 41.499498
+    lon2 = -81.695391
+
+    expected = 864456.76162966
+
+    npt.assert_allclose(gcd(lat1, lon1, lat2, lon2), expected)
diff --git a/pandana/utils.py b/pandana/utils.py
index ca78628d..97086573 100644
--- a/pandana/utils.py
+++ b/pandana/utils.py
@@ -1,51 +1,38 @@
-from shapely.geometry import Point
-from fiona import crs
-import geopandas as gpd
-import pandas as pd
-import numpy as np
-
-
-def bbox_convert(bbox, from_epsg, to_epsg):
-    bbox = gpd.GeoSeries([Point(bbox[0], bbox[1]),
-                          Point(bbox[2], bbox[3])],
-                         crs=crs.from_epsg(from_epsg))
-    bbox = bbox.to_crs(epsg=to_epsg)
-    bbox = [bbox[0].x, bbox[0].y, bbox[1].x, bbox[1].y]
-    return bbox
-
-
-def get_nodes_from_osm(bbox, query, to_epsg=3740):
-    gdf = gpd.io.osm.query_osm('node',
-                               bbox=bbox,
-                               tags=query)
-    gdf = gdf[gdf.type == 'Point'].to_crs(epsg=to_epsg)
-    print "Found %d nodes" % len(gdf)
-    x, y = zip(*[(p.x, p.y) for (i, p)
-               in gdf.geometry.iteritems()])
-    x = pd.Series(x)
-    y = pd.Series(y)
-    return x, y
-
-
-def anything_score(net, config, max_distance, decay, bbox):
-    score = pd.Series(np.zeros(len(net.node_ids)), index=net.node_ids)
-    for query, weights in config.iteritems():
-        print "Computing for query: %s" % query
-        print "Fetching nodes from OSM"
-        x, y = get_nodes_from_osm(bbox, query)
-        print "Done"
-        net.set_pois(query, x, y)
-        print "Computing nearest"
-        df = net.nearest_pois(max_distance, query, num_pois=len(weights))
-        print "Done"
-
-        for idx, weight in enumerate(weights):
-            # want the 1st not the 0th
-            idx += 1
-            print "Adding contribution %f for number %d nearest" % \
-                  (weight, idx)
-            score += decay(df[idx])*weight
-            # print score.describe()
-
-    assert score.min() > 0
-    return score/score.max()*100
+from __future__ import division
+
+import math
+
+
+def great_circle_dist(lat1, lon1, lat2, lon2):
+    """
+    Get the distance (in meters) between two lat/lon points
+    via the Haversine formula.
+
+    Parameters
+    ----------
+    lat1, lon1, lat2, lon2 : float
+        Latitude and longitude in degrees.
+
+    Returns
+    -------
+    dist : float
+        Distance in meters.
+
+    """
+    radius = 6372795  # meters
+
+    lat1 = math.radians(lat1)
+    lon1 = math.radians(lon1)
+    lat2 = math.radians(lat2)
+    lon2 = math.radians(lon2)
+
+    dlat = lat2 - lat1
+    dlon = lon2 - lon1
+
+    # formula from:
+    # http://en.wikipedia.org/wiki/Haversine_formula#The_haversine_formula
+    a = math.pow(math.sin(dlat / 2), 2)
+    b = math.cos(lat1) * math.cos(lat2) * math.pow(math.sin(dlon / 2), 2)
+    d = 2 * radius * math.asin(math.sqrt(a + b))
+
+    return d