Joindre des données de points spatiaux à des polygones dans R

J'essaie d'effectuer une jointure spatiale entre les données de points et les données de polygones.

J'ai des données qui indiquent les coordonnées spatiales d'un événement dans mon fichier csv A et j'ai un autre fichier, le fichier de formes B, qui contient les limites d'une zone sous forme de polygones.

head(A)
  month   longitude latitude lsoa_code                   crime_type
1 2014-09 -1.550626 53.59740 E01007359        Anti-social behaviour
2 2014-09 -1.550626 53.59740 E01007359                 Public order
3 2014-09 -1.865236 53.93678 E01010646        Anti-social behaviour

head(B@data)
  code      name                                  altname
0 E05004934 Longfield, New Barn and Southfleet    <NA>
1 E05000448                   Lewisham Central    <NA>
2 E05003149                            Hawcoat    <NA>

Je souhaite joindre les données de criminalité A à mon fichier de formes B pour cartographier les événements de criminalité qui se produisent dans ma région A. Malheureusement, je ne peux pas effectuer de jointure d'attribut codecar le code en A fait référence à des unités différentes de celles du code en B.

J'ai lu un certain nombre de didacticiels et d'articles, mais je n'ai pas trouvé de réponse. J'ai essayé:

joined = over(A, B)

et overlay , mais n'a pas accompli ce que je voulais.

Existe-t-il un moyen de faire cette jointure directement ou une transformation intermédiaire de A vers un autre format serait-elle nécessaire?

Conceptuellement, je veux sélectionner les points de A qui tombent dans le code zones de B (similaires à «joindre en fonction de l'emplacement spatial dans ArcGIS»).

Quelqu'un a-t-il eu ce problème et l'a-t-il résolu?

La fonction point.in.poly du package spatialEco renvoie un objet SpatialPointsDataFrame des points qui coupent un objet polygone sp et ajoute éventuellement les attributs de polygone.

Ajoutons d'abord les packages requis et créons des exemples de données.

require(spatialEco)
require(sp)
data(meuse)
coordinates(meuse) = ~x+y
sr1=Polygons(list(Polygon(cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409,
  180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676,
  332618, 332413, 332349)))),'1')
sr2=Polygons(list(Polygon(cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437,
  179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
  331133, 331623, 332152, 332357, 332373)))),'2')
sr3=Polygons(list(Polygon(cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875,
  179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
  c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
  329783, 329665, 329720, 329933, 330478, 331062, 331086)))),'3')
sr4=Polygons(list(Polygon(cbind(c(180304, 180403,179632,179420,180304),
  c(332791, 333204, 333635, 333058, 332791)))),'4')
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(row.names=c('1','2','3','4'), PIDS=1:4, y=runif(4)))

Maintenant, jetons un coup d'œil aux données et traçons-les.

head(srdf@data)  # polygons
head(meuse@data) # points
plot(srdf)
points(meuse, pch=20)

Enfin, nous pouvons intersecter les points avec les polygones. Le résultat sera un objet SpatialPointsDataFrame avec, dans ce cas, deux attributs supplémentaires (PIDS, y) qui étaient contenus dans les données du polygone srdf.

  pts.poly <- point.in.poly(meuse, srdf)
    head(pts.poly@data)

S'il n'y a pas de colonne d'identification unique dans les données du polygone, vous pouvez facilement en ajouter une.

srdf@data$poly.ids <- 1:nrow(srdf) 

Une fois que nous avons intersecté les points et les polygones, nous pouvons agréger les points en utilisant les ID de polygone uniques qui étaient un attribut dans les données de polygone.

# Number of points in each polygon
tapply(pts.poly@data$lead, pts.poly@data$PIDS, FUN=length)

# Mean lead in each polygon
tapply(pts.poly@data$lead, pts.poly@data$PIDS, FUN=mean)

over()à partir du paquet sppeut être un peu déroutant mais fonctionne bien. Je suppose que vous avez déjà créé l'espace "A" avec coordinates(A) <- ~longitude+latitude:

# Overlay points and extract just the code column: 
a.data <- over(A, B[,"code"])

Au lieu d'un objet géographique ponctuel, cela vous donne simplement un bloc de données, avec le même non. lignes comme A, et une seule variable "code" de chaque polygone qui se croisent de B.

# Add that data back to A:
A$bcode <- a.data$code

Voici une solution similaire à dplyr:

library(spdplyr)

ukcounties <- geojsonio::geojson_read("data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/uk_country.geojson",
                                      what = "sp")
pop <- read_excel("data/SAPE20DT7-mid-2017-parlicon-syoa-estimates-unformatted.xls",sheet = "data")
pop <- janitor::clean_names(pop)

ukcounties_pop <- ukcounties %>% inner_join(pop, by = c("pcon18nm" = "pcon11nm"))

Les données démographiques proviennent de: https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/parliamentaryconstituencymidyearpopulationestimates

J'ai dû convertir les fichiers de forme téléchargés depuis en geoJson: https://geoportal.statistics.gov.uk/datasets/westminster-parliamentary-constituencies-december-2018-uk-bgc/data?page=1

Vous pouvez le faire en:

uk_constituencies <- readOGR("data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC.shp")
uk_constituencies # this is in tmerc format. we need to convert it to WGS84 required by geoJson format.

# First Convert to Longitude / Latitude with WGS84 Coordinate System
wgs84 = '+proj=longlat +datum=WGS84'
uk_constituencies_trans <- spTransform(uk_constituencies, CRS(wgs84))

# Convert from Spatial Dataframe to GeoJSON
uk_constituencies_json <- geojson_json(uk_constituencies_trans)

# Save as GeoJSON file on the file system.
geojson_write(uk_constituencies_json, file = "data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/uk_country.geojson")

#read back in:
ukcounties <- geojsonio::geojson_read("data/Westminster_Parliamentary_Constituencies_December_2018_UK_BGC/uk_country.geojson",
                                      what = "sp")