gis r – Calcular la densidad de cada capa con lidR

Pregunta:

Estoy intentando realizar un índice de perfil 3D como función. Para esto, necesito cortar las capas y calcular la densidad de puntos para cada capa cortada. Escribo un bucle for;

El problema es

density=grid_density(subset,res=0.1)

y

pt<-as.raster(pt)

partes da siempre un error. Cuando ejecuto as.raster siempre obtengo este error:

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

as.raster() raster() lugar de as.raster() , funciona pero no sé cuál es la principal diferencia.

Filtro los puntos con filter_poi() . En este filtro, especifiqué el valor z, los puntos deben calcularse más bajos que la siguiente capa cortada y más altos que la capa en sí. Entonces quiero calcular el grid_density() da este 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
} 

Respuesta:

No entiendo su problema porque no está claro y su pregunta es confusa. ¿Por qué está utilizando as.raster en un objeto que ya es un RasterLayer ? ¿Dónde ocurrieron los errores? ¿Qué datos usaste? Entonces, hice un ejemplo que parece ser lo que estás buscando pero sin ninguna explicación porque no sé cuál es tu problema:

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)

ingrese la descripción de la imagen aquí

Leave a Comment

Your email address will not be published.

Scroll to Top

istanbul avukat

-

web tasarım