# Computing the density for each layer with lidR

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
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)

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
}
``````

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

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)
``````

• I am sorry but I am pretty new in r even programming. Why are you using as.raster on an object that is already a RasterLayer? -because I wanted to create a base layer to use in the calculation. Where the errors occurred? -creating as.raster and applying grid density in function gives error which is written above. Which data did you used? I used las data which is already ground classified. @JRR Commented Feb 23, 2021 at 15:03
• Please marked the question as answered (click on the checkmark) if it is the case. If not please edit your question to clarify where you are still struggled.
– JRR
Commented Feb 23, 2021 at 15:14
• This code works without any problem! Just want to be sure, when I change the input to my las file, all the output rasters have min and max values as zero. Is that because of the input or should I tune anything? Commented Feb 23, 2021 at 15:23
• I can't know. It depends on your data. If you used 5 cm slices on a 2-3 pt/m2 point cloud I'd say yes many slices are likely to be empty. But you can't ask if the output is correct without describing the input.
– JRR
Commented Feb 23, 2021 at 15:26