Zliczanie punktów w wielokącie za pomocą R?

Mam dwie klasy dzielące ten sam CRS (Latitit i Longitude):

1. bolognaQuartieriMap: SpatialPolygonDataFramezawierające dane gmin miejskich.
2. crashPoints: SpatialPointsDataFramezawierający dane o wypadkach.

Są dobrze wykreślone przy użyciu:

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

Potrzebuję uzyskać liczbę punktów ( crashPoints) w każdym wielokącie, które stanowią bolognaQuartieriMap. Zaproponowano mi użycie, over()ale mi się nie udało.

Ponieważ nie podałeś odtwarzalnego przykładu ani komunikatu o błędzie, sprawdź, czy ten fragment kodu pozwala rozpocząć:

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

Chcę zostawić inną opcję. Można to osiągnąć za pomocą zadania poly.counts()w GISToolspakiecie. Używając przykładowych danych przez rcs, możesz wykonać następujące czynności. Jeśli spojrzysz na funkcję, zdasz sobie sprawę, że funkcja jest zapisana jako colSums(gContains(polys, pts, byid = TRUE)). Możesz więc po prostu użyć gContains()w rgeospakiecie i colSums().

library(GISTools)

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

Lub

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

Rezultat to:

#              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 

Możesz to samo osiągnąć za pomocą sfpakietu. Sprawdź odtwarzalny i skomentowany kod poniżej. Pakiet sfsłuży do obsługi obiektów przestrzennych jako prostych obiektów obiektów. W tej odpowiedzi pakiet rastersłuży tylko do pobrania przykładowych danych wielokąta, a pakiet dplyrdo transformacji danych na końcu.

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