Expanded example of joining spatial datasets, counting intersecting points across shapes, and using as an aesthetic

st_join_example.R

# refer to https://gis.stackexchange.com/questions/110117/counting-number-of-points-in-polygon-using-r/270479#270479
# for original code example

# Load libraries ----------------------------------------------------------

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

# Get sample data ---------------------------------------------------------

# Get polygon
polygon <- getData('GADM', country='URY', level = 1)[,1] # Download polygon of country admin level 1 
polygon <- st_as_sf(polygon) # convert to sf object
colnames(polygon) <- c("id_polygons", "geometry") # change colnames
polygon$id_polygons <- paste0("poly_", LETTERS[1:19]) #  change polygon ID

# Get sample random points from polygon bbox
set.seed(4)
bbox <- st_as_sfc(st_bbox(polygon))
points <- st_sample(x = bbox, size = 100, type = "random")
points <- st_as_sf(data.frame(id_points = as.character(1:100)), points) # add points ID

# Plot data ---------------------------------------------------------------
# plot.new()

# Plot polygon + points
# plot(points, pch = 19, col = "black", graticule = st_crs(4326), add = TRUE)
# plot(polygon, graticule = st_crs(4326), key.pos = 1)

ggplot() +
  geom_sf(data = polygon) +
  geom_sf(data = points)

# Intersection between polygon and points ---------------------------------

intersectionPolygonPoints <- st_intersection(x = polygon, y = points)

# Plot intersection

# plot(polygon, graticule = st_crs(4326), key.pos = 1)
# plot(intersectionPolygonPoints[1], col = "black", pch = 19, add = TRUE)

ggplot() +
  geom_sf(data = polygon) +
  geom_sf(data = points) +
  geom_sf(data = intersectionPolygonPoints[1])

# View result using base R cross tabulation function
(crosstab1 <- table(intersectionPolygonPoints$id_polygons)) # using table

# using dplyr
(int_result <- intersectionPolygonPoints %>% 
  group_by(id_polygons) %>% 
  count())

as.data.frame(int_result)[,-3]

(crosstab2 <- polygon %>% st_join(int_result))

ggplot() +
  geom_sf(data = crosstab2, aes(fill = n)) +
  geom_sf(data = points)