¿Contando el número de puntos en el polígono usando R?

Tengo dos clases que comparten el mismo CRS (Latitute y Longitude):

1. bolognaQuartieriMap: SpatialPolygonDataFramedatos que contienen los municipios de una ciudad.
2. crashPoints: a que SpatialPointsDataFramecontiene datos de accidentes.

Están bien trazados usando:

plot(bolognaQuartieriMap)
title("Crash per quartiere")
plot(crashPoints, col="red",add=TRUE)

Lo que necesito es obtener el número de puntos ( crashPoints) en cada polígono que constituyen bolognaQuartieriMap. Me sugirieron usar over()pero no tuve éxito.

Como no proporcionó un ejemplo reproducible ni un mensaje de error, vea si este fragmento de código le ayuda a comenzar:

library("raster")
library("sp")

x <- getData('GADM', country='ITA', level=1)
class(x)
# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"

set.seed(1)
# sample random points
p <- spsample(x, n=300, type="random")
p <- SpatialPointsDataFrame(p, data.frame(id=1:300))

proj4string(x)
# [1] " +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
proj4string(p)
# [1] " +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

plot(x)
plot(p, col="red" , add=TRUE)

res <- over(p, x)
table(res$NAME_1) # count points
#               Abruzzo                Apulia            Basilicata
#                    11                    20                     9
#              Calabria              Campania        Emilia-Romagna
#                    16                     8                    25
# Friuli-Venezia Giulia                 Lazio               Liguria
#                     7                    14                     7
#             Lombardia                Marche                Molise
#                    22                     4                     3
#              Piemonte              Sardegna                Sicily
#                    35                    18                    21
#               Toscana   Trentino-Alto Adige                Umbria
#                    33                    15                     6
#         Valle d'Aosta                Veneto
#                     4                    22

Quiero dejar otra opción Puede lograr la tarea utilizando poly.counts()en el GISToolspaquete. Usando los datos de muestra de rcs, puede hacer lo siguiente. Si analizara la función, se daría cuenta de que la función está escrita como colSums(gContains(polys, pts, byid = TRUE)). Por lo tanto, puede usar gContains()en el rgeospaquete y colSums().

library(GISTools)

poly.counts(p, x) -> res
setNames(res, x@data$NAME_1)

O

colSums(gContains(x, p, byid = TRUE)) -> res
setNames(res, x@data$NAME_1)

Y el resultado es:

#              Abruzzo                Apulia            Basilicata 
#                   11                    20                     9 
#             Calabria              Campania        Emilia-Romagna 
#                   16                     8                    25 
#Friuli-Venezia Giulia                 Lazio               Liguria 
#                    7                    14                     7 
#            Lombardia                Marche                Molise 
#                   22                     4                     3 
#             Piemonte              Sardegna                Sicily 
#                   35                    18                    21 
#              Toscana   Trentino-Alto Adige                Umbria 
#                   33                    15                     6 
#        Valle d'Aosta                Veneto 
#                    4                    22 

Puede lograr lo mismo usando el sfpaquete. Verifique el código reproducible y comentado a continuación. El paquete sfse usa para manejar objetos espaciales como objetos de características simples. En esta respuesta, el paquete rasterse usa solo para descargar datos de polígonos de ejemplo y el paquete dplyrpara la transformación de datos al final.

# 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 poins 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 polygon + points
plot(polygon, graticule = st_crs(4326), key.pos = 1)
plot(points, pch = 19, col = "black", add = TRUE)

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

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

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

# View result
table(intersection$id_polygons) # using table

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

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