Recently needed to calculate the overlap area of many circles and thought that Numpy would be perfect.

Imagine we have a long CSV that looks a bit like this (units are Miles):

distance radiusA radiusB id
7.70 3.40 4.70 555

Where the example filled in looks like this:

stores_circles

We can simply wrap up this formula in a function using Numpy:

inline41

import numpy as np

in_csv = '.../test_circle_overlaps_py.csv'
out_csv = '.../out_circle_overlaps_py.csv'

def circle_area_overlap(dist, rada, radb):
    """
    http://mathworld.wolfram.com/Circle-CircleIntersection.html
    Given two circles in 2D euclidean space find their intersection area
    """
    a = (rada**2)*np.arccos(((dist**2)+(rada**2)-(radb**2))/(2*dist*rada))
    b = (radb**2)*np.arccos(((dist**2)+(radb**2)-(rada**2))/(2*dist*radb))
    c = (0.5*((-dist+rada+radb)*(dist+rada-radb)*(dist-rada+radb)*(dist+rada+radb))**0.5)
    return a + b - c

# Import data
try_calcs = np.genfromtxt(in_csv,
                          delimiter=',',
                          names=True,
                          usecols=(0, 1, 2, 3),
                          dtype=(float, float, float, float))
# Run
py_calc = circle_area_overlap(try_calcs['distance'],
                              try_calcs['radiusA'],
                              try_calcs['radiusB'])
# Convert NaN to 0 and combine with ID
out = np.column_stack((try_calcs['id'],  # ID
                       np.nan_to_num(py_calc)))
# Save the calculations
np.savetxt(out_csv, out, fmt='%8i,%8f')

To find the figure to be around 0.65 miles squared. If we wanted to do a rough-sense check (and make sure we aren’t 100x off) we could draw the intersection and check the area using Google Maps (which gives us 1.92 km^2 = 0.74 M^2)

inst.png

Edit:

Without doing anything too complicated (Fortran/C) below are some timings for the three basic cases (python list, python numpy, numexpr):

import numpy as np
from math import acos
import numexpr as ne

def circle_area_overlap(dist, rada, radb):
    """
    http://mathworld.wolfram.com/Circle-CircleIntersection.html
    Given two circles in 2D euclidean space find their intersection area
    """
    try:
        a = (rada**2)*acos(((dist**2)+(rada**2)-(radb**2))/(2*dist*rada))
        b = (radb**2)*acos(((dist**2)+(radb**2)-(rada**2))/(2*dist*radb))
        c = (0.5*((-dist+rada+radb)*(dist+rada-radb)*(dist-rada+radb)*(dist+rada+radb))**0.5)
    except ValueError:
        # Circles don't overlap
        return 0
    return a + b - c

def np_circle_area_overlap(dist, rada, radb):
    a = (rada**2)*np.arccos(((dist**2)+(rada**2)-(radb**2))/(2*dist*rada))
    b = (radb**2)*np.arccos(((dist**2)+(radb**2)-(rada**2))/(2*dist*radb))
    c = (0.5*((-dist+rada+radb)*(dist+rada-radb)*(dist-rada+radb)*(dist+rada+radb))**0.5)
    return a + b - c

def numpexpr_overlap(dist, rada, radb):
    a = "(rada**2)*arccos(((dist**2)+(rada**2)-(radb**2))/(2*dist*rada))"
    b = "(radb**2)*arccos(((dist**2)+(radb**2)-(rada**2))/(2*dist*radb))"
    c = "(0.5*((-dist+rada+radb)*(dist+rada-radb)*(dist-rada+radb)*(dist+rada+radb))**0.5)"
    str_eval = "%s + %s - %s" % (a,b,c)
    print(str_eval)
    return ne.evaluate(str_eval)

# Create test data
# Assume columns are: distance, radius A, radius B
try_calcs = np.random.rand(100*1000*1000,3) # Numpy array
py_try_calcs = try_calcs.tolist() # Python list

"""
Timings
"""
# Python list comprehension with math
# 4 min 13 seconds
%time py_calc_a = [circle_area_overlap(row[0], row[1], row[2]) for row in py_try_calcs]
print(py_calc_a[:10])

# Numpy
# 31.8 seconds
%time py_calc_b = np_circle_area_overlap(try_calcs[:,0],try_calcs[:,1],try_calcs[:,2])
print(py_calc_b[:10])

# Numexpr
# 9.28 seconds
%time py_calc_c = numpexpr_overlap(try_calcs[:,0],try_calcs[:,1],try_calcs[:,2])
print(py_calc_c[:10])

Curious to try wrapping a Fortran function …