create hillshade for LAScatalog

[wiesehahn]

Hey, I would appreciate help for the following problem.

I have many tiles of lidar data for which I want to create a hillshade of the terrain.
To not run into the problem of edge effects the idea is to apply a hillshade algorithm on a LAScatalog (with buffer).
Since I try to get high resolution shades for a broad area mosaicing them for the entire area of the catalog is not a good ideas and I decided to write them as tiles likes the point cloud data.

For this I tried to follow your advices on the LAScatalog processing engine. However, I am still struggling. I guess there is still some things which I misunderstood, here is my code:

library(here)
library(lidR)


create_hillshade <- function(chunk)
{
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)
  
  dtm <- rasterize_terrain(las, 0.25, tin(), pkg = "terra")
  
  dtm_prod <- terra::terrain(dtm, v = c("slope", "aspect"), unit = "radians")
  dtm_hillshade <- terra::shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
  
  rastershade <- crop(dtm_hillshade, ext(las))
  return(rastershade)  
}


hillshade <- function(las)
{
  if (is(las, "LAScatalog"))  {
    options <- list(automerge = FALSE, need_buffer = TRUE)
    catalog_apply(las, create_hillshade, .options = options)
  }
  else if (is(las, "LAScluster")) {
    las <- readLAS(las)
    if (is.empty(las)) return(NULL) 
    create_hillshade(las)
  }
  else if (is(las, "LAS")) {
    dtm <- rasterize_terrain(las, 0.25, tin(), pkg = "terra")
    dtm_prod <- terra::terrain(dtm, v = c("slope", "aspect"), unit = "radians")
    rastershade <- terra::shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
    return(rastershade)
  }
  else {
    stop("This type is not supported.")
  }
}

# set up parallel processing
library(future)
plan(multisession) 

load("lidr-catalog.RData")

temp_path <- here("data")
opt_output_files(ctg) <- paste0(temp_path, "/{ORIGINALFILENAME}_hillshade")
opt_chunk_buffer(ctg) <- 30
opt_chunk_size(ctg) <- 0 #  Set 0 use the tiling pattern (processing by file)
opt_select(ctg) <- "xyzc"
opt_filter(ctg) <- "-keep_class 2"
opt_merge(ctg) <- FALSE

summary(ctg)

# test function for selected tiles
ctg$processed <- FALSE
ctg$processed[1:3] <- TRUE

hillshade(ctg)

which results in Error: Trying to write an object of class SpatRaster but this type is not supported. (although this driver should be supported)

What is the problem(s)?

[Jean-Romain]

(although this driver should be supported)

Sure. Your issue is unexpected. Your code looks pretty good at first glance. I spotted one error but that is not related to your problem

  rastershade <- crop(dtm_hillshade, ext(chunk))

[Jean-Romain]

The following minimal reproducible example works for me

library(lidR)

create_hillshade <- function(chunk)
{
  las <- readLAS(chunk)
  if (is.empty(las)) return(NULL)
  
  dtm <- rasterize_terrain(las, 0.25, tin(), pkg = "terra")
  
  dtm_prod <- terra::terrain(dtm, v = c("slope", "aspect"), unit = "radians")
  dtm_hillshade <- terra::shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
  
  rastershade <- crop(dtm_hillshade, ext(chunk))
  return(rastershade)  
}


hillshade <- function(las)
{
  if (is(las, "LAScatalog"))  {
    options <- list(automerge = FALSE, need_buffer = TRUE)
    catalog_apply(las, create_hillshade, .options = options)
  }
  else if (is(las, "LAScluster")) {
    las <- readLAS(las)
    if (is.empty(las)) return(NULL) 
    create_hillshade(las)
  }
  else if (is(las, "LAS")) {
    dtm <- rasterize_terrain(las, 0.25, tin(), pkg = "terra")
    dtm_prod <- terra::terrain(dtm, v = c("slope", "aspect"), unit = "radians")
    rastershade <- terra::shade(slope = dtm_prod$slope, aspect = dtm_prod$aspect)
    return(rastershade)
  }
  else {
    stop("This type is not supported.")
  }
}

# set up parallel processing
library(future)
plan(multisession) 

LASfile <- system.file("extdata", "Topography.laz", package="lidR")
ctg = readLAScatalog(LASfile)


opt_output_files(ctg) <- paste0(tempdir(), "/{ID}_hillshade")
opt_chunk_buffer(ctg) <- 30
opt_chunk_size(ctg) <- 150
opt_select(ctg) <- "xyzc"
opt_filter(ctg) <- "-keep_class 2"
opt_merge(ctg) <- FALSE

plot(ctg, chunk = T)


sdtm = hillshade(ctg)
plot(terra::rast(sdtm[[1]]), col = gray(0:50/50))

[wiesehahn]

thanks again for your quick reply.

Your example works for me as well. I will try to find out whats the problem in my case but I guess it has something to do with the catalog.
Running it on part of the collection such as the following works:

LASfile <- ctg@data$filename[some_tilenumber]
ctg = readLAScatalog(LASfile)

However, using the entire catalog and the setting some tile to test on such as the following, results in the error:

# run function for selected tiles
ctg$processed <- FALSE
ctg$processed[some_tilenumber] <- TRUE

I think its strange that th error occours since if some tiles are errorneous its at least not the one I tested on.

Are there any recommended methods to check the validity of a catalog or further restrictions which should be added in the code above?

[Jean-Romain]

I don't know. Is there a traceback with your error that may help me to figure out where does the error come from ? Are you able to identify which file failed ?

[wiesehahn]

OK, I found the problem.
Since its a huge catalog I created the catalog once and stored it as RData object. I created it with lidR <4.0.0 and now read it with lidR >=4.0.0. I wasnt aware there were changes introduced related to the catalog creation.

Recreating the entire catalog solved the problem!
This also explains why my first attempt from above worked were I recreated the catalog for only part of the tile collection to test it.

with the old catalog the following drivers were supported:

drivers     :
 - Raster : format = GTiff  NAflag = -999999  
 - LAS : no parameter
 - Spatial : overwrite = FALSE  
 - SimpleFeature : quiet = TRUE  
 - DataFrame : no parameter

with the new catalog this looks as following:

drivers     :
 - Raster : format = GTiff  NAflag = -999999  
 - stars : NA_value = -999999  
 - SpatRaster : overwrite = FALSE  NAflag = -999999  
 - SpatVector : overwrite = FALSE  
 - LAS : no parameter
 - Spatial : overwrite = FALSE  
 - sf : quiet = TRUE  
 - data.frame : no parameter

[Jean-Romain]

Good ! A collction v3 is expected to be repaired automatically in v4 but I missed to manage the case of drivers.