Regions of Interest

Relevant Resources

Overview

We demonstrate the selection of regions of interest (ROIs) from LiDAR data. Geometries like circles and rectangles are selected based on coordinates. Complex geometries are extracted from shapefiles to clip specific areas.

Environment

# Clear environment
rm(list = ls(globalenv()))

# Load packages
library(lidR)
library(sf)

Simple Geometries

Load LiDAR Data and Inspect

We start by loading some LiDAR data and inspecting its header and number of point records.

las <- readLAS(files = 'data/zrh_norm.laz')

# Inspect the header and the number of point records
las@header
#> File signature:           LASF 
#> File source ID:           0 
#> Global encoding:
#>  - GPS Time Type: GPS Week Time 
#>  - Synthetic Return Numbers: no 
#>  - Well Know Text: CRS is GeoTIFF 
#>  - Aggregate Model: false 
#> Project ID - GUID:        00000000-0000-0000-5053-004d00000000 
#> Version:                  1.1
#> System identifier:         
#> Generating software:      rlas R package 
#> File creation d/y:        0/0
#> header size:              227 
#> Offset to point data:     305 
#> Num. var. length record:  1 
#> Point data format:        1 
#> Point data record length: 28 
#> Num. of point records:    1270080 
#> Num. of points by return: 471651 342346 241371 135987 56797 
#> Scale factor X Y Z:       0.01 0.01 0.01 
#> Offset X Y Z:             0 0 0 
#> min X Y Z:                2670505 1258734 -0.88 
#> max X Y Z:                2670755 1258984 48.31 
#> Variable Length Records (VLR):
#>    Variable Length Record 1 of 1 
#>        Description: by LAStools of rapidlasso GmbH 
#>        Tags:
#>           Key 1024 value 1 
#>           Key 3072 value 2056 
#> Extended Variable Length Records (EVLR):  void
las@header$`Number of point records`
#> [1] 1270080

Select Circular and Rectangular Areas

We can select circular and rectangular areas from the LiDAR data based on specified coordinates and radii or dimensions.

# Establish coordinates
x <- 2670592
y <- 1258890

# Select a circular area
circle <- clip_circle(las = las, xcenter = x, ycenter = y, radius = 30)

# Inspect the circular area and the number of point records
circle
#> class        : LAS (v1.1 format 1)
#> memory       : 3.1 Mb 
#> extent       : 2670562, 2670622, 1258860, 1258920 (xmin, xmax, ymin, ymax)
#> coord. ref.  : CH1903+ / LV95 
#> area         : 2903 m²
#> points       : 58.4 thousand points
#> type         : airborne
#> density      : 20.13 points/m²
#> density      : 7.55 pulses/m²
circle@header$`Number of point records`
#> [1] 58447

# Plot the circular area
plot(circle)

We can do the same with a rectangular area by defining corner coordinates.

# Select a rectangular area
rect <- clip_rectangle(las = las, xleft = x, ybottom = y, xright = x + 40, ytop = y + 30)

# Plot the rectangular area
plot(rect)

We can also supply multiple coordinate pairs to clip multiple ROIs.

# Select multiple random circular areas
x <- runif(2, x, x)
y <- runif(2, 1258840, 1258890)

plots <- clip_circle(las = las, xcenter = x, ycenter = y, radius = 10)

# Plot each of the multiple circular areas
plot(plots[[1]])

# Plot each of the multiple circular areas
plot(plots[[2]])

Extraction of Complex Geometries from Shapefiles

We demonstrate how to extract complex geometries from shapefiles using the clip_roi() function from the lidR package.

We use the sf package to load an ROI and then clip to its extents.

# Load the shapefile using sf
stand_bdy <- sf::st_read(dsn = "data/roi/roi.gpkg", quiet = TRUE)

# Plot the LiDAR header information without the map
plot(las@header, map = FALSE)

# Plot the stand boundary areas on top of the LiDAR header plot
plot(stand_bdy, add = TRUE, col = "#08B5FF39")


# Extract points within the stand boundary using clip_roi()
stand <- clip_roi(las = las, geometry = stand_bdy)

# Plot the extracted points within the planting areas
plot(stand)

Clipping ROIs with a catalog

We clip the LAS data in the catalog using specified coordinate groups.


ctg <- catalog("data/ctg_norm")

# Set coordinate groups
x <- c(2670578, 2671234, 2671499, 2671755, 2671122)
y <- c(1258601, 1259050, 1259450, 1259900, 1258750)

# Visualize coordinate groups
plot(ctg)
points(x, y)


# Clip plots
rois <- clip_circle(las = ctg, xcenter = x, ycenter = y, radius = 30)

plot(rois[[1]])

plot(rois[[3]])

Validate clipped data

We validate the clipped LAS data using the las_check function.

las_check(rois[[1]])
#> 
#>  Checking the data
#>   - Checking coordinates... ✓
#>   - Checking coordinates type... ✓
#>   - Checking coordinates range... ✓
#>   - Checking coordinates quantization... ✓
#>   - Checking attributes type... ✓
#>   - Checking ReturnNumber validity... ✓
#>   - Checking NumberOfReturns validity... ✓
#>   - Checking ReturnNumber vs. NumberOfReturns... ✓
#>   - Checking RGB validity... ✓
#>   - Checking absence of NAs... ✓
#>   - Checking duplicated points... ✓
#>   - Checking degenerated ground points... ✓
#>   - Checking attribute population...
#>     🛈 'ScanDirectionFlag' attribute is not populated
#>     🛈 'EdgeOfFlightline' attribute is not populated
#>   - Checking gpstime incoherances ✓
#>   - Checking flag attributes... ✓
#>   - Checking user data attribute... ✓
#>  Checking the header
#>   - Checking header completeness... ✓
#>   - Checking scale factor validity... ✓
#>   - Checking point data format ID validity... ✓
#>   - Checking extra bytes attributes validity... ✓
#>   - Checking the bounding box validity... ✓
#>   - Checking coordinate reference system... ✓
#>  Checking header vs data adequacy
#>   - Checking attributes vs. point format... ✓
#>   - Checking header bbox vs. actual content... ✓
#>   - Checking header number of points vs. actual content... ✓
#>   - Checking header return number vs. actual content... ✓
#>  Checking coordinate reference system...
#>   - Checking if the CRS was understood by R... ✓
#>  Checking preprocessing already done 
#>   - Checking ground classification... yes
#>   - Checking normalization... yes
#>   - Checking negative outliers...
#>     ⚠ 6 points below 0
#>   - Checking flightline classification... yes
#>  Checking compression
#>   - Checking attribute compression...
#>    -  ScanDirectionFlag is compressed
#>    -  EdgeOfFlightline is compressed
#>    -  Synthetic_flag is compressed
#>    -  Keypoint_flag is compressed
#>    -  Withheld_flag is compressed
#>    -  UserData is compressed
las_check(rois[[3]])
#> 
#>  Checking the data
#>   - Checking coordinates... ✓
#>   - Checking coordinates type... ✓
#>   - Checking coordinates range... ✓
#>   - Checking coordinates quantization... ✓
#>   - Checking attributes type... ✓
#>   - Checking ReturnNumber validity... ✓
#>   - Checking NumberOfReturns validity... ✓
#>   - Checking ReturnNumber vs. NumberOfReturns... ✓
#>   - Checking RGB validity... ✓
#>   - Checking absence of NAs... ✓
#>   - Checking duplicated points... ✓
#>   - Checking degenerated ground points... ✓
#>   - Checking attribute population...
#>     🛈 'ScanDirectionFlag' attribute is not populated
#>     🛈 'EdgeOfFlightline' attribute is not populated
#>   - Checking gpstime incoherances
#>     ✗ 5 pulses (points with the same gpstime) have points with identical ReturnNumber
#>   - Checking flag attributes... ✓
#>   - Checking user data attribute... ✓
#>  Checking the header
#>   - Checking header completeness... ✓
#>   - Checking scale factor validity... ✓
#>   - Checking point data format ID validity... ✓
#>   - Checking extra bytes attributes validity... ✓
#>   - Checking the bounding box validity... ✓
#>   - Checking coordinate reference system... ✓
#>  Checking header vs data adequacy
#>   - Checking attributes vs. point format... ✓
#>   - Checking header bbox vs. actual content... ✓
#>   - Checking header number of points vs. actual content... ✓
#>   - Checking header return number vs. actual content... ✓
#>  Checking coordinate reference system...
#>   - Checking if the CRS was understood by R... ✓
#>  Checking preprocessing already done 
#>   - Checking ground classification... yes
#>   - Checking normalization... yes
#>   - Checking negative outliers...
#>     ⚠ 5 points below 0
#>   - Checking flightline classification... yes
#>  Checking compression
#>   - Checking attribute compression...
#>    -  ScanDirectionFlag is compressed
#>    -  EdgeOfFlightline is compressed
#>    -  Synthetic_flag is compressed
#>    -  Keypoint_flag is compressed
#>    -  Withheld_flag is compressed
#>    -  UserData is compressed

Independent files (e.g. plots) as catalogs

We read an individual LAS file as a catalog and perform operations on it.

# Read single file as catalog
ctg <- readLAScatalog(folder = "data/zrh_norm.laz")

# Set options for output files
opt_output_files(ctg) <- paste0(tempdir(),"/{XCENTER}_{XCENTER}")

# Write file as .laz
opt_laz_compression(ctg) <- TRUE

# Get random plot locations and clip
x <- runif(n = 4, min = ctg$Min.X, max = ctg$Max.X)
y <- runif(n = 4, min = ctg$Min.Y, max = ctg$Max.Y)
rois <- clip_circle(las = ctg, xcenter = x, ycenter = y, radius = 10)

#> Chunk 1 of 4 (25%): state ✓
#> 
                                                                                
Chunk 2 of 4 (50%): state ✓
#> 
                                                                                
Chunk 3 of 4 (75%): state ✓
#> 
                                                                                
Chunk 4 of 4 (100%): state ✓
#>

# Read catalog of plots
ctg_plots <- readLAScatalog(tempdir())

# Set independent files option
opt_independent_files(ctg_plots) <- TRUE
opt_output_files(ctg_plots) <- paste0(tempdir(),"/{XCENTER}_{XCENTER}")

# Generate plot-level terrain models
rasterize_terrain(las = ctg_plots, res = 1, algorithm = tin())

#> Chunk 1 of 4 (25%): state ✓
#> 
                                                                                
Chunk 2 of 4 (50%): state ✓
#> 
                                                                                
Chunk 3 of 4 (75%): state ✓
#> 
                                                                                
Chunk 4 of 4 (100%): state ✓
#> 
                                                                                
#> class       : SpatRaster 
#> dimensions  : 92, 202, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 2670512, 2670714, 1258842, 1258934  (xmin, xmax, ymin, ymax)
#> coord. ref. : CH1903+ / LV95 (EPSG:2056) 
#> source      : rasterize_terrain.vrt 
#> name        : Z 
#> min value   : 0 
#> max value   : 0

# Check files
path <- paste0(tempdir())
file_list <- list.files(path, full.names = TRUE)
file <- file_list[grep("\\.tif$", file_list)][[1]]

# plot dtm
plot(terra::rast(file))

Conclusion

This concludes our tutorial on selecting simple geometries and extracting complex geometries from geopackage files (or shapefiles etc…) using the lidR package in R.