python - How to find intersection of a line with a mesh?

Question

I have trajectory data, where each trajectory consists of a sequence of coordinates(x, y points) and each trajectory is identified by a unique ID.

These trajectories are in?x - y?plane, and I want to divide the whole plane into equal sized grid (square grid). This grid is obviously invisible but is used to divide trajectories into sub-segments. Whenever a trajectory intersects with a grid line, it is segmented there and becomes a new sub-trajectory with new_id.

I have included a simple handmade graph to make clear what I am expecting.

It can be seen how the trajectory is divided at the intersections of the grid lines, and each of these segments has new unique id.

I am working on Python, and seek some python implementation links, suggestions, algorithms, or even a pseudocode for the same.

Please let me know if anything is unclear.

UPDATE

In order to divide the plane into grid , cell indexing is done as following:

#finding cell id for each coordinate
#cellid = (coord / cellSize).astype(int)
cellid = (coord / 0.5).astype(int)
cellid
Out[] : array([[1, 1],
              [3, 1],
              [4, 2],
              [4, 4],
              [5, 5],
              [6, 5]])
#Getting x-cell id and y-cell id separately 
x_cellid = cellid[:,0]
y_cellid = cellid[:,1]

#finding total number of cells
xmax = df.xcoord.max()
xmin = df.xcoord.min()
ymax = df.ycoord.max()
ymin = df.ycoord.min()
no_of_xcells = math.floor((xmax-xmin)/ 0.5)
no_of_ycells = math.floor((ymax-ymin)/ 0.5)
total_cells = no_of_xcells * no_of_ycells
total_cells
Out[] : 25 

Since the plane is now divided into 25 cells each with a cellid. In order to find intersections, maybe I could check the next coordinate in the trajectory, if the cellid remains the same, then that segment of the trajectory is in the same cell and has no intersection with grid. Say, if x_cellid[2] is greater than x_cellid[0], then segment intersects vertical grid lines. Though, I am still unsure how to find the intersections with the grid lines and segment the trajectory on intersections giving them new id.

Answer 1

This can be solved by shapely:

%matplotlib inline
import pylab as pl
from shapely.geometry import MultiLineString, LineString
import numpy as np
from matplotlib.collections import LineCollection

x0, y0, x1, y1 = -10, -10, 10, 10
n = 11

lines = []
for x in np.linspace(x0, x1, n):
    lines.append(((x, y0), (x, y1)))

for y in np.linspace(y0, y1, n):
    lines.append(((x0, y), (x1, y)))

grid = MultiLineString(lines)

x = np.linspace(-9, 9, 200)
y = np.sin(x)*x
line = LineString(np.c_[x, y])

fig, ax = pl.subplots()
for i, segment in enumerate(line.difference(grid)):
    x, y = segment.xy
    pl.plot(x, y)
    pl.text(np.mean(x), np.mean(y), str(i))

lc = LineCollection(lines, color="gray", lw=1, alpha=0.5)
ax.add_collection(lc);

The result:

To not use shapely, and do it yourself:

import pylab as pl
import numpy as np
from matplotlib.collections import LineCollection

x0, y0, x1, y1 = -10, -10, 10, 10
n = 11
xgrid = np.linspace(x0, x1, n)
ygrid = np.linspace(y0, y1, n)
x = np.linspace(-9, 9, 200)
y = np.sin(x)*x
t = np.arange(len(x))

idx_grid, idx_t = np.where((xgrid[:, None] - x[None, :-1]) * (xgrid[:, None] - x[None, 1:]) <= 0)
tx = idx_t + (xgrid[idx_grid] - x[idx_t]) / (x[idx_t+1] - x[idx_t])

idx_grid, idx_t = np.where((ygrid[:, None] - y[None, :-1]) * (ygrid[:, None] - y[None, 1:]) <= 0)
ty = idx_t + (ygrid[idx_grid] - y[idx_t]) / (y[idx_t+1] - y[idx_t])

t2 = np.sort(np.r_[t, tx, tx, ty, ty])

x2 = np.interp(t2, t, x)
y2 = np.interp(t2, t, y)

loc = np.where(np.diff(t2) == 0)[0] + 1

xlist = np.split(x2, loc)
ylist = np.split(y2, loc)


fig, ax = pl.subplots()
for i, (xp, yp) in enumerate(zip(xlist, ylist)):
    pl.plot(xp, yp)
    pl.text(np.mean(xp), np.mean(yp), str(i))


lines = []
for x in np.linspace(x0, x1, n):
    lines.append(((x, y0), (x, y1)))

for y in np.linspace(y0, y1, n):
    lines.append(((x0, y), (x1, y)))

lc = LineCollection(lines, color="gray", lw=1, alpha=0.5)
ax.add_collection(lc);