如何找到线与网格的交点？
How to find intersection of a line with a mesh?
我有轨迹数据，其中每个轨迹都由一系列坐标（x，y 点）组成，每个轨迹都由唯一的 ID 标识。
这些轨迹在 x - y 平面内，我想将整个平面分成大小相等的网格（正方形网格）。该网格显然是不可见的，但用于将轨迹划分为子段。每当轨迹与网格线相交时，它就会在那里分段并成为一个新的子轨迹 new_id.
我附上了一个简单的手工图表来说明我的期望。
可以看出轨迹是如何在网格线的交点处进行划分的，并且每段都有新的唯一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。
你问的太多了。一旦你有了通用的方法，你就应该自己攻击大部分的设计和编码。算法识别对于Stack Overflow来说是合理的；要求设计和参考链接 is not.
我建议您将点坐标放入列表中。使用 NumPy 和 SciKit 功能来插入网格交叉点。您可以将段存储在列表中（在数据设计中定义段的任何内容）。考虑制作一个字典，允许您通过网格坐标检索线段。例如，如果线段仅由端点表示，并且点是您的 class，您可能会得到类似这样的东西，使用每个正方形的左下角作为其定义点：
grid_seg = {
(0.5, 0.5): [p0, p1],
(1.0, 0.5): [p1, p2],
(1.0, 1.0): [p2, p3],
...
}
其中 p0、p1 等是插值交叉点。
每条轨迹都由一系列直线段组成。因此，您需要一个例程来将每条线段分成完全位于网格单元格内的部分。这种例程的基础是 Digital Differential Analyzer (DDA) 算法，但您需要修改基本算法，因为您需要每个单元格内直线的端点，而不仅仅是访问哪些单元格。
您必须注意的几件事：
1) 如果您使用的是浮点数，请注意步长值计算中的舍入错误，因为这些可能会导致算法失败。出于这个原因，许多人选择转换为整数网格，显然会损失精度。 This 是对问题的很好的讨论，有一些工作代码（虽然不是 python）。
2) 您需要决定单元格周围的 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）可以通过将点的坐标除以单元格的维度，然后舍入为整数来找到。如果当前点的单元格索引与最后一个点的单元格索引相同，则这两个点在同一单元格中。列表适合插入新点，但它的内存效率不如数组。
为轨迹创建 class 可能是个好主意。或者如果坐标列表浪费太多内存，则使用内存缓冲区和稀疏数据结构而不是列表和列表以及 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.
为了检索，制作一个包含坐标数组中起点索引的稀疏数组。
如果像元大小很大，线性插值可能看起来很糟糕。
进一步说明：
网格描述
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。delta_cell_ids[i] = cell_id_1d[i + 1] - cell_id[i]，当且仅当第 ith 和 (i + 1)th 时为零 个点在同一个单元格中。
这可以通过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);
结果：
不要用得体，自己动手：
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);
我有轨迹数据，其中每个轨迹都由一系列坐标（x，y 点）组成，每个轨迹都由唯一的 ID 标识。
这些轨迹在 x - y 平面内，我想将整个平面分成大小相等的网格（正方形网格）。该网格显然是不可见的，但用于将轨迹划分为子段。每当轨迹与网格线相交时，它就会在那里分段并成为一个新的子轨迹 new_id.
我附上了一个简单的手工图表来说明我的期望。
可以看出轨迹是如何在网格线的交点处进行划分的，并且每段都有新的唯一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。
你问的太多了。一旦你有了通用的方法，你就应该自己攻击大部分的设计和编码。算法识别对于Stack Overflow来说是合理的；要求设计和参考链接 is not.
我建议您将点坐标放入列表中。使用 NumPy 和 SciKit 功能来插入网格交叉点。您可以将段存储在列表中（在数据设计中定义段的任何内容）。考虑制作一个字典，允许您通过网格坐标检索线段。例如，如果线段仅由端点表示，并且点是您的 class，您可能会得到类似这样的东西，使用每个正方形的左下角作为其定义点：
grid_seg = {
(0.5, 0.5): [p0, p1],
(1.0, 0.5): [p1, p2],
(1.0, 1.0): [p2, p3],
...
}
其中 p0、p1 等是插值交叉点。
每条轨迹都由一系列直线段组成。因此，您需要一个例程来将每条线段分成完全位于网格单元格内的部分。这种例程的基础是 Digital Differential Analyzer (DDA) 算法，但您需要修改基本算法，因为您需要每个单元格内直线的端点，而不仅仅是访问哪些单元格。
您必须注意的几件事：
1) 如果您使用的是浮点数，请注意步长值计算中的舍入错误，因为这些可能会导致算法失败。出于这个原因，许多人选择转换为整数网格，显然会损失精度。 This 是对问题的很好的讨论，有一些工作代码（虽然不是 python）。
2) 您需要决定单元格周围的 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）可以通过将点的坐标除以单元格的维度，然后舍入为整数来找到。如果当前点的单元格索引与最后一个点的单元格索引相同，则这两个点在同一单元格中。列表适合插入新点，但它的内存效率不如数组。
为轨迹创建 class 可能是个好主意。或者如果坐标列表浪费太多内存，则使用内存缓冲区和稀疏数据结构而不是列表和列表以及 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.
为了检索，制作一个包含坐标数组中起点索引的稀疏数组。
如果像元大小很大，线性插值可能看起来很糟糕。
进一步说明：
网格描述
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。delta_cell_ids[i] = cell_id_1d[i + 1] - cell_id[i]，当且仅当第 ith 和 (i + 1)th 时为零 个点在同一个单元格中。
这可以通过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);
结果：
不要用得体，自己动手：
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);