Unir datos de puntos espaciales a polígonos en R

21

Estoy tratando de 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 del crimen A a mi archivo de forma B para mapear los eventos del crimen que suceden en mi área A. Lamentablemente, no puedo realizar una unión de atributos basada en codeque el código en A se refiere a unidades diferentes al código en B.

He leído varios tutoriales y publicaciones, pero no pude encontrar una respuesta. Lo intenté:

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 las codeáreas de B (similar a "unirse según la ubicación espacial en ArcGIS").

¿Alguien tuvo este problema y lo resolvió?

ben_aaron
fuente
¿Has mirado point.in.polygon()en el paquete sp?
@ arvi1000 Tengo e intentaré esto de nuevo. Pensé en point.in.polygonsi esto preservaría las variables monthy crime_type. ¿Sabe usted acerca de eso?
ben_aaron
He intentado un poco más point.in.polyy finalmente he seleccionado los puntos que se encuentran en los polígonos relevantes. Gracias.
ben_aaron
Entonces quizás deberías responder tu propia pregunta con tu solución. Recuerde, las buenas respuestas son de lo que se trata este sitio.
SlowLearner

Respuestas:

8

La función point.in.poly en el paquete spatialEco devuelve un objeto SpatialPointsDataFrame de los puntos que se cruzan con un objeto de polígono sp y opcionalmente agrega los atributos de polígono.

Primero agreguemos los paquetes requeridos y creemos 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, echemos un vistazo rápido a los datos y tracémoslos.

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

Finalmente, podemos intersecar los puntos con los polígonos. Los resultados serán un objeto SpatialPointsDataFrame con, en este caso, dos atributos adicionales (PIDS, y) 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 del polígono, puede agregar fácilmente una.

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

Una vez que tenemos los puntos y los polígonos intersectados, podemos agregar los puntos utilizando las ID de polígono únicas que eran un atributo en los datos del polígono.

# 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)
Jeffrey Evans
fuente
@ arvi1000, sí, pero sp :: point.in.polygon produce una lógica. SpatialEco: point.in.poly es un contenedor para over pero devuelve un sp SpatialPointsDataFrame y ataja algunos pasos para relacionar los atributos del polígono, de forma muy similar a raster: intersect para rgeos :: gIntersect.
Jeffrey Evans
sp::point.in.polygonen realidad devuelve un valor numérico (0 = el punto está afuera, 1 = adentro, 2 = en el borde, 3 = en el vértice). Podría ser lo correcto para algunas circunstancias. Pensé que era útil tener en cuenta aquí, ya que este es un resultado superior de Google para términos relacionados
arvi1000
27

over()del paquete sppuede ser un poco confuso pero funciona bien. Supongo que ya has hecho "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 puntual, esto simplemente le da un marco de datos, con el mismo no. filas como A y un único "código" variable de cada polígono de intersección de B.

# Add that data back to A:
A$bcode <- a.data$code
Simbamangu
fuente
Descubrí over()que tengo problemas con los puntos en los vértices de los polígonos, aunque creo que esta es la solución más fácil que he encontrado hasta ahora.
JMT2080AD
¿Qué problemas has tenido?
Simbamangu
Exclusión. Necesito explorarlo más a fondo. Te enviaré algunos datos más tarde hoy y podemos verlos juntos si estás interesado. Podría estar equivocado, pero estoy bastante seguro de que hay algunas degeneraciones en el algoritmo que deben ser atendidas, al menos para mis datos.
JMT2080AD
No importa. Debe ser algo con mis datos. Este conjunto experimental funciona bien. r-fiddle.org/#/fiddle?id=m5sTjE4N&version=1
JMT2080AD
1
Este es un enfoque mucho más directo que la respuesta aceptada, y no requiere instalar paquetes adicionales que no sean rgdal.
Bryce Frank
0

Aquí hay una solución como 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 provienen de: https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/parliamentaryconstituencymidyearpopulationestimates

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

Puedes hacerlo de la siguiente manera:

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