很抱歉。我试着缩短它，但我认为所有的东西都会对那些愿意阅读它的人有所帮助。
我有xyz点云，我正试图将它们网格化到dem(数字高程模型，如果你不熟悉-只是二维高程数组)。我的DEM需要比我的点云分辨率低得多(更确切地说，它们不需要像我的点云那样高的分辨率)，所以我的点云中每个DEM单元大约有10个点。我的代码遍历最终DEM的行和列，并根据提供的XYZ坐标计算每个单元格的高程值。为此，我指定网格大小和间距，然后代码计算每个网格单元的x和y最大值和最小值。然后，它找到x和y在该单元格的最大值和最小值范围内的所有z值，拒绝异常值，并取其余值的平均值，以确定该单元格的最终值。Xyz值存储在如下所示的数组中:
xyz = np.array([[x1, y1, z1],
[x2, y2, z2],
[x3, y3, z3]]) # with tens of thousands of xyz cominations
我目前的方法包括在X和Y中创建一个单元格边界值列表，然后将边界内具有X和Y值的所有z值添加到海拔列表:
dx = 0.5 # cell size in meters
xmin = np.min(xyz[:,0])
xmax = np.max(xyz[:,0])
ymin = np.min(xyz[:,1])
ymax = np.max(xyz[:,1])
xoffset = ((xmax-xmin) % dx)/2.0
yoffset = ((ymax-ymin) % dx)/2.0
xlims = np.arange(xmin+xoffset, xmax, dx) # list of grid cell limits in x
ylims = np.arange(ymin+yoffset, ymax, dx) # list of grid cell limits in y
DEM = np.empty((len(ylims)-1, len(xlims)-1), 'float') # declares output array
for i in range(DEM.shape[0]): # iterate over rows of final DEM
for j in range(DEM.shape[1]): # iterate over columns of final DEM
bottom = ylims[i]
top = ylims[i+1]
left = xlims[j] # these rows just pick minimum and
right = xlims[j+1] # maximum of cell [i,j]
elevations = xyz[np.where(((xyz[:,0] > left) &
(xyz[:,0] < right)) &
((xyz[:,1] > bottom) &
(xyz[:,1] < top)))[0]][:,2]
elevations = reject_outliers(elevations)
if len(elevations) == 0:
elevation = np.nan
else:
elevation = np.mean(elevations)
DEM[i,j] = elevation
这是有效的，但我必须做数百个dem，每个dem都有数十万个点，所以如果我这样做的话，我要等一个星期才能让我的电脑供电。我觉得它也很笨重。有什么办法可以简化吗?
一种可能性是找到这些点适合一次的段，然后循环遍历这些组，这样你就不必一直掩盖相同的元素，即使它们只属于一个段。
一种方法是使用numpy内置的searchsorted:xmin = np.min(xyz[:,0])
xmax = np.max(xyz[:,0])
ymin = np.min(xyz[:,1])
ymax = np.max(xyz[:,1])
xoffset = ((xmax-xmin) % dx)/2.0
yoffset = ((ymax-ymin) % dx)/2.0
xlims = np.arange(xmin+xoffset, xmax, dx) # list of grid cell limits in x
ylims = np.arange(ymin+yoffset, ymax, dx) # list of grid cell limits in y
DEM = np.empty((len(ylims) - 1, len(xlims) - 1), dtype=float) # declares output array
# Find the bins that each point fit into
x_bins = np.searchsorted(xlims, xyz[:, 0]) - 1
y_bins = np.searchsorted(ylims, xyz[:, 1]) - 1
for i in range(DEM.shape[0]): # iterate over rows of final DEM
y_mask = y_bins == i
for j in range(DEM.shape[1]):
elevations = xyz[y_mask & (x_bins == j), 2]
elevations = reject_outliers(elevations)
if len(elevations) == 0:
elevations = np.nan
else:
elevations = np.mean(elevations)
DEM[i, j] = elevations
当我使用timeit与xyz = randn(100000, 3)和dx = 0.1分析已经列出的替代方案时(在定义reject_outliers为lambda x: x之后)，我得到了以下时间:
- (dermon's)逐步法:10次循环，最佳3:16.7 s/循环
- 搜索排序方法:10个循环，最佳3:2.11 s每循环
然而，如果你愿意使用Pandas，你可以修改代码，使双for循环被Pandas出色的groupby功能所取代:
elevation_df = pd.DataFrame({'x_bins': x_bins, 'y_bins': y_bins, 'z': xyz[:, 2]})
for x_y_bins, data in elevation_df.groupby(['x_bins', 'y_bins']):
elevations = reject_outliers(data['z'])
elevations = data['z']
if len(elevations) == 0:
elevations = np.nan
else:
elevations = np.mean(elevations)
if 0 <= x_y_bins[1] < DEM.shape[0] and 0 <= x_y_bins[0] < DEM.shape[1]:
DEM[x_y_bins[1], x_y_bins[0]] = elevations
这几乎将时间降低了2倍(10个循环，最好为3:1.11 s/循环)。
我还应该注意到，由于你的np，你似乎已经排除了范围中的一些点。不等命令。在上面我假设你的意思是排除这些点，但如果你想包括所有的数据，你可以使用:
xlims = np.arange(xmin-xoffset, xmax+dx, dx)
ylims = np.arange(ymin-yoffset, ymax+dx, dx)
如果你使用这些范围，你可以修改我之前的for循环:
for i in np.unique(y_bins):
y_mask = y_bins == i
for j in np.unique(x_bins[y_mask]):
将我之前的搜索排序示例的结果时间降低到10个循环，最好的是3:1.56 s每个循环，这至少更接近Pandas组
你可以试试
elev = xyz[ xyz[:,0] > left ]
elev = elev[ elev[:,0] < right ]
elev = elev[ elev[:,1] > bottom]
elev = elev[ elev[:,1] < top ]
这样，每个条件语句都减少了elev的大小，因此每个下一个条件语句需要考虑的对象更少。不过，我怀疑还有更好的方法。您也可以尝试从xyz搜索最近的x,y到xlims, ylims的x,y，然后使用dx, dy返回相邻坐标。如果你能找到一种快速的搜索方式，那么这可能是你最好的选择(但是当你接近边界时，你需要小心
你可以尝试scipy.interpolate.RectBivariateSpline，如果x,y在xyz形成一个规则的细网格。
否则可以使用scipy.interpolate.interp2d:
from scipy.interpolate import interp2d
interped = interp2d( x=xyz[:,0], y=xyz[:,1], z=xyz[:,2], kind='linear', fill_value=xyz[:,2].mean() )
DEM = interped( xlims, ylims)
(它将比RectBivariateSpline慢得多，但可能比你正在做的快得多)
您可以将'kind'参数更改为'cubic'，但它会减慢速度。