r-使用terra包对点进行光栅化，以生成多层SpatRaster

Question

我有一些使用raster:：rasterize的长期包代码，我正试图将其更新为terra:：raster。该代码获取点数据，其中每个点都有两个可能的整数ID值中的一个。输出是一个具有两层的光栅，每个可能的点ID对应一层，其中单元值是计数。相关位为：

# r0 is template raster to define extent and resolution
r <- raster::rasterize(dat[, c("X", "Y")],
r0,
field = dat$flightlineID,
fun = f,
background = 0)

这里，f是一个函数，它获取点ID的矢量，并返回计数的两元素矢量，从而产生所需的两层输出光栅。

我第一次尝试将其移植到terra:：rasterize(包版本1.6-17(是。。。

r <- terra::rasterize(cbind(dat$X, dat$Y), # seem to need a matrix rather than a data frame
r0,  # template SpatRaster
values = dat$flightlineID, 
fun = f, 
background = 0)

此操作失败，错误为：

w[vv[，1]，]<-中的错误vv[，-1]：要更换的物品数量不是更换长度的倍数

深入研究terra:::rasterize_points的代码，输出光栅的层数似乎是通过将"values"参数视为数据帧并检查列数来确定的。这有点令人困惑，因为包文档声明values参数应该是一个长度为1或nrow(x(的数字向量，其中x是输入点数据。此外，用户提供的汇总函数返回的矢量长度似乎在确定输出光栅层数方面没有任何作用。

目前，我只是保留了旧的光栅：：光栅化代码并将输出光栅转换为SpatRaster，但我认为我一定缺少了一些明显的东西。有没有一种方法可以只使用terra:：光栅化来完成这项任务？

编辑：根据注释中的要求，这里有一个输入点数据的小样本来显示格式。典型的输入数据大小在2到4000万个点之间。

structure(list(X = c(420094, 420067, 420017, 420050, 420058, 
420090, 420038, 420040, 420081, 420097, 420075, 420041, 420039, 
420062, 420050, 420083, 420019, 420019, 420044, 420087, 420099, 
420077, 420030, 420014, 420015, 420051, 420033, 420056, 420041, 
420030, 420027, 420024, 420058, 420042, 420063, 420028, 420073, 
420053, 420010, 420100, 420048, 420062, 420056, 420080, 420053, 
420068, 420074, 420004, 420010, 420078), Y = c(6676049, 6676029, 
6676034, 6676019, 6676096, 6676010, 6676003, 6676048, 6676073, 
6676023, 6676089, 6676082, 6676010, 6676051, 6676039, 6676099, 
6676024, 6676073, 6676040, 6676056, 6676072, 6676086, 6676030, 
6676042, 6676002, 6676033, 6676078, 6676073, 6676013, 6676056, 
6676055, 6676069, 6676072, 6676089, 6676069, 6676058, 6676023, 
6676039, 6676043, 6676017, 6676011, 6676054, 6676095, 6676068, 
6676098, 6676077, 6676049, 6676073, 6676097, 6676057), flightlineID = c(2L, 
1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 
2L)), row.names = c(NA, -50L), class = "data.frame")

EDIT:在光栅包代码中，专用.pointsToRaster函数有一行(请参见此处(，其中使用一些任意测试值检查用户提供的摘要函数的输出长度，以确定输出光栅中的层数。这似乎在terra包代码中不存在。

Answer 1

您可能不希望将其作为一个光栅中的两层，尽管这与所提供的数据很难区分，因为它似乎都"在"重叠范围内。我注意到在你的包中，有一种尝试节流/减少瓦片边缘点，可能只需要将其设置为低于1K。

当rasterize(ing时，terra的工作原理与raster不同，这可能是一个决定，即在terra下，应该通过使每个层都成为add<-ing或<- c(ing来达到两层，而在raster中，它是通过"字段"one_answers"值"的难以遵循的逻辑来假设的。使用以上数据(并保留两个光栅(：

library(terra) 
#las_df <- structure(...)
las_df1 <- las_df[which(las_df$flightlineID == 1L), ]
las_df2 <- las_df[which(las_df$flightlineID == 2L), ]
las_vect1 <- vect(las_df1, geom = c('X', 'Y'), crs = 'EPSG:32755')
las_vect2 <- vect(las_df2, geom = c('X', 'Y'), crs = 'EPSG:32755')
las_rast <- rast(xmin=0, nrow = length(unique(las_df$X)), ncol = length(unique(las_df$Y)), crs='EPSG:32755')
set.ext(las_rast, c(min(las_df$X), max(las_df$X), min(las_df$Y), max(las_df$Y)))
pts1_rast <- rasterize(las_vect1, las_rast, fun = length)
pts2_rast <- rasterize(las_vect2, las_rast, fun = length)
pts1_pts2_rast <- c(pts1_rast, pts2_rast)
names(pts1_pts2_rast) <- c('lyr.1', 'lyr.2') # have to attend to this as both lyr.1 after `c(`
plot(pts1_pts2_rast$lyr.1, col = 'red')
plot(pts1_pts2_rast$lyr.2, col = 'blue', alpha=.75, add = TRUE)
# there is 1 cell that contains points from both pts1_rast and pts2_rast
cells(pts1_rast) %in% cells(pts2_rast)
[1] FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
cells(pts2_rast) %in% cells(pts1_rast)
[1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE FALSE FALSE

可以建议一致的merge策略，其中pts1或pts2总是受到青睐。最后，如果这是为了优化稀缺资源的分配，请清理您拥有最佳数据的灌木丛，进行检查，然后再次清理。但似乎最好还是在上游的las级别解决这个问题。