Computing the density for each layer with lidR

Question

I am trying to perform a 3D profile index as a function. For this, I need to slice the layers and calculate the point density for each sliced layer. I write a for loop;

The problem is

density=grid_density(subset,res=0.1)

and

pt<-as.raster(pt)

parts gives always an error. When I run as.raster I always get this error:

Error in UseMethod("as.raster") : 
  no applicable method for 'as.raster' (applied to an object of class "c ('LAS', 'Spatial')")

I tried raster() instead of as.raster() it works but I don't know what is the main difference.

I filter the points with filter_poi() . In this filter I specified z value, the points must be calculated lower than the next sliced layer and higher than the layer itself. Then I want to calculate the grid_density() it gives this error:

Error in validityMethod (object): invalid extent: xmin> = xmax
Additionally: Warning messages:
1: In min (x @ data $ X): No missing arguments for min; Returning inf
2: In max (x @ data $ X):
   No complete arguments for max; Returning inf
3: In min (x @ data $ Y): No missing arguments for min; Returning inf
4: In max (x @ data $ Y):
   No missing arguments for max; Returning inf

library(lidR)
library(raster)
library(sp)
library(rgdal)
library(gridExtra)

#set the parameters
sz=0.05 #slice size in meters
s=seq(0,2,sz) #number of layer 
k=seq(-1, 2.00,0.05) #correction factor (k)

las <- readLAS(x)
pt=grid_density(las,res=0.1) #pt is the total LiDAR points for all layers  
pt<-as.raster(pt) #convert to raster
crs(pt) <- "+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs"

# #make base layer, used to maintain the extent of the Pi raster layers
base_layer<-pt
values(base_layer)<-0

density_rasters=c()
for (i in s){
  subset=filter_poi(las, Z>=i & Z<(i+sz))
  if (!is.null(subset)){
    density=grid_density(subset,res=0.1)
    ras=as.raster(density)
    setExtent(ras, ext=extent(base_layer),keepres = TRUE, snap = TRUE) # set extent to baselayer, otherwise the extent will be to small to be able to stack the rasters
    ras=merge(ras,base_layer)
    density_rasters=c(density_rasters,ras)
   
  }
  else { #used to insert empty raster if no points are within a given pi layer
    density_rasters=c(density_rasters,base_layer)
    
  }
  density=0
  ras=0
  subset=0
} 

Answer 1

I don't understand you problem because it is unclear and your question is messy. Why are you using as.raster on an object that is already a RasterLayer? Where the errors occurred? Which data did you used? So, I made an example that seems to be what you are looking for but without any explanation because I don't know what your problem is:

library(lidR)

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

sz = 2 # slice size in meters
s = seq(0, 20, sz) # number of layer 

layout = grid_density(las, res = 4)
layout[] <- 0

density_rasters = vector("list", length(s))
for (i in seq_along(s))
{
  subset = filter_poi(las, Z >= s[i] & Z < (s[i]+sz))
  
  if (!is.empty(subset))
    density_rasters[[i]] = grid_density(subset, res = layout)
  else
    density_rasters[[i]] = layout
}

density_rasters <- lapply(density_rasters, function(x) { extend(x, extent(layout))})
density_rasters <- stack(density_rasters)
names(density_rasters) <- paste0("Layer ", s)
plot(density_rasters)