2

I wrote a code using sf objects in order to count the number of points per category within cells. The code works well, but I would like to iterate over factor levels.

Here is a reproducible example :

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

The above works well to count the number of points with the "A" label inside each cell. In my original code/data I have this block as many time as I have categories, and I have a lot of categories.

Is it possible to automate the above to obtain the count of "A", "B", "C" and "D" without repeating the last block?

2
  • 1
    Just use tapply with length as the specified function. – Jeffrey Evans Sep 25 '20 at 0:22
  • Thx, I gave tapply some tries but with no luck so far – ePoQ Sep 25 '20 at 9:15
3

The following reprex should present three solutions to your problem. The last one is probably unnecessarily complicated but I think that the visualisation is nice.

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

Count the number of points per cell using 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

Count the number of points per cell using 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 to create a map displaying the counts

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"

Created on 2020-09-25 by the reprex package (v0.3.0)

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.