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

【问题标题】：How to efficiently count the number of spatial points within a certain distance around raster cells in R?如何有效计算R中栅格单元周围一定距离内的空间点数？
【发布时间】：2020-05-18 09:45:56
【问题描述】：

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

# 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 组成。有没有一种方法可以更有效地生成简单地计算每个栅格单元周围一定距离内有多少点的栅格？

【问题讨论】：

标签： r vector distance spatial raster

【解决方案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-format（光栅类型格式——可以用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 做一些优雅而快速的事情，但我不适合在那里问。

【讨论】：