2

我arr用模块创建了一个 NumPy 结构化数组（称为 ）arcpy：

arr = arcpy.da.FeatureClassToNumPyArray('MPtest','SHAPE@XYZ',explode_to_points=True)

数组看起来像（仅显示第一行）：

array([([309243.1420999998, 6143470.0293, 470.43829999999434],),
       ([309252.14080000017, 6143461.0857, 470.43829999999434],),
       ([309246.0201000003, 6143459.2754, 470.43829999999434],),
       ........................................................,
       ([309252.14080000017, 6143461.0857, 464.6000000000058],)], 
      dtype=[('SHAPE@XYZ', '<f8', (3,))])

它表示取自'MPtest'ArcGIS 中构建的 3d 对象 ( ) 的顶点的 XYZ 坐标（multipatch几何）。

我有另一个 NumPy 数组（通过使用模块读取.las文件生成laspy），名为point_cloud. 这个数组看起来像：

[((30922713, 614349673, 46679, 154, 17, 1, -10, 79, 5,  11570.850892),)
 ((30922712, 614349659, 46674, 112, 17, 1, -10, 78, 5,  11570.850897),)
 ((30922709, 614349645, 46663, 161, 17, 1, -10, 77, 5,  11570.850902),),
 ..................................................................,)],
[('point', [('X', '<i4'), ('Y', '<i4'), ('Z', '<i4'), ('intensity', '<u2'), ('flag_byte', 'u1'), ('raw_classification', 'u1'), ('scan_angle_rank', 'i1'), ('user_data', 'u1'), ('pt_src_id', '<u2'), ('gps_time', '<f8')])]

我希望能够获得point cloud属于arr. 这甚至可能吗？

我一直在尝试使用诸如np.where、np.intersect1d、np.logical_and和 finally之类的功能np.vstack，但到目前为止，我无法做到这一点。此外，我在 Python 方面有相当扎实的背景，但 NumPy 对我来说仍然是一种新的和非常复杂的东西（至少乍一看......）。

4

1 回答 1

2

一旦你得到非结构化数组（正如我看到你已经在评论中实现的那样），你可以使用scipy.spatial.Delanuay如下：

我只是创建了一个示例框和一些要点来阐明示例：

import numpy as np
from itertools import product
from scipy.spatial import Delaunay

arr = np.array(list(product([0,1], repeat=3)))
arr
array([[0, 0, 0],
       [0, 0, 1],
       [0, 1, 0],
       [0, 1, 1],
       [1, 0, 0],
       [1, 0, 1],
       [1, 1, 0],
       [1, 1, 1]])
point_cloud = np.array([[0.5,0.5,0.5], [2,2,2]])

然后创建 Delanuay 三角剖分：

dt = Delaunay(arr)

并找到 point_cloud 的哪些点在 dt 内（arr 的 Delanuay 三角剖分）：

points_inside_arr = dt.find_simplex(point_cloud) >=0
points_inside_arr
array([ True, False], dtype=bool)

这会产生一个 numpy 布尔数组，指示点云中的哪些点在 arr 内。

于 2017-03-22T18:11:07.557 回答