Rasterize a point cloud using different approaches. This stage does not modify the point cloud. It produces a derived product in raster format.
Usage
rasterize(res, operators = "max", filter = "", ofile = temptif(), ...)
Arguments
- res
numeric. The resolution of the raster. Can be a vector with two resolutions. In this case it does not correspond to the x and y resolution but to a buffered rasterization. (see section 'Buffered' and examples)
- operators
Can be a character vector. "min", "max" and "count" are accepted as well as many others (see section 'Operators'). Can also rasterize a triangulation if the input is a LASRalgorithm for triangulation (see examples). Can also be a user-defined expression (see example and section 'Operators').
- filter
the 'filter' argument allows filtering of the point-cloud to work with points of interest. For a given stage when a filter is applied, only the points that meet the criteria are processed. The most common strings are
Classification == 2"
,"Z > 2"
,"Intensity < 100"
. For more details see filters.- ofile
character. Full outputs are always stored on disk. If
ofile = ""
then the stage will not store the result on disk and will return nothing. It will however hold partial output results temporarily in memory. This is useful for stage that are only intermediate stage.- ...
default_value
numeric. When rasterizing with an operator and a filter (e.g.-keep_z_above 2
) some pixels that are covered by points may no longer contain any point that pass the filter criteria and are assigned NA. To differentiate NAs from non covered pixels and NAs from covered pixels without point that pass the filter, the later case can be assigned another value such as 0.
Value
This stage produces a raster. The path provided to `ofile` is expected to be `.tif` or any other format supported by GDAL.
Operators
If operators
is a string or a vector of strings: read metric_engine to see the possible strings
Below are some examples of valid calls:
rasterize(10, c("max", "count", "i_mean", "z_p95"))
rasterize(10, c("z_max", "c_count", "intensity_mean", "p95"))
If operators
is an R user-defined expression, the function should return either a vector of numbers
or a list
containing atomic numbers. To assign a band name to the raster, the vector or the list
must be named accordingly. The following are valid operators:
Buffered
If the argument res
is a vector with two numbers, the first number represents the resolution of
the output raster, and the second number represents the size of the windows used to compute the
metrics. This approach is called Buffered Area Based Approach (BABA).
In classical rasterization, the metrics are computed independently for each pixel. For example,
predicting a resource typically involves computing metrics with a 400 square meter pixel, resulting
in a raster with a resolution of 20 meters. It is not possible to achieve a finer granularity with
this method.
However, with buffered rasterization, it is possible to compute the raster at a resolution of 10
meters (i.e., computing metrics every 10 meters) while using 20 x 20 windows for metric computation.
In this case, the windows overlap, essentially creating a moving window effect.
This option does not apply when rasterizing a triangulation, and the second value is not considered
in this case.
Examples
f <- system.file("extdata", "Topography.las", package="lasR")
read <- reader()
tri <- triangulate(filter = keep_ground())
dtm <- rasterize(1, tri) # input is a triangulation stage
avgi <- rasterize(10, mean(Intensity)) # input is a user expression
chm <- rasterize(2, "max") # input is a character vector
pipeline <- read + tri + dtm + avgi + chm
ans <- exec(pipeline, on = f)
ans[[1]]
#> class : SpatRaster
#> dimensions : 286, 286, 1 (nrow, ncol, nlyr)
#> resolution : 1, 1 (x, y)
#> extent : 273357, 273643, 5274357, 5274643 (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83(CSRS) / MTM zone 7 (EPSG:2949)
#> source : file1f165560a8fc.tif
#> name : file1f165560a8fc
ans[[2]]
#> class : SpatRaster
#> dimensions : 30, 30, 1 (nrow, ncol, nlyr)
#> resolution : 10, 10 (x, y)
#> extent : 273350, 273650, 5274350, 5274650 (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83(CSRS) / MTM zone 7 (EPSG:2949)
#> source : file1f167181c618.tif
#> name : file1f167181c618
ans[[3]]
#> class : SpatRaster
#> dimensions : 144, 144, 1 (nrow, ncol, nlyr)
#> resolution : 2, 2 (x, y)
#> extent : 273356, 273644, 5274356, 5274644 (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83(CSRS) / MTM zone 7 (EPSG:2949)
#> source : file1f167d772f13.tif
#> name : max
# Demonstration of buffered rasterization
# A good resolution for computing point density is 5 meters.
c0 <- rasterize(5, "count")
# Computing point density at too fine a resolution doesn't make sense since there is
# either zero or one point per pixel. Therefore, producing a point density raster with
# a 2 m resolution is not feasible with classical rasterization.
c1 <- rasterize(2, "count")
# Using a buffered approach, we can produce a raster with a 2-meter resolution where
# the metrics for each pixel are computed using a 5-meter window.
c2 <- rasterize(c(2,5), "count")
pipeline = read + c0 + c1 + c2
res <- exec(pipeline, on = f)
terra::plot(res[[1]]/25) # divide by 25 to get the density
terra::plot(res[[2]]/4) # divide by 4 to get the density
terra::plot(res[[3]]/25) # divide by 25 to get the density