Compter le nombre de points dans un polygone en utilisant R?

J'ai deux classes partageant le même CRS (Latitute et Longitude):

1. bolognaQuartieriMap: un SpatialPolygonDataFramecontenant des données d'un arrondissement de la ville.
2. crashPoints: un SpatialPointsDataFramecontenant des données d'accidents.

Ils sont bien tracés en utilisant:

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

Ce dont j'ai besoin, c'est d'obtenir le nombre de points ( crashPoints) dans chaque polygone qui le constitue bolognaQuartieriMap. On m'a suggéré d'utiliser over()mais je n'ai pas réussi.

Étant donné que vous n'avez pas fourni d'exemple reproductible ni de message d'erreur, voyez si cet extrait de code vous permet de démarrer:

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

Je veux laisser une autre option. Vous pouvez réaliser la tâche en utilisant poly.counts()dans le GISToolspackage. En utilisant les exemples de données par rcs, vous pouvez effectuer les opérations suivantes. Si vous regardiez la fonction, vous vous rendriez compte que la fonction est écrite comme colSums(gContains(polys, pts, byid = TRUE)). Donc, vous pouvez simplement utiliser gContains()dans le rgeospackage et colSums().

library(GISTools)

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

Ou

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

Et le résultat est:

#              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 

Vous pouvez obtenir le même résultat en utilisant le sfpackage. Vérifiez le code reproductible et commenté ci-dessous. Le package sfest utilisé pour gérer les objets spatiaux comme de simples objets d'entités. Dans cette réponse, le package rasterest utilisé uniquement pour télécharger des exemples de données de polygone et le package dplyrpour la transformation des données à la fin.

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