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 :

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.

it extracts all points because min, max are the bounding box of the polygon. — zwnk
– zwnk, Commented Feb 12, 2020 at 15:09
Hi, did you ever get a solution for this problem? @User18981898198119 — Fletcher_john
– Fletcher_john, Commented Aug 8, 2022 at 12:02

Marc · Accepted Answer · 2022-08-12 09:36:29Z

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

Bera · Accepted Answer · 2025-12-08 06:11:52Z

You can clip with pdals filters.crop:

The crop filter removes points that fall outside or inside a cropping bounding box (2D or 3D), polygon, or point+distance.

I read the shapefile and convert the geometry to WKT with Geopandas

import geopandas as gpd
import pdal
import json

las = r"J:/data/Laserdata_Skog/orginaldata/2020/20A012/20A012_62750_3825_25.laz"
shp = r"C:/Users/bera/Desktop/gistest/cliplas/zone.shp"

#Read the shapefile as a Geopandas dataframe
df_shape = gpd.read_file(shp)
#Extract the first polygon as as WKT geometry
polygon1 = df_shape.geometry.to_wkt().iloc[0]

#Create a pdal pipeline from a dictionary
d = {"pipeline": [las,
                  {"type":"filters.crop",
                   "polygon":polygon1}]}
json_str = json.dumps(d) #dict to json string

pipeline = pdal.Pipeline(json_str)

#Execute the pipeline and read the clipped LAS as an array
count = pipeline.execute()
array = pipeline.arrays[0]
metadata = pipeline.metadata
print(f"Point count: {count:,}")
#Point count: 2,303

#Create a point dataframe using the arrays X and Y coordinates
df_las = gpd.GeoDataFrame(
    geometry=gpd.points_from_xy(x=array["X"], y=array["Y"], crs=df_shape.crs))

df_las["height"] = array["Z"]

import matplotlib.pyplot as plt
import contextily as cx
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(20, 5), width_ratios=(2, 4))

#Plot the shapefile with a satellite image basemap
df_shape.plot(ax=ax[0], color="none", edgecolor="red", linewidth=3)
ax[0].set_title("Shapefile geometry")
cx.add_basemap(ax=ax[0], source=cx.providers["Esri"]["WorldImagery"], crs=3006)

#Plot the clipped las
ax[1].scatter(df_las.geometry.x, df_las.height, c=df_las.height)
ax[1].set_ylabel("Height")
ax[1].set_xlabel("X coordinate")
ax[1].set_title("Clipped LAS")

Stack Exchange Network

Extract point cloud data from shapefile with Python

2 Answers 2

Your Answer

Hot Network Questions

Extract point cloud data from shapefile with Python

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Related

Hot Network Questions