lidR (version 4.1.2)

rasterize: Rasterize a point cloud

Description

Rasterize a point cloud in different ways to compute a DTM, a CHM or a density map. Most raster products can be computed with pixel_metrics but some are more complex and require dedicated and optimized functions. See Details and Examples.

Usage

rasterize_canopy(las, res = 1, algorithm = p2r(), ...)
rasterize_density(las, res = 4, ...)
rasterize_terrain(
  las,
  res = 1,
  algorithm = tin(),
  use_class = c(2L, 9L),
  shape = "convex",
  ...
)

Value

RasterLayer or a stars or a SpatRaster depending on the settings.

Arguments

las

An object of class LAS or LAScatalog.

res

numeric. The size of a grid cell in point cloud coordinates units. Can also be RasterLayer or a stars or a SpatRaster used as layout.

algorithm

function. A function that implements an algorithm to compute a digital surface model or a digital terrain model. lidR implements p2r, dsmtin, pitfree for digital surface models, and knnidw, tin, and kriging for digital terrain models (see respective documentation and examples).

...

Use pkg = "terra|raster|stars" to get an output in SpatRaster, RasterLayer or stars format. Default is getOption("lidR.raster.default").

use_class

integer vector. By default the terrain is computed by using ground points (class 2) and water points (class 9).

shape

By default the interpolation is made only within the "convex" hull of the point cloud to get a DTM with the shape of the point cloud. This prevents meaningless interpolations where there is no data. It can also be "concave" or "bbox". It can also be an sfc to define a polygon in which to perform the interpolation.

Non-supported LAScatalog options

The option select is not supported and not respected in rasterize_* because it is internally known what is best to select.
The option chunk_buffer is not supported and not respected in rasterize_canopy and rasterize_density because it is not necessary.

Details

rasterize_terrain

Interpolates the ground points and creates a rasterized digital terrain model. The algorithm uses the points classified as "ground" and "water" (Classification = 2 and 9, respectively, according to LAS file format specifications) to compute the interpolation. How well the edges of the dataset are interpolated depends on the interpolation method used. A buffer around the region of interest is always recommended to avoid edge effects.

rasterize_canopy

Creates a digital surface model (DSM) using several possible algorithms. If the user provides a normalized point cloud, the output is indeed a canopy height model (CHM).

rasterize_density

Creates a map of the point density. If a "pulseID" attribute is found, also returns a map of the pulse density.

Examples

# =====================
# Digital Terrain Model
# =====================

LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las = readLAS(LASfile, filter = "-inside 273450 5274350 273550 5274450")
#plot(las)

dtm1 = rasterize_terrain(las, algorithm = knnidw(k = 6L, p = 2))
dtm2 = rasterize_terrain(las, algorithm = tin())

if (FALSE) {
dtm3 = rasterize_terrain(las, algorithm = kriging(k = 10L))

plot(dtm1, col = gray(0:25/25))
plot(dtm2, col = gray(0:25/25))
plot(dtm3, col = gray(0:25/25))
plot_dtm3d(dtm1)
plot_dtm3d(dtm2)
plot_dtm3d(dtm3)
}

# =====================
# Digital Surface Model
# =====================

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile, filter = "-inside 481280 3812940 481330 3812990")
col <- height.colors(15)

# Points-to-raster algorithm with a resolution of 1 meter
chm <- rasterize_canopy(las, res = 1, p2r())
plot(chm, col = col)

# Points-to-raster algorithm with a resolution of 0.5 meters replacing each
# point by a 20-cm radius circle of 8 points
chm <- rasterize_canopy(las, res = 0.5, p2r(0.2))
plot(chm, col = col)

# Basic triangulation and rasterization of first returns
chm <- rasterize_canopy(las, res = 0.5, dsmtin())
plot(chm, col = col)

# Khosravipour et al. pitfree algorithm
chm <- rasterize_canopy(las, res = 0.5, pitfree(c(0,2,5,10,15), c(0, 1.5)))
plot(chm, col = col)

# ====================
# Digital Density Map
# ====================

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile,  filter = "-inside 684800 5017800 684900 5017900")

d <- rasterize_density(las, 5)
plot(d)

las <- retrieve_pulses(las)
d <- rasterize_density(las)
plot(d)