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

Question

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

从概念上讲，该操作与叠加光栅和计算与每个像素相交的点相同。 细胞网格可以形成正方形或圆盘，这无关紧要。在

经过大量的搜索和阅读，我现在的想法是用一些numpy.linspace以及numpy.meshgrid生成存储在数组中的每个单元格的顶点，并针对每个点测试每个单元格，以查看它是否为“in”。这似乎效率低下，尤其是在处理数千个点时。在

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

我不确定我完全理解你在找什么，但既然每个“细胞”似乎都有一个1米的侧面，可以向各个方向移动，你能不能：

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

• 从这些整数坐标创建一个双投影到更方便的位置，例如：

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

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

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

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

也许我不太明白，但万一有帮助。。。在

Answer 2

谢谢@barucher。原来@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补丁提取看起来是另一个可能的解决方案。在