In this very short post I want to point you to some code for calculating the centroid and distance to that centroid for a set of points in numpy. As bonus I also included some profiling. Note that I already blogged about the centroid function in a previous post.

from math import atan2, sqrt, degrees
import numpy as np
from math import radians, sin, cos
RADIUS = 6371.009
def get_centroid(points):
xy = np.asarray(points)
xy = np.radians(xy)
lon, lat = xy[:, 0], xy[:, 1]
avg_x = np.sum(np.cos(lat) * np.cos(lon)) / xy.shape[0]
avg_y = np.sum(np.cos(lat) * np.sin(lon)) / xy.shape[0]
avg_z = np.sum(np.sin(lat)) / xy.shape[0]
center_lon = atan2(avg_y, avg_x)
hyp = sqrt(avg_x * avg_x + avg_y * avg_y)
center_lat = atan2(avg_z, hyp)
return degrees(center_lon), degrees(center_lat)
def gc_distance_points(a, points):
b = np.asarray(points)
lat1, lng1 = radians(a[1]), radians(a[0])
lat2, lng2 = np.radians(b[:, 1]), np.radians(b[:, 0])
sin_lat1, cos_lat1 = sin(lat1), cos(lat1)
sin_lat2, cos_lat2 = np.sin(lat2), np.cos(lat2)
delta_lng = np.subtract(lng2, lng1)
cos_delta_lng, sin_delta_lng = np.cos(delta_lng), np.sin(delta_lng)
d = np.arctan2(np.sqrt((cos_lat2 * sin_delta_lng) ** 2 +
(cos_lat1 * sin_lat2 -
sin_lat1 * cos_lat2 * cos_delta_lng) ** 2),
sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng)
return RADIUS * d
def gc_dist(a, b):
from math import radians, sin, cos, sqrt, atan2
lat1, lng1 = radians(a[1]), radians(a[0])
lat2, lng2 = radians(b[1]), radians(b[0])
sin_lat1, cos_lat1 = sin(lat1), cos(lat1)
sin_lat2, cos_lat2 = sin(lat2), cos(lat2)
delta_lng = lng2 - lng1
cos_delta_lng, sin_delta_lng = cos(delta_lng), sin(delta_lng)
d = atan2(sqrt((cos_lat2 * sin_delta_lng) ** 2 +
(cos_lat1 * sin_lat2 -
sin_lat1 * cos_lat2 * cos_delta_lng) ** 2),
sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng)
return RADIUS * d
if __name__ == '__main__':
import random
random.seed(42)
pts = [(random.uniform(-180, 180), random.uniform(-90, 90)) for _ in range(1000000)]
centr = get_centroid(pts)
import cProfile
cProfile.runctx('gc_distance_points(centr, pts)', globals(), locals())
cProfile.runctx('dists = [gc_dist(centr, b) for b in pts]', globals(), locals())

All code is also available at Github. Feel free to post a comment or email in case you have any questions.