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

Question

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

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

经过大量的搜索和阅读，我目前的想法是使用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

我不知道我完全理解你在找什么，但因为每一个“细胞”似乎有四面八方1m的一面，不能你：

• 轮所有值整数（光栅化您的数据）可能与一些floor函数;如果您想更多easely专注于高度后

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

  你可以把z，而开头：

• 创建这些整数双射坐标的东西更方便的东西，如

  计算每个“细胞”的直方图（numpy/scipy或numpy的几个函数可以做到这一点）;

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

也许是我没有很好地理解，但在情况下，它可以帮助...

Answer 2

谢谢@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补丁提取看起来像另一种可能的解决方案。