Extract point cloud data from shapefile with Python

Question

I am trying to plot profiles from LAS data. So my first step is to extract the data using a shapefile. I found this piece of code on the laspy website and thought I could use it to do this :

def explode(coords):
 for e in coords:
 if isinstance(e, (int, float,complex)):
 yield coords
 break
 else:
 for f in explode(e):
 yield f
def bbox(f):
 x, y = zip(*list(explode(f['geometry']['coordinates'])))
 return min(x), min(y), max(x), max(y)
def clip_extract(las_file , shp_file ):
 import laspy, os
 import numpy as np
 inFile = laspy.file.File(las_file, mode = "r")
 shp = fiona.open(shp_file, "r")
 min_x, min_y, max_x, max_y = acav.bbox(shp)
 X_invalid = np.logical_and((min_x <= inFile.x),
 (max_x >= inFile.x))
 Y_invalid = np.logical_and((min_y <= inFile.y),
 (max_y >= inFile.y))
 Z_invalid = np.logical_and((inFile.header.min[2] <= inFile.z),
 (inFile.header.max[2] >= inFile.z))
 good_indices = np.where(np.logical_and(X_invalid, Y_invalid, Z_invalid))
 good_points = inFile.points[good_indices]

But of course, it extracts the full min(x), min(y), max(x), max(y) rectangle while ideally it would only extract the points within the shapefile : Points selected vs shapefile

I have read many other posts on how to mask a numpy array with a shapefile but I can't seem to find one that matches my case... I am trying to do this with Fiona or shapely, and Numpy.

• zwnk: it extracts all points because min, max are the bounding box of the polygon.
• Fletcher_john: Hi, did you ever get a solution for this problem? @User18981898198119

Answer 1

you should look at pdal.

https://pdal.io/en/latest/stages/filters.crop.html

import pdal
polygon_draw_bounds = polygon_draw.bounds
cropper = {
 "pipeline": [str(datasource/file),
 { "type":"filters.crop",
 'bounds':str(([polygon_draw_bounds[0], polygon_draw_bounds[2]],[polygon_draw_bounds[1], polygon_draw_bounds[3]]))},
 { "type":"filters.crop",
 'polygon':polygon_draw.wkt
 },
 { "type":"writers.las",
 "filename": str(output)+"\\"+"cloud_crop.laz"
 },
 ]}
pipeline = pdal.Pipeline(json.dumps(cropper))
pipeline.execute()