pypoints2grid 0.2.2

implemention of points2grid algorithm

pypoints2grid

Python library wrapping the points2grid algorithm for generate Digital Elevation Models (DEMs) from pointclouds. The original points2grid software was developed by the San Diego Supercomputer Center.

Acknowledgements

• Original C++ implementation
• Information about points2grid algorithm

Installation

pip install pypoints2grid

Usage

Note unlike the original C++ version of points2grid, pypoints2grid isn't a command line tool and 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 pypoints2grid import points2grid

dem = points2grid(pts, cell_size, bounds=None, radius=0, window_size=3, grid_data=['idw'], verbose=False)

Parameters

• 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 same units as points coordinates.

• 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

Returns

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.

Example

A complete script for generating a DEM from a point cloud using laspy and rasterio:

import laspy
import rasterio
import numpy as np
from time import time

from 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 bounds: ({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)

License

pypoins2grid is licensed under the MIT license. The original points2grid software is licensed under the BSD 4-clause license.