使用 numpy 数组计算规则网格单元格内的点数

Question

我正在处理大量 3D 点，每个点的 x、y、z 值都存储在 numpy 数组中。 对于背景，这些点将始终落在一个固定半径的圆柱体中，并且高度 = 点的最大 z 值。 我的目标是将边界圆柱体（或列，如果它更容易）分成例如1 m高地层，然后计算每个单元格内的点数 覆盖在每个地层上的规则网格（例如 1 m x 1 m）。

从概念上讲，该操作与覆盖栅格并计算与每个像素相交的点数相同。 网格单元格可以形成正方形或圆盘，没关系。

经过大量搜索和阅读，我目前的想法是使用 numpy.linspace 和 numpy.meshgrid 的某种组合来生成存储在数组中的每个单元格的顶点，并针对每个点测试每个单元格，看看它是否是'在'。这似乎效率低下，尤其是在处理数千个点时。

numpy / scipy 套件似乎很适合这个问题，但我还没有找到解决方案。任何建议将不胜感激。 我已经包含了一些示例点和一些代码来可视化数据。

# Setup
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Load in X,Y,Z values from a sub-sample of 10 points for testing
# XY Values are scaled to a reasonable point of origin
z_vals = np.array([3.08,4.46,0.27,2.40,0.48,0.21,0.31,3.28,4.09,1.75])
x_vals = np.array([22.88,20.00,20.36,24.11,40.48,29.08,36.02,29.14,32.20,18.96])
y_vals = np.array([31.31,25.04,31.86,41.81,38.23,31.57,42.65,18.09,35.78,31.78])

# This plot is instructive to visualize the problem
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x_vals, y_vals, z_vals, c='b', marker='o')
plt.show()

Answer 1

感谢@Baruchel。事实证明，@DilithiumMatrix 建议的 n 维直方图为我发布的问题提供了一个相当简单的解决方案。经过一番阅读，这是我目前为面临类似问题的其他人提供的解决方案。 由于这是我第一次尝试 Python/Numpy，因此欢迎任何改进/建议，尤其是关于性能的改进/建议。谢谢。

# Setup
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Load in X,Y,Z values from a sub-sample of 10 points for testing
# XY Values are scaled to a reasonable point of origin
z_vals = np.array([3.08,4.46,0.27,2.40,0.48,0.21,0.31,3.28,4.09,1.75])
x_vals = np.array([22.88,20.00,20.36,24.11,40.48,29.08,36.02,29.14,32.20,18.96])
y_vals = np.array([31.31,25.04,31.86,41.81,38.23,31.57,42.65,18.09,35.78,31.78])

# Updated code below
# Variables needed for 2D,3D histograms
xmax, ymax, zmax = int(x_vals.max())+1, int(y_vals.max())+1, int(z_vals.max())+1
xmin, ymin, zmin = int(x_vals.min()), int(y_vals.min()), int(z_vals.min())
xrange, yrange, zrange = xmax-xmin, ymax-ymin, zmax-zmin
xedges = np.linspace(xmin, xmax, (xrange + 1), dtype=int)
yedges = np.linspace(ymin, ymax, (yrange + 1), dtype=int)
zedges = np.linspace(zmin, zmax, (zrange + 1), dtype=int)

# Make the 2D histogram
h2d, xedges, yedges = np.histogram2d(x_vals, y_vals, bins=(xedges, yedges))
assert np.count_nonzero(h2d) == len(x_vals), "Unclassified points in the array"
extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
plt.imshow(h2d.transpose(), extent=extent,  interpolation='none', origin='low')
# Transpose and origin must be used to make the array line up when using imshow, unsure why
# Plot settings, not sure yet on matplotlib update/override objects
plt.grid(b=True, which='both')
plt.xticks(xedges)
plt.yticks(yedges)
plt.xlabel('X-Axis')
plt.ylabel('Y-Axis')
plt.plot(x_vals, y_vals, 'ro')
plt.show()

# 3-dimensional histogram with 1 x 1 x 1 m bins. Produces point counts in each 1m3 cell.
xyzstack = np.stack([x_vals,y_vals,z_vals], axis=1)
h3d, Hedges = np.histogramdd(xyzstack, bins=(xedges, yedges, zedges))
assert np.count_nonzero(h3d) == len(x_vals), "Unclassified points in the array"
h3d.shape  # Shape of the array should be same as the edge dimensions
testzbin = np.sum(np.logical_and(z_vals >= 1, z_vals < 2))  # Slice to test with
np.sum(h3d[:,:,1]) == testzbin  # Test num points in second bins
np.sum(h3d, axis=2)  # Sum of all vertical points above each x,y 'pixel'
# only in this example the h2d and np.sum(h3d,axis=2) arrays will match as no z bins have >1 points

# Remaining issue - how to get a r x c count of empty z bins.
# i.e. for each 'pixel'  how many z bins contained no points?
# Possible solution is to reshape to use logical operators
count2d = h3d.reshape(xrange * yrange, zrange)  # Maintain dimensions per num 3D cells defined
zerobins = (count2d == 0).sum(1)
zerobins.shape
# Get back to x,y grid with counts - ready for output as image with counts=pixel digital number
bincount_pixels = zerobins.reshape(xrange,yrange)
# Appears to work, perhaps there is a way without reshapeing?

PS 如果您面临类似的问题，scikit 补丁提取似乎是另一种可能的解决方案。

Answer 2

我不确定我是否完全理解您在寻找什么，但是由于每个“单元”似乎都有一个 1m 的边用于所有方向，您可以吗：

• 可能使用floor 函数将所有值四舍五入为整数（栅格化数据）；

• 从这些整数坐标创建一个双射到更方便的东西，例如：

  (64**2)*x + (64)*y + z # assuming all values are in [0,63]

  如果您想稍后更轻松地关注高度，可以将z 放在开头

• 计算每个“单元格”的直方图（numpy/scipy 或 numpy 的几个函数都可以做到）；

• 如果需要，恢复双射（即，一旦知道计数，就知道每个单元格的“真实”坐标）

也许我不太了解，但如果它有帮助......