CODEWITHSUNDEEP
numpyinterpolation
flotzilla
flotzilla

Reputation: 1252

Numpy: find mean coordinate of points along line

I have a bunch of points in a 2D space which all reside on a line (polygon). How can I compute the mean coordinate of these points on the line?

I don't mean the centroid of the points in the 2D space (as @rth initially proposed in his answer), but the mean location of the points along the line on which they reside. So basically, I could transform the line to a 1D axis, compute the mean location in 1D, and transform the location of the mean back into the 2D space.

Maybe these are exactly the necessary steps, but I think (or hope) that there is a function in numpy/scipy which allows me to do this in one step.

Upvotes: 1

Views: 1684

Answers (2)

flotzilla
flotzilla

Reputation: 1252

Meanwhile I've found the answer to my question, although using Shapely instead of Numpy.

from shapely.geometry import LineString, Point

# lists of points as (x,y) tuples
path_xy = [...]
points_xy = [...] # should be on or near path

path = LineString(path_xy)            # create path object
pts = [Point(p) for p in points_xy]   # create point objects
dist = [path.project(p) for p in pts] # distances along path
mean_dist = np.mean(dist)             # mean distance along path
mean = path.interpolate(mean_dist)    # mean point

mean_xy = (mean.x,mean.y)

This works perfectly!

(That's is also why I have to accept it as the answer, though I highly appreciate @rth's help!)

Upvotes: 1

rth
rth

Reputation: 11201

Edit: The approach you describe in the question is indeed probably the simplest way for solving this problem.

Here is an implementation that calculates the positions of vertices along the line in 1D, takes their mean, and finally calculates the corresponding 2D position with parametric interpolation,

import numpy as np
from scipy.interpolate import splprep, splev

vert = np.random.randn(1000, 2) # vertices definition here

# calculate the Euclidean distances between consecutive vertices
# equivalent to a for loop with
# dl[i] = ((vert[i+1, 0] - vert[i, 0])**2 + (vert[i+1,1] - vert[i,1])**2)**0.5
dl = (np.diff(vert, axis=0)**2).sum(axis=1)**0.5 

# pad with 0, so dl.shape[0] == vert.shape[0] for convenience
dl = np.insert(dl, 0, 0.0)
l = np.cumsum(dl) # 1D coordinates along the line
l_mean = np.mean(l) # mean in the line coordinates

# calculate the coordinate of l_mean in 2D space
# with parametric B-spline interpolation
tck, _ = splprep(x=vert.T,  u=l, k=3)
res = splev(l_mean, tck)
print(res)

Edit2: Assuming now that you have a high resolution set of points for your path vert_full and some approximate measurements vert_1, vert_2, etc, what you could do is the following.

  • Project each points of vert_1, etc. onto the exact path. Assuming that vert_full has much more datapoints than vert_1, we can simply look for the nearest neighbours of vert_1 in vert_full:

    from scipy.spatial import cKDTree
tr = cKDTree(vert_full)
d, idx = tr.query(vert_1, k=1)
vert_1_proj = vert_full[idx] # this gives the projected corrdinates onto vert_full
# I have not actually run this, so it might require minor changes

  • Use the above mean calculation with the new vert_1_proj vector.

Upvotes: 1

Related Questions