2

I have created a NumPy structured array (called arr) with arcpy module doing:

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

The array looks like (only first rows are shown):

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,))])

It represents XYZ coordinates taken from the vertices of a 3d object ('MPtest') built in ArcGIS (a multipatch geometry).

I have another NumPy array (generated from reading a .las file using the laspy module), called point_cloud. This array looks like:

[((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')])]

I'd like to be able to possibly get the indices of the points of point cloud which fall within arr. Is this even possible?

I've been trying to play around with functions like np.where, np.intersect1d, np.logical_and, and finally np.vstack, but so far, I wasn't able to do that. Also, I have a quite robust background in Python, but still NumPy is kind of new and very complex to my eyes (at least at first sight...).

umbe1987
  • 2,894
  • 6
  • 35
  • 63
  • 1
    Define `fall inside` clearly. What kinds of things have you tried? Using structured arrays may be making the task harder. – hpaulj Mar 22 '17 at 15:50
  • @hpaulj With `fall inside` I mean "the points of the point cloud that are inside the volume occupied by the 3d object". Also, I noticed that I can access the structured array like `arr['SHAPE@XYZ']` to obtain a "non-structured array", which might be easy to work with. Finally, I'll edit my question to include some tests I made as soon as possible. Thanks for the response. – umbe1987 Mar 22 '17 at 16:14
  • Similarly `point_cloud['point']`, though its three fields should be changed to 3 columns of a 2d array. I forget the correct use of `astype` or `view` for that task. Check `scipy` for `qhull` functions. I'm aware of them but haven't used them. – hpaulj Mar 22 '17 at 16:26
  • 1
    `pointcloud['point'][['X','Y','Z']].copy().view('3i')` to get a 2d 3 column array. – hpaulj Mar 22 '17 at 16:49
  • Thanks, the `pointcloud` becomes a 2d 3 columns array. I am now checking the [`scipy.spatial.ConvexHull`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.ConvexHull.html) function to see if I can implement it in my code. – umbe1987 Mar 22 '17 at 17:06

1 Answers1

2

Once you get the unstructured arrays (as I see that you already achieve in comments) you can use scipy.spatial.Delanuay as follows:

I just create a sample box and some points to clarify the example:

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]])

You then create the Delanuay Triangulation:

dt = Delaunay(arr)

And find wich points of point_cloud are inside dt (The Delanuay triangulation of arr):

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

This results in a numpy boolean array indicating which of the points in point cloud are inside arr.

David de la Iglesia
  • 2,436
  • 14
  • 29