如何有效地计算R中栅格像元周围一定距离内的空间点数？

Question

我想计算R中SpatialPointsDataFrame的每个像元一定距离内（RasterLayer对象的空间点的数量。结果值应替换该特定栅格像元的原始值。这是一个可重现的示例：

# load library
library(raster)

# generate raster
ras <- raster(nrow=18, ncol=36)
values(ras) <- NA

# create SpatialPointsDataFrame
x <- c(-160,-155,-153,-150, 30, -45, -44, -42, -40, 100, 110, 130)
y <- c(-75,-73,-71,-60, 0, 30, 35, 40, 41, 10, -10, 60)
z <- c(seq(1, 12, 1))
df <- data.frame(x,y,z)
spdf <- SpatialPointsDataFrame(coords=df[,c(1,2)],
                               data=as.data.frame(df[,3]), 
                               proj4string=CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
# visualize
plot(ras)
plot(spdf, add=T)

# loop over all raster cells
for(r in 1:nrow(ras)){
  for(c in 1:ncol(ras)){
    # duplicate raster for subsequent modification
    ras_x <- ras
    # define cell for which to count the number of surrounding points
    ras_x[r,c] <- nrow(spdf) # some value that is impossible to be true, this is only a temporary placeholder
    ras_x[ras_x != nrow(spdf)] <- NA
    # convert raster cell to spatial point
    spatial_point <- rasterToPoints(ras_x, spatial=T)
    # calculate distance around raster cell
    ras_dist <- distanceFromPoints(ras_x, spatial_point)
    ras_dist <- ras_dist / 1000000 # scale values
    # define circular zone by setting distance threshold (raster only with values 1 or NA)
    ras_dist[ras_dist > 2] <- NA
    ras_dist[ras_dist <= 2] <- 1

    # create empty vector to count number of spatial points located within zone around the particular raster cell
    empty_vec <- c()
    # loop to check which value every point of SpatialPointsDataFrame corresponds to 
    for (i in 1:nrow(spdf)){
      point <- extract(ras_dist, spdf[i,])
      empty_vec[i] <- point
    }
    # sum of resulting vector is the number of points within surrounding zone around predefined raster cell
    val <- sum(na.omit(empty_vec))
    val
    ras[r,c] <- val

    # print for progress monitoring
    print(paste0("sum of points within radius around cell row ", r, " and column ", c, " is ", val))
    print(paste0("finished ", r, " out of ", nrow(ras)))
    print(paste0("finished ", c, " out of ", ncol(ras)))
    # both plots are just for visualization and progress monitoring
    plot(ras)
    plot(spdf, add=T)
  }
}

plot(ras)
plot(spdf, add=T)

生成的栅格正是我想要的，但是我检查SpatialPointsDataFrame的每个点的基础栅格值的方式似乎效率很低。我的真实数据包括具有2160、4320、9331200（nrow，ncol，ncell）的RasterLayer和具有2664特征的SpatialPointsDataFrame。有没有一种方法可以简单地更有效地计算每个栅格像元周围某个距离内的点数，从而生成栅格？

Answer 1

如果您可以使用投影坐标，则可以使用spatstat包轻松完成。这要求您使用例如来投影点（和网格）。 sf::st_transform()，将无法使用在全球范围内。

加载spatstat并获得2000个随机点来进行测试：

library(spatstat)
W <- square(1)
set.seed(42)
Y <- runifpoint(2000) # Random points in the unit square
plot(Y, main = "Random points in unit square")

获得3000x3000的点数网格（900万点）：
xy <- gridcenters(W, 3000, 3000) # Grid of points in the unit square
X <- ppp(xy$x, xy$y, window = W, check = FALSE, checkdup = FALSE)

对于900万个网格点中的每一个，计算其中的其他点数半径0.01（在具有16GB RAM的快速笔记本电脑上计时）：
system.time(counts <- crosspaircounts(X, Y, r = .01))
#>    user  system elapsed 
#>   1.700   0.228   1.928

转换为spatstat的im格式（栅格类型格式–可以使用maptools转换）并绘制：
rslt <- as.im(data.frame(x = xy$x, y = xy$y, counts))
plot(rslt, main = "Point counts in raster cells")

叠加在计数上的点表示我们已做正确的事：
plot(rslt, main = "Point counts in raster cells")
plot(Y, add = TRUE, col = rgb(1,1,1,.7), pch = 3)

我敢肯定，您也可以使用raster做一些优雅而又快速的事情，但我不适合在这里问这个问题。