如何找到一条直线与网格的交点？

Question

我有轨迹数据，其中每个轨迹由一系列坐标（x，y点）组成，每个轨迹都有一个唯一的ID。

这些轨迹在x-y平面上，我想将整个平面划分为相等大小的网格（正方形网格）。这个网格显然是看不见的，但用于将轨迹划分为子段。每当轨迹与网格线相交时，它就在那里被分割，并成为具有新ID的新子轨迹。

我已经包含了一个简单的手工图表，以清楚说明我的期望。

可以看到轨迹在网格线的交点处被分割，每个分段都有一个新的唯一标识符。我正在使用Python进行工作，并寻求关于此的Python实现链接、建议、算法或伪代码。如果有任何不清楚的地方，请告诉我。更新：为了将平面划分为网格，单元格索引如下所示：

#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 

由于飞机现在被分成了25个单元格，每个单元格都有一个cellid。为了找到交点，也许我可以检查轨迹中的下一个坐标，如果cellid保持不变，则该轨迹段位于同一单元格中，并且与网格没有交点。例如，如果x_cellid[2]大于x_cellid[0]，则该段与垂直网格线相交。尽管如此，我仍然不确定如何找到与网格线相交并将轨迹在交点处分段的方法并给它们新的id。

- Liza

---

段1和12怎么办..段的开始和结束可以在任何地方吗？ - Hariom Singh

是的，轨迹可以在任何地方，因此第一个和最后一个段落也可能是如此。 - Liza

Answer 1

这可以通过使用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);

结果:

如果不使用shapely，可以自己实现：

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);

- HYRY

---

非常感谢，这很完美，但我有一个小问题，有没有一种方法可以不使用Shapely来完成相同的操作？我正在使用另一个与Shapely不兼容的库。 - Liza

与哪个库不兼容？编写新代码是最后的手段。 - hamster on wheels

@Liza，我添加了一个只使用numpy的解决方案。 - HYRY

Answer 2

每个轨迹由一系列直线段组成。因此，您需要编写一种例程，将每个线段分解为完全位于网格单元内的部分。这样一个例程的基础将是数字差分分析（DDA）算法，但您需要修改基本算法，因为您需要每个单元格内的线段端点，而不仅仅是访问哪些单元格。

有几件事情需要注意：

1）如果您使用浮点数，请注意计算步骤值时的舍入误差，因为这可能会导致算法失败。出于这个原因，许多人选择转换为整数网格，显然会失去精度。这 是关于该问题的很好的讨论，并附有一些工作代码（但不是python）。

2）您需要决定围绕单元格的4条网格线中哪些属于该单元格。一种约定是使用底部和左侧边缘。如果考虑水平线段落在网格线上的情况，则可以看到该问题 - 它的线段属于上面的单元格还是下面的单元格？

干杯

- RaffleBuffle

Answer 3

你提出了很多要求。一旦你有了一个大致的方案，你应该自己攻击设计和编码的大部分。在Stack Overflow上进行算法识别是合理的；但请求设计和参考链接不是。

我建议您将点坐标放入列表中。使用NumPy和SciKit的能力来插值网格交叉点。您可以将线段存储在列表中（任何定义数据设计中线段的方式都可以）。考虑制作一个字典，使您可以通过网格坐标检索线段。例如，如果线段仅由端点表示，而且点是您的类别，您可能会像这样做，使用每个正方形的左下角作为其定义点：

grid_seg = {
    (0.5, 0.5): [p0, p1],
    (1.0, 0.5): [p1, p2],
    (1.0, 1.0): [p2, p3],
    ...
}

p0、p1等是插值交叉点。

- Prune

Answer 4

data = list of list of coordinates
For point_id, point_coord in enumerate(point_coord_list):
   if current point & last point stayed in same cell:
        append point's index to last list of data
   else:
        append a new empty list to data
        interpolate the two points and add a new point
        that is on the grid lines.

数据存储所有轨迹。数据中的每个列表都是一条轨迹。

通过将点的坐标除以单元格的维度，然后四舍五入为整数，可以找到沿x和y轴的单元格索引（x_cell_id，y_cell_id）。如果当前点的单元格索引与上一个点的单元格索引相同，则这两个点在同一个单元格中。列表适合插入新点，但不如数组内存效率高。

可能创建一个轨迹类是个好主意。或者如果坐标列表浪费太多内存，可以使用内存缓冲区和稀疏数据结构代替列表和数组来存储x-y坐标。向数组中插入新点很慢，因此我们可以使用另一个数组来存储新点。

警告：我没有深思熟虑下面的事情。它可能存在错误，需要有人填补空白。

# coord       n x 2 numpy array. 
#             columns 0, 1 are x and y coordinate. 
#             row n is for point n
# cell_size   length of one side of the square cell.
# n_ycells    number of cells along the y axis

import numpy as np
cell_id_2d = (coord / cell_size).astype(int)
x_cell_id = cell_id_2d[:,0]
y_cell_id = cell_id_2d[:,1]
cell_id_1d = x_cell_id + y_cell_id*n_x_cells

# if the trajectory exits a cell, its cell id changes
# and the delta_cell_id is not zero.
delta_cell_id = cell_id_1d[1:] - cell_id_1d[:-1]

# The nth trajectory should contains the points from
# the (crossing_id[n])th to the (crossing_id[n + 1] - 1)th
w = np.where(delta_cell_id != 0)[0]
crossing_ids = np.empty(w.size + 1)
crossing_ids[1:] = w
crossing_ids[0] = 0

# need to interpolate when the trajectory cross cell boundary.
# probably can replace this loop with numpy functions/indexing
new_points = np.empty((w.size, 2))
for i in range(1, n):
    st = coord[crossing_ids[i]]
    en = coord[crossing_ids[i+1]]
    # 1. check which boundary of the cell is crossed
    # 2. interpolate
    # 3. put points into new_points

# Each trajectory contains some points from coord array and 2 points 
# in the new_points array.

为了检索，创建一个稀疏数组，其中包含coord数组中起始点的索引。

如果单元格大小较大，线性插值可能会看起来不好。

进一步解释：

网格描述

For n_xcells = 4, n_ycells = 3, the grid is:

   0   1   2   3   4
0 [  ][  ][  ][  ][  ]
1 [  ][  ][  ][* ][  ]
2 [  ][  ][  ][  ][  ]

[* ] has an x_index of 3 and a y_index of 1.

网格中有 (n_x_cells * n_y_cells) 个单元格。

点和单元格之间的关系

包含轨迹第i个点的单元格具有 x_index 为 x_cell_id[i]，y_index 为 x_cell_id[i]。这是通过将点的 xy 坐标除以单元格长度进行离散化，然后截断为整数得到的。

The cell_id_1d of the cells are the number in [  ]

   0   1   2   3   4
0 [0 ][1 ][2 ][3 ][4 ]
1 [5 ][6 ][7 ][8 ][9 ]
2 [10][11][12][13][14]

cell_id_1d[i] = x_cell_id[i] + y_cell_id[i]*n_x_cells

我将第i个点的细胞索引对(x_cell_id[i], y_cell_id[i])转换为一个称为cell_id_1d的单一索引。

如何找出轨迹在第i个点是否退出细胞

现在，如果且仅当(x_cell_id[i], y_cell_id[i]) == (x_cell_id[i + 1], y_cell_id[i + 1])且cell_id_1d[i] == cell_id[i + 1]，以及cell_id[i + 1] - cell_id[i] == 0时，第i个点和第(i + 1)个点位于同一单元格中。delta_cell_ids[i] = cell_id_1d[i + 1] - cell_id[i]，当且仅当第i个点和(i + 1)个点在同一单元格中时，它为零。

- hamster on wheels

---

非常感谢，让我来处理它，如果需要进一步的帮助或澄清，我会尽快更新给您。 - Liza

你真的应该学习一些关于numpy库的基础知识。唉...你需要这个做什么？ - hamster on wheels