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
}