使用terra软件包对点进行遥感以生成多层SpatRaster

Question

我有一些使用raster::rasterize的长期的包代码，我试图将它更新为terra::rasterize。代码接受点数据，其中每个点都有两个可能的整数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”参数作为数据框架并检查列数来确定的。这有点令人困惑，因为包docs声明，值参数应该是一个数值向量，长度为1或nrow(x)，其中x是输入点数据。此外，用户提供的摘要函数返回的向量的长度在确定输出光栅层数方面似乎没有任何作用。

目前，我只保留了旧的raster::rasterize代码，并将输出光栅转换为SpatRaster，但我想我肯定缺少了一些显而易见的东西。有没有一种使用terra::rasterize来完成这一任务的方法？

编辑:按照的要求在注释中，这里是一个小样本的输入点数据，以显示格式。典型的输入数据大小从200万到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")

编辑：在光栅包代码中，私有.pointsToRaster函数有一行(请看这里)，其中检查用户提供的摘要函数输出的长度和一些任意的测试值，以确定输出光栅中的层数。这似乎没有在地球包裹代码。

Answer 1

这可能是你不希望这两个层在一个栅格，虽然这是很难说提供的数据，因为它似乎都‘在’内的重叠。我注意到在你的包，有一个尝试节流/减少瓦边缘点，可能只是需要设置低于1K。

当terra的工作方式与raster不同时，rasterize(可能是一种决定，即在terra下，人们应该通过将每一层设置为add<-或<- c(来实现两层，而对于raster，则是通过很难遵循“字段”和“值”的逻辑来假定的。使用上面的数据(并保持两个栅格)：

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级别解决这一问题。