Manually load large raster to perform lidR::locate_trees()

Question

I am attempting perform individual tree detection on a lidar-derived canopy height model raster using the lidR package in R. The raster is approximately 2.5 Gb. When I attempt to run the locate_trees() function, I get this error: Error: Large on-disk rasters are not supported by locate_tree. Load the raster manually. What is meant by "loading the raster manually?" and how can I do this? I made a reproducible example of my workflow. However, making a fake >3GB raster may crash your R session or take forever. Use caution before running:

library(lidR)
library(terra)

chm<-rast(xmin=-121.5, xmax =-120, ymin=44.5, ymax = 45.5, crs="EPSG:2992", nrows=150000, ncols=200000) #Fake 3 billion cell raster
values(chm)<-rnorm(3000000000, 20, 5) #giving random values to the raster cells

writeRaster(chm, paste(tempdir(), "CHM_Example.tif", sep="/"), filetype="GTiff")#writing the raster to disk

canopy<-rast(paste(tempdir(), "CHM_Example.tif", sep="/")) # reading the raster into R using the terra package framework

taos<-locate_trees(canopy, lmf(ws=5))#Error is thrown here attempting to perform individual tree detection

Answer 1

It dawned on me that I had run into this issue before...and solved it. I dug up and dusted off an old script from a previous project, and voila the answer was there all along! "Loading the raster manually" means actually using the readStart(), readValues(), and readStop() functions in terra to manually read all values within predetermined block of the raster to then feed to lidR::locate_trees(). I adapted the chunk and buffer scheme that @JRR in lidR into a loop to ensure I am getting full coverage without duplicates. Admittedly, this is probably not the most efficient solution, but it works.

N.B. in this reproducible example, I just use the Mixed Conifer dataset that comes with the lidR package

library(terra)
library(sf)
library(lidR)

LASFile<-system.file("extdata", "MixedConifer.laz", package="lidR")
las<-readLAS(LASFile)

dem<-rasterize_terrain(las)
norm<-las-dem
chm<-rasterize_canopy(norm)

writeRaster(chm, paste(tempdir(), "CHM_Example.tif", sep="/"), filetype="GTiff")#writing the raster to disk

Canopy<-rast(paste(tempdir(), "CHM_Example.tif", sep="/")) # reading the raster into R using the terra package framework

rws<-seq(1, nrow(Canopy), 5)
rws<-c(rws, nrow(Canopy))
cls<-seq(1, ncol(Canopy), 5)
cls<-c(cls, ncol(Canopy))
chunks_r<-length(rws)
chunks_c<-length(cls)
buff<-2
K<-1
cref<-crs(Canopy)


readStart(Canopy)

for(i in 1:(chunks_c-1)){
  
  xmn<-xFromCol(Canopy, cls[i])
  xmx<-xFromCol(Canopy, cls[i+1])
  
  if (i == 1){
    ncls<-(cls[i+1]+buff-cls[i])
    cread<-cls[i]
    xmnb<-xFromCol(Canopy, cls[i])
    xmxb<-xFromCol(Canopy, cls[i+1]+buff)
  }
  else{
    
    if (i == (chunks_c-1)){
      ncls<-(cls[i+1]+buff-cls[i])
      cread<-cls[i]-buff
      xmnb<-xFromCol(Canopy, cls[i]-buff)
      xmxb<-xFromCol(Canopy, cls[i+1])
    }
    else{
      
      ncls<-(cls[i+1]+buff+buff-cls[i])
      cread<-cls[i]-buff
      xmnb<-xFromCol(Canopy, cls[i]-buff)
      xmxb<-xFromCol(Canopy, cls[i+1]+buff)
    }
  }  
  for(j in 1:(chunks_r-1)){
    
    ymx<-yFromRow(Canopy, rws[j])
    ymn<-yFromRow(Canopy, rws[j+1])
    
    if (j == 1){
      nrws<-(rws[j+1]+buff-rws[j])
      rread<-rws[j]
      ymxb<-yFromRow(Canopy, rws[j])
      ymnb<-yFromRow(Canopy, rws[j+1]+buff)
    } 
    else{
      
      if (j == (chunks_r-1)){
        
        nrws<-(rws[j+1]+buff-rws[j])
        rread<-rws[j]-buff
        ymxb<-yFromRow(Canopy, rws[j]-buff)
        ymnb<-yFromRow(Canopy, rws[j+1])
      }
      else{
        
        nrws<-(rws[j+1]+buff+buff-rws[j])
        rread<-rws[j]-buff
        ymxb<-yFromRow(Canopy, rws[j]-buff)
        ymnb<-yFromRow(Canopy, rws[j+1]+buff)
      }
    }
    
    vals<-terra::readValues(Canopy, start=1, row=rread, nrows=nrws, col=cread, ncols=ncls)
    
    buff_ras<-rast(ymin=ymnb, ymax=ymxb, xmin=xmnb, xmax=xmxb, nrows=nrws, ncols=ncls, crs=cref, vals=vals)
    blank<-rast(ymin=ymn, ymax=ymx, xmin=xmn, xmax=xmx, nrows=nrws, ncols=ncls, crs=cref, vals=NA)
    blank_sf<-st_as_sf(st_as_sfc(st_bbox(blank)))
    
    if(any(!is.nan(vals))){
      out<-locate_trees(buff_ras, lmf(3, hmin = 10), uniqueness = "gpstime")
      ttops<-st_crop(out, blank_sf)
      fname<-paste(paste(tempdir(),"Gilchrist_itd_", sep="/"), K, ".gpkg", sep="")
      st_write(obj=ttops, dsn=fname, driver="GPKG", layer="tree_approx_obj", append=FALSE)
      K<-K+1
    }else
    {
      fname<-NULL
    }
    print(c(i, j,K))
  }
}

readStop(Canopy)

tmp<-list.files(tempdir(), pattern="*.gpkg", full.names = TRUE)

fxn<-function(X, fname){
  gpkg<-st_read(X)
  st_write(obj=gpkg, dsn=fname, layer="tree_approx_obj", append=TRUE)
}

taos<-lapply(tmp, fxn, fname=paste(tempdir(), "Example.gpkg", sep="/))

Answer 2

The message means that lidR determined that it cannot load all the raster data into memory. One way to work around it is to use "tiles" (subsets), for example with terra::makeTiles

Example data

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)

Make tiles. Do not use "30" for a big dataset like yours!

tiles <- makeTiles(r, 30)

If the output data is not very large you can do

out <- sapply(tiles, \(tile) {
        x <- rast(tile)
        taos <- lidR::locate_trees(x, lidR::lmf(ws=5))
        vect(taos)
    })

out <- vect(out)

Otherwise, it may be wise to write the output for each tile to disk. For example:

fout <- gsub(".tif", "_trees.shp", tiles)

for (i in 1:length(tiles)) {
    x <- rast(tiles[i])
    taos <- lidR::locate_trees(x, lidR::lmf(ws=5))
    writeVector(vect(taos), fout[i])
}