Python library wrapping the points2grid algorithm for generate Digital Elevation Models (DEMs) from pointclouds.
Note unlike the original C++ version of points2grid, pypoints2grid doesn't read or write data from/to file or handle point cloud classification filtering. Use other libraries, like laspy and rasterio for handeling IO and preparing your data.
from src.pypoints2grid import points2grid
dem = points2grid(pts, cell_size, bounds=None, radius=0, window_size=3, grid_data=['idw'], verbose=False)-
pts: list of lists or numpy array of shape (n, 3) containing x, y, z coordinates of points
-
cell_size: size of the grid cells in meters
-
bounds: list of 4 floats containing the bounds of the grid in the form (xmin, ymin, xmax, ymax). points outside these bounds will be ignored. If None, all points will be included.
-
radius: radius of the search circle in meters. If 0, the radius will be computed from the cell size.
-
window_size: size of the window used for the moving average filter. If 0, no moving average filter will be applied. For more information about radius and window_size see the original points2grid documentation
-
grid_data: list of strings containing the data interpolations you want returned. Possible values are
- 'idw' (inverse distance weighted interpolation)
- 'min' (minimum of the points in the cell)
- 'max' (maximum of the points in the cell)
- 'mean' (mean of the points in the cell)
- 'std' (standard deviation of the points in the cell)
-
verbose: if True, print progress information
An (n,m) or (n, m, k) dimensional numpy array containing the interpolated data. n and m are the number of cells in the x and y directions, and k is the number of data interpolations requested. The order of the data interpolations is the same as the order in the grid_data parameter. If k = 1 then the returned array will be 2 dimensional.
import laspy
import rasterio
from rasterio.transform import Affine
import numpy as np
from time import time
from src.pypoints2grid import points2grid
las = laspy.read("pointcloud.las")
crs = las.header.parse_crs()
# filter out only ground (classififation 2) and water (classification 9)
ground = las.points[(las.classification == 2) | (las.classification == 9)]
pts = np.vstack((ground.x, ground.y, ground.z)).T
x_min, y_min, z_min = pts.min(axis=0)
x_max, y_max, z_max = pts.max(axis=0)
print(f"loaded {pts.shape[0]} points")
print(f"with bounnds: ({x_min}, {y_min}), ({x_max}, {y_max})")
cell_size = 0.5
print("creating grid")
start_time = time()
dem = points2grid(pts, cell_size)
print(f"grid created in {round(time() - start_time, 2)} seconds")
transform = rasterio.transform.from_origin(x_min, y_max, cell_size, cell_size)
with rasterio.open(
"dem.tif",
"w",
driver="GTiff",
height=dem.shape[0],
width=dem.shape[1],
count=1,
dtype=dem.dtype,
crs=crs,
transform=transform,
) as dst:
dst.write(dem, 1)