Calculate average spacing per direction I am interested in finding the average spacing between objects on a directional basis. I did my best to minimally implement this in Python and to conform to PEP8.

I try to explain what I'm doing in the code comments. The problem I want to solve is to find the average spacing between circles of known radius and locations (in my application, these locations are irregular) on a directional basis. For example, what is the average spacing between these circles in the x, diagonal, and y directions?



import numpy as np

# I have a set of x and y coordinates describing the
# locations of objects
px = np.array([0, 0, 0, 10, 10, 10, 20, 20, 20])
py = np.array([0, 10, 20, 0, 10, 20, 0, 10, 20])

# Each object is a circle with known radius size
point_diameter = 2

# I'll put my answers here
spacing = {}

# The challenge is to calculate the average spacing per point on a directional
#  basis, e.g., what is the average distance between points in the direction
#  of 45 degrees?


# To do this, I define a function to check if points are alligned in a
#  direction
def check_alligned(direction, p1x, p1y, p2x, p2y):

    # Since the points have some wisdth, I need to check the "edge cases"
    # as well as the centerlines.
    tolerance = 0.2
    pr = point_diameter / 2
    for offsetx1 in [-1, 0, 1]:
        for offsetx2 in [-1, 0, 1]:
            for offsety1 in [-1, 0, 1]:
                for offsety2 in [-1, 0, 1]:

                    # The scheme is to check if the direction leads to the object x or y location. 
                    # I chose or to handle cases where the x or y difference is zero.
                    distance = distance_between([p1x + pr*offsetx1,
                                                     p1y + pr*offsety1],
                                                    [p2x + pr*offsetx2,
                                                     p2y + pr*offsety2])
                    if abs(p1x + distance * np.cos(np.radians(direction)) - p2x) &lt;= tolerance or abs(p1y + distance * np.sin(np.radians(direction)) - p2y) &lt;= tolerance:
                         return True
    return False


def distance_between(v1, v2):
    return np.sqrt( (v1[0] - v2[0])**2 + (v1[1] - v2[1])**2)

# Initial Smallest Spacing - use to find if we get to a spacing in the loop
# and record the smallest one
ISP = 9e4

# Iterate through directions to calculate average distance in each direction.
# Directions go counterclockwise and x points in the direction of 0 degrees.
for direction in [0, 60, 45, 90]:

    # These are a running sum and counter of 1st nearest neighbor spacings
    # that are less than an initial smallest spacing value. They are used
    # to calculate the average spacing in each direction.
    total_spacing = 0.
    space_count = 0

    # Look through every point
    for ii in range(px.size):

        smallest_spacing = ISP
        # and every other point except ones we've already checked
        for jj in range(px.size - ii - 1):

            # If the jj point is not in the same line as the
            # ii point along this diretion, skip it
            idx2 = ii + jj + 1

            if check_alligned(direction, px[ii], py[ii], px[idx2], py[idx2]):
                space = np.sqrt((px[ii]-px[idx2])**2 +
                                          (py[ii] - py[idx2])**2)
                if space &lt; smallest_spacing:
                    smallest_spacing = space

        # Record a distance if it is less than the initial
        # smallest spacing
        if smallest_spacing != ISP:
            total_spacing += smallest_spacing
            space_count += 1

    # record the average spacing in this direction
    if space_count:
        spacing[direction] = total_spacing / space_count
    else:
        spacing[direction] = 0.

print spacing
I don't fully understand your problem statement, could you elaborate? I asked about your definition of distance because unless the distance you are interested in is some projection of the euclidian distance, it should be constant regardless which direction you are looking.

In any case, I did your problem manually and found a different result than your program above reports. So here's my implementation with inline comments:

#!/usr/local/bin/python3

from sklearn.neighbors import NearestNeighbors
import numpy as np
import matplotlib.pyplot as plt

class Turbine():

    def __init__(self, x, y):
        self.initialx = x
        self.initialy = y
        self.x = x
        self.y = y

    def rotateBy(self, theta):
        newX = self.x*np.cos(theta)-self.y*np.sin(theta)
        newY = self.x*np.sin(theta)+self.y*np.cos(theta)
        self.x = newX
        self.y = newY

    def point(self):
        return [self.x, self.y]

class Farm():
    def __init__(self, turbineList):
        self.turbineList = turbineList

    def rotateFarmBy(self, theta):
        for turbine in self.turbineList:
            turbine.rotateBy(theta)

    def turbineCoordinates(self):
        return np.array([turbine.point() for turbine in self.turbineList])

    def turbineCount(self):
        return len(self.turbineList)


for theta in [0, np.pi/4., np.pi/3., np.pi/2.]:

    # establish the coordinates,
    # but first test with a small subset and compare with hand calculations
    coordinates = [[0, 0], [0, 10], [0, 20]]
    # coordinates = [[0, 0], [10, 10], [20, 20]]
    # coordinates = [[0, 0], [0, 10], [0, 20], [10, 0], [10, 10], [10, 20], [20, 0], [20, 10], [20, 20]]

    # instantiate the turbines and the container farm
    turbines = [Turbine(coord[0], coord[1]) for coord in coordinates]
    farm = Farm(turbines)
    farm.rotateFarmBy(theta)

    # get the N nearest neighbors, where N in this case is the number of turbines in the farm
    coordinates = farm.turbineCoordinates()
    neighbors = NearestNeighbors(n_neighbors=farm.turbineCount(), algorithm='ball_tree').fit(coordinates)
    distances, indices = neighbors.kneighbors(coordinates)

    # remove any zero distances so that you dont count each turbine twice
    nonzeros = distances[distances &gt; 0]

    print("at {0:0.1f} degrees, the average distance is {1:0.3f}".format(theta*180/np.pi, sum(nonzeros)/farm.turbineCount()))

    # print(indices)
    # print(distances)
    # plt.scatter([c[0] for c in coordinates], [c[1] for c in coordinates])
    # plt.show()