Unir datos espaciales de puntos a polígonos en R

Question

Estoy intentando realizar una unión espacial entre datos de puntos y datos de polígonos.

Tengo datos que indican las coordenadas espaciales de un evento en mi archivo csv A y tengo otro archivo, shapefile B, que contiene los límites de un área como polígonos.

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>

Quiero unir los datos de delincuencia A a mi shapefile B para mapear los eventos de delincuencia que ocurren en mi área A. Lamentablemente, no puedo realizar una unión de atributos basada en code ya que el código en A se refiere a unidades diferentes que el código en B.

He leído varios tutoriales y posts pero no he podido encontrar una respuesta. Lo he intentado:

joined = over(A, B)

y overlay Pero no logré lo que quería.

¿Hay alguna manera de hacer esta unión directamente o sería necesaria una transformación intermedia de A a otro formato?

Conceptualmente quiero seleccionar aquellos puntos de A que caen en el code áreas de B (similar a la "unión basada en la ubicación espacial en ArcGIS").

¿Alguien ha tenido este problema y lo ha resuelto?

Answer 1

over() del paquete sp puede ser un poco confuso pero funciona bien. Asumo que ya has hecho la "A" espacial con coordinates(A) <- ~longitude+latitude :

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

En lugar de un objeto espacial de punto, esto simplemente le da un marco de datos, con el mismo número de filas que A, y una sola variable "código" de cada polígono de intersección de B.

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

Answer 2

La función point.in.poly del paquete spatialEco devuelve un objeto SpatialPointsDataFrame de los puntos que intersectan un objeto poligonal sp y opcionalmente añade los atributos del polígono.

Primero vamos a añadir los paquetes requeridos y crear algunos datos de ejemplo.

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)))

Ahora, vamos a echar un vistazo a los datos y a trazarlos.

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

Por último, podemos intersecar los puntos con los polígonos. El resultado será un objeto SpatialPointsDataFrame con, en este caso, dos atributos extra (PIDS, y) que estaban contenidos en los datos del polígono srdf.

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

Si no hay una columna de identificación única en los datos de los polígonos, podría añadir una fácilmente.

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

Una vez que tenemos los puntos y los polígonos intersectados, podemos agregar los puntos usando los ID's únicos de los polígonos que eran un atributo en los datos de los polígonos.

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

Answer 3

Aquí hay una solución similar a la de 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"))

Los datos de población proceden de: https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/parliamentaryconstituencymidyearpopulationestimates

Tuve que convertir los archivos shape descargados a geoJson: https://geoportal.statistics.gov.uk/datasets/westminster-parliamentary-constituencies-december-2018-uk-bgc/data?page=1

Puedes hacerlo:

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")