요인 수준에 대해 st_intersection (sf) 반복

Sep 25 2020

셀 내 범주 당 포인트 수를 계산하기 위해 sf 객체를 사용하여 코드를 작성했습니다. 코드는 잘 작동하지만 요인 수준을 반복하고 싶습니다.

다음은 재현 가능한 예입니다.

library(sf)
library(raster)
library(dplyr)

region <- st_read(system.file("external/lux.shp", package="raster"), quiet = TRUE)
PT <- st_sf(st_sample(region, size =100, type='random'))
PT$Class <- sample(LETTERS[1:4]) Mesh <- st_sf(st_make_grid(PT, cellsize = 0.2)) Mesh$ID <- seq.int(nrow(Mesh))

intersection <- st_intersection(Mesh, PT[PT$Class == "A",])
PT_result_A <- intersection %>% 
  group_by(ID) %>% 
  count()

위의 방법은 각 셀 내부에 "A"레이블이있는 포인트 수를 계산하는 데 효과적입니다. 내 원래 코드 / 데이터에는 카테고리가있는만큼이 블록이 있고 카테고리도 많이 있습니다.

마지막 블록을 반복하지 않고 "A", "B", "C"및 "D"의 개수를 얻기 위해 위의 작업을 자동화 할 수 있습니까?

답변

3 agila Sep 25 2020 at 21:17

다음 reprex는 문제에 대한 세 가지 해결책을 제시해야합니다. 마지막은 아마도 불필요하게 복잡 할 수 있지만 시각화가 좋다고 생각합니다.

# packages
suppressPackageStartupMessages({
  library(sf)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
})

# build data
region <- st_read(system.file("external/lux.shp", package="raster"), quiet = TRUE)
points <- st_sf(
  data.frame(class = sample(LETTERS[1:4], size = 100, replace = TRUE)), 
  geometry = st_sample(region, size = 100, type='random')
)
mesh <- st_make_grid(points, cellsize = 0.2)
mesh <- st_sf(
  data.frame(ID = seq_along(mesh)), 
  geometry = mesh
)

# plot data
par(mar = rep(0, 4))
plot(mesh, reset = FALSE, col = sf.colors(nrow(mesh), alpha = 0.25))
plot(points, pch = 16, add = TRUE)

다음을 사용하여 셀당 포인트 수를 계산합니다. tapply

tapply(st_geometry(points), points$class, FUN = function(x) lengths(st_intersects(mesh, x))) #> $A
#>  [1] 1 4 2 0 4 3 0 1 1 6 0 0 0 3 0 0
#> 
#> $B #> [1] 0 3 1 0 1 4 3 0 5 1 0 0 0 0 0 0 #> #> $C
#>  [1] 2 4 3 0 5 1 4 2 5 5 0 0 1 1 0 0
#> 
#> $D
#>  [1] 3 3 1 1 3 1 5 2 2 1 1 0 0 1 0 0

다음을 사용하여 셀당 포인트 수를 계산합니다. dplyr

points %>% 
  group_by(class) %>% 
  group_map(~ lengths(st_intersects(mesh, .x)))
#> [[1]]
#>  [1] 1 4 2 0 4 3 0 1 1 6 0 0 0 3 0 0
#> 
#> [[2]]
#>  [1] 0 3 1 0 1 4 3 0 5 1 0 0 0 0 0 0
#> 
#> [[3]]
#>  [1] 2 4 3 0 5 1 4 2 5 5 0 0 1 1 0 0
#> 
#> [[4]]
#>  [1] 3 3 1 1 3 1 5 2 2 1 1 0 0 1 0 0

dplyr+ group_nest카운트를 표시하는지도를 만들려면

counts_per_cell_and_class <- points %>% 
  group_by(class) %>% 
  group_nest() %>% 
  mutate(counts = lapply(data, function(df) as.character(lengths(st_intersects(mesh, df))))) %>% 
  select(-data) %>% 
  unnest(cols = counts) %>% 
  mutate(ID = rep(seq_len(nrow(mesh)), 4))

ggplot(st_as_sf(right_join(counts_per_cell_and_class, mesh))) +  
  geom_sf(aes(fill = counts)) + 
  geom_sf(data = points) + 
  facet_wrap(~ class, nrow = 1, ncol = 4) + 
  scale_fill_brewer() + 
  theme(legend.position = "bottom")
#> Joining, by = "ID"

reprex 패키지 (v0.3.0)로 2020-09-25에 생성됨