normalize intensities across and between tiles for LAScatalog

[wiesehahn]

My aim is to create some intensity metrics as raster. I want to normalize values within tiles but also try to normalize the range bewteen tiles from different flight campaigns (at least approximately).

I use the following code but have the problem that no files are written to the output location.

library(here)
library(lidR)

load(here::here("data/interim/lidr-catalog-new.RData"))

storage_path <- here("data/processed")
opt_output_files(ctg) <- paste0(storage_path, "/{XLEFT}_{YBOTTOM}_iupper")
ctg@output_options$drivers$SpatRaster$param$gdal = c("COMPRESS=DEFLATE", "PREDICTOR=2")
ctg@output_options$drivers$SpatRaster$param$datatype = "INT1U"
ctg@output_options$drivers$SpatRaster$param$NAflag = 0


opt_chunk_buffer(ctg) <- 200
opt_chunk_size(ctg) <- 0
opt_merge(ctg) <- FALSE
opt_stop_early(ctg) <- FALSE



myMetric<-function(i, z){

  # select upper 25% and 
  q75=quantile(z,probs=c(0.75))
  aboveq75= z>q75

  # calculate mean intensity
  zq75 = i[aboveq75]
  imeanq75=mean(zq75, na.rm=TRUE)

  return(imeanq75)
  }


create_imean <- function(las, resolution = 1)
{
  if (is(las, "LAScatalog"))  {
    options <- list(automerge = FALSE, need_buffer = TRUE)
    opt_select(las) <- "xyzci"
    # readLAS(filter = "-h")
    opt_filter(las) <- "-keep_class 2 -keep_class 13 -keep_class 20 -drop_z_above 1200 -drop_z_below -3 -drop_number_of_returns 0 -drop_gps_time_below 0"
    return(catalog_apply(las, create_imean, resolution=resolution, .options = options))
  }
  else if (is(las, "LAScluster")) {
    las <- readLAS(las)
    if (is.empty(las)) return(NULL) 
    imean <- create_imean(las, resolution=resolution)
    imean <- terra::crop(imean, ext(chunk))
    return(imean)
  }
  else if (is(las, "LAS")) {
    # calculate sensor positions
    sensor <- track_sensor(las, thin_pulse_with_time = 0.001, Roussel2020(interval = 0.5, pmin = 50))
    
    # normalize intensity
    lasnorm <- normalize_intensity(las, range_correction(sensor, Rs = 500, f = 2.3))
    
    # reduce to 1-99% interval
    lasnorm$Intensity <- ifelse(lasnorm$Intensity < quantile(lasnorm$Intensity, probs= 0.01), 
                   as.integer(quantile(lasnorm$Intensity, probs= 0.01)),
                   
                   ifelse(lasnorm$Intensity > quantile(lasnorm$Intensity, probs= 0.99), 
                          as.integer(quantile(lasnorm$Intensity, probs= 0.99)),
                          lasnorm$Intensity))
                   
    # normalize to 8bit range
    lasnorm$Intensity = as.integer(lasnorm$Intensity/max(lasnorm$Intensity)*255L)
    
    # normalize height
    nlas <- normalize_height(lasnorm, tin())
    # select above DBH
    noground <- filter_poi(nlas, Classification != 2, Z >1.3)
    
    imean <- pixel_metrics(noground, ~myMetric(Intensity, Z), res = resolution)
    names(imean) <- "Intensity"
    
    return(imean)
  }
  else {
    stop("This type is not supported.")
  }
}


create_imean(ctg, 1)

I tried to find the error for quite a long time until I just recognized the Non-supported LAScatalog options of track sensor.

I dont really understand what is meant with

The option 'output_files' is not supported and not respected because the output must be post-processed as a whole
but I am pretty sure that this is the case why I cant find the files.

I guess the files are not written at all, or are they just written somewhere else where I cant find them?

And if they are not written what would I have to do to process a catalog to normalize intensities?

[Jean-Romain]

I don't really understand what is meant with

When applying a method on a collection of file by chunks with a buffer you necessarily compute some stuff twice or more in some regions. For example, for a DTM, you compute the DTM twice or more in the buffered regions because the buffered region are loaded more than once (they overlap). For rasters the buffer is removed on-the-fly based on easy spatial criterion. If a pixel is in the buffer we remove it. Job done!

For track_sensor this is different. For a given chunk + buffer there is no way to remove the buffer. Actually it is meaningless. We split the dataset in spatial chunks with spatial buffers but in theory we should better split the dataset in temporal chunks with temporal buffer (based on gpstime). A valid sensor position within the temporal bounding box of a chunk can be out of the spatial bounding box of the chunk because of the different angles involved.

After processing, for each chunk we have some sensor position that are duplicated with other chunks. And each sensor position is not estimated with the same score of accuracy. There is no way to produce a clean and continuous output by processing each chunk independently due to the fact we are processing spatial chunk instead of temporal chunk. Consequently, at the end, all the results are merged and duplicated sensor position are removed by keeping the best one.

This means we cannot process each spatial chunk independently. We need to see what is in the other chunks to produce a clean an continuous output. This is why it is not possible to write the result of each chunk independently.

In your case, if you have a buffer large enough (200m is not large enough), your code looks legit. You compute the sensor, the sensor is poor on the edge, consequently the normalization is poor on the edge but who care because you rasterize later and clip the spatial buffer on the edges. You don't keep the sensor positions anyway.

If your intensity raster is not written it comes from somewhere else.

[wiesehahn]

Ok, I understand why We need to see what is in the other chunks to produce a clean and continuous output.

As you say I just calculate sensor positions as an intermediate step to normalize intensities and then output the clipped raster. I am aware that many points are used multiple times with tile by tile calculations and a large buffer, and I am also aware that the calculated sensor positions may vary a little between tiles if we dont process all points of a certain gpstime. But my assumption would be that the unperfectly normalized points are a better indicator for intensities than the unnormalized points (maybe I am wrong).

Do I understand your answer correct that my code and thus also the intermediate usage of track_sensor() applied on a LAScatalog should work under normal circumstances?

I have no idea what is going on under the hood with

The option 'output_files' is not supported and not respected

But running track_sensor() inside the function prevents the writing of output to 'output_files', although the result of track_sensor() is not even used.

Changing part of the function above to something without track_sensor() works as expected, e.g.:

else if (is(las, "LAS")) {
       
    nlas <- normalize_height(las, tin())
    noground <- filter_poi(nlas, Classification != 2, Z >1.3)
    
    imean <- pixel_metrics(noground, ~myMetric(Intensity, Z), res = resolution)
    names(imean) <- "Intensity"
    
    return(imean)
  }

While it does not work as expected when introducing track_sensor():

else if (is(las, "LAS")) {
     # try with sensor tracking although not using the result
     sensor <- track_sensor(las, thin_pulse_with_time = 0.001, Roussel2020(interval = 0.5, pmin = 50))
  
    nlas <- normalize_height(las, tin())
    noground <- filter_poi(nlas, Classification != 2, Z >1.3)
    
    imean <- pixel_metrics(noground, ~myMetric(Intensity, Z), res = resolution)
    names(imean) <- "Intensity"
    
    return(imean)
  }

Is this the expected behavior?
And if yes, how would I normalize intensities for a LAScatalog?

[Jean-Romain]

Absolutely not. This is weird and deserve a bug report with a mininal reproducible exemple. In your code, what could be written in files is the imean for each chunk. The output of track_sensor() is only an intermediate object in memory.

[wiesehahn]

I just tried to make a reproducible example and figured out that the problem is not the fact that 'output_files' is not supported.
Instead it is (again) the fact that my data is imperfect. I tried to run the intensity normalization on several separate LASfiles before and it worked, but now I figured out that some of them are problematic and when applied on a catalog the errorneous pulses are taken into account for the othwerwise functioning LASfile as well because of the buffer.
This leads to errors from the track_sensor() function which is preventing the result from being written to disk although the sensor position is not even used in the example.

The error I get is
Error: After keeping only first and last returns of multiple returns pulses, 26 pulses still have more than 2 points. This dataset is corrupted and gpstime is likely to be invalid.

Since it is only a relatively small amount of pulses I think It might be reasonable to simply filter these pulses. Is there a decent method to do this?

[Jean-Romain]

It may or may not be only 26 pulses. I mean, why 26 pulses triggered an error but all other looks good? Where are these pulses? Is there a spatial pattern? Is your dataset corrupted? Have those pulse other issues? Is is a specific flight-line? Is is random? Is it a special file sent by the data provider after the main project? So many questions before to simply get rid of the points

To get rid of those point you must redo the following steps: identify pulses based on gpstime, keep first and last, count the number of points per pulse, and remove the points for which the pulse ID is in the list of ID where the number of points is more than 2. No built in tools.