使用 R 中的非空间方法将空间点 ID 分配给网格单元答案

【问题标题】：Assign spatial point ID to grid cell using non-spatial method in R使用 R 中的非空间方法将空间点 ID 分配给网格单元
【发布时间】：2021-06-05 22:35:19
【问题描述】：

我有很多点，想根据它们相交的网格单元分配一个 ID。我知道使用sf 包你可以使用st_intersections 或st_intersects 来做到这一点。但是，对于具有数百万个特征/点的大型数据集，这可能需要很长时间（或者只会导致 R 崩溃）。

在投影的 CRS（在我的示例中为英国国家网格）上有一个网格系统意味着存在一致的形状，并且坐标（xmin、ymin、xmax 和 ymax）将始终提供每个单元格的完整范围。是否有使用 tidyverse 或 data.table 中的某些东西的非空间方法，例如可以利用它来分配网格 ID。

这是我的示例数据（如果有人可以向我展示一种更简洁的方法来提取每个网格单元的 xmin、ymin、xmax 和 ymax，我也将不胜感激）：

library(sf)
library(dplyr)

BBox <- st_bbox(c(xmin = 0, xmax = 10000, ymax = 10000, ymin = 0), crs = st_crs(27700))

Grid  <- st_as_sfc(BBox) %>%
         st_make_grid(square = TRUE, cellsize = c(1e3, 1e3)) %>%
         cbind(data.frame(ID = sprintf(paste("GID%0",nchar(length(.)),"d",sep=""), 1:length(.)))) %>%
         st_sf()
             
Points <- st_sample(st_as_sfc(BBox), 3000, exact = TRUE) %>%
                    st_sf('ID' = seq(length(.)), 'geometry' = .) %>%
                    mutate(X = st_coordinates(.)[,1],
                    Y = st_coordinates(.)[,2])      
             
Table <- NULL
for(i in 1:nrow(Grid)) { 
Row <- cbind(as.numeric(st_bbox(Grid[i,])[1]),
    as.numeric(st_bbox(Grid[i,])[2]),
    as.numeric(st_bbox(Grid[i,])[3]),
    as.numeric(st_bbox(Grid[i,])[4]))
Table <- as.data.frame(rbind(Table, Row))
}
names(Table) <- c("xmin","ymin","xmax","ymax")

Grid <- cbind(Grid,Table)

【问题讨论】：

我建议坚持使用空间方法。如果您反对内存限制，您可以将数据卸载到 PostgreSQL 数据库（它非常支持空间数据格式）并在数据库中运行空间连接 = 将内存约束转换为磁盘空间约束。这可能需要一段时间，但您不太可能用完磁盘...

标签： r data.table tidyverse spatial sf

【解决方案1】：

我会（强烈）建议使用像“默认”sf::st_join() 这样的空间方法，您可以轻松地将网格中的 ID 连接到点上

可重现的样本数据

library(sf)
library(dplyr)
set.seed(123)  #<<-- !!

BBox <- st_bbox(c(xmin = 0, xmax = 10000, ymax = 10000, ymin = 0), crs = st_crs(27700))

Grid  <- st_as_sfc(BBox) %>%
  st_make_grid(square = TRUE, cellsize = c(1e3, 1e3)) %>%
  cbind(data.frame(ID = sprintf(paste("GID%0",nchar(length(.)),"d",sep=""), 1:length(.)))) %>%
  st_sf()

Points <- st_sample(st_as_sfc(BBox), 3000, exact = TRUE) %>%
  st_sf('ID' = seq(length(.)), 'geometry' = .) %>%
  mutate(X = st_coordinates(.)[,1],
         Y = st_coordinates(.)[,2])

sf解决方案

ans.sf <- Points %>% st_join( Grid, join = st_intersects )
# Simple feature collection with 6 features and 4 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 455.565 ymin: 1835.024 xmax: 9404.673 ymax: 9425.39
# projected CRS:  OSGB 1936 / British National Grid
#   ID.x        X        Y   ID.y                  geometry
# 1    1 2875.775 2058.269 GID023 POINT (2875.775 2058.269)
# 2    2 7883.051 9425.390 GID098  POINT (7883.051 9425.39)
# 3    3 4089.769 3793.238 GID035 POINT (4089.769 3793.238)
# 4    4 8830.174 6262.401 GID069 POINT (8830.174 6262.401)
# 5    5 9404.673 1835.024 GID020 POINT (9404.673 1835.024)
# 6    6  455.565 6592.076 GID061  POINT (455.565 6592.076)

但是，如果您遇到内存问题，并且正在使用简单的多边形（如网格），则以下方法可能有效（同样，空间解决方案（imo）始终是首选！！）：

data.table 解决方案非空间！！！！

library( data.table )
#easy for ppints, only one xY per row
DT.points <- as.data.table( st_set_geometry( Points, NULL ) )
#    ID        X        Y
# 1:  1 2875.775 2058.269
# 2:  2 7883.051 9425.390
# 3:  3 4089.769 3793.238
# 4:  4 8830.174 6262.401
# 5:  5 9404.673 1835.024
# 6:  6  455.565 6592.076
# ...

#a bit more complcated for polygons
Grid.coords <- Grid %>% 
  st_coordinates() %>%
  as.data.frame() %>%
  group_by( L2 ) %>%
  summarise( minx = min(X, na.rm = TRUE), 
             maxx = max(X, na.rm = TRUE),
             miny = min(Y, na.rm = TRUE),
             maxy = max(Y, na.rm = TRUE) ) 
DT.grid   <- cbind(
  as.data.table( st_set_geometry( Grid, NULL ) ),
  Grid.coords )
#        ID L2 minx maxx miny maxy
# 1: GID001  1    0 1000    0 1000
# 2: GID002  2 1000 2000    0 1000
# 3: GID003  3 2000 3000    0 1000
# 4: GID004  4 3000 4000    0 1000
# 5: GID005  5 4000 5000    0 1000
# 6: GID006  6 5000 6000    0 1000
# ...

#perform non-equi join (assuming a point only cvna fall into 1 polygon!)
DT.points[ DT.grid, Grid.id := i.ID, on = .(X >= minx, X < maxx, Y >= miny, Y < maxy)][]
#         ID         X        Y Grid.id
#    1:    1 2875.7752 2058.269  GID023
#    2:    2 7883.0514 9425.390  GID098
#    3:    3 4089.7692 3793.238  GID035
#    4:    4 8830.1740 6262.401  GID069
#    5:    5 9404.6728 1835.024  GID020
# ---                                
# 2996: 2996 1370.2189 8770.564  GID082
# 2997: 2997 8619.4776 8827.825  GID089
# 2998: 2998  671.2062 9011.131  GID091
# 2999: 2999 4352.0844 7979.424  GID075
# 3000: 3000 6801.3496 9014.823  GID097
#

【讨论】：