Iterar st_intersection (sf) sobre niveles de factor

Question

Escribí un código usando objetos sf para contar el número de puntos por categoría dentro de las celdas. El código funciona bien, pero me gustaría iterar sobre los niveles de los factores.

Aquí hay un ejemplo reproducible:

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()

Lo anterior funciona bien para contar el número de puntos con la etiqueta "A" dentro de cada celda. En mi código / datos originales, tengo este bloque tantas veces como categorías, y tengo muchas categorías.

¿Es posible automatizar lo anterior para obtener el recuento de "A", "B", "C" y "D" sin repetir el último bloque?

Answer 1

El siguiente reprex debería presentar tres soluciones a su problema. El último es probablemente innecesariamente complicado, pero creo que la visualización es agradable.

# 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)

Cuente el número de puntos por celda usando 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

Cuente el número de puntos por celda usando 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_nestpara crear un mapa que muestre los conteos

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"

^{Creado el 2020-09-25 por el paquete reprex (v0.3.0)}