5

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

1 Answer 1

4

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)

enter image description here

4
  • 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
    – RS_girl
    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?
    – RS_girl
    Commented Feb 23, 2021 at 15:23
  • 4
    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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.