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
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="/))