요인 수준에 대해 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에 생성됨