Unir dados de pontos espaciais a polígonos em R

21

Estou tentando executar uma junção espacial entre dados de ponto e dados de polígono.

Eu tenho dados que indicam as coordenadas espaciais de um evento no meu arquivo csv A e outro arquivo, shapefile B, que contém os limites de uma á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>

Quero associar os dados de crime A ao meu shapefile B para mapear os eventos de crime que ocorrem na minha área A. Infelizmente, não consigo executar uma associação de atributo com base codeno código em A se referir a unidades diferentes do código em B.

Eu li vários tutoriais e postagens, mas não consegui encontrar uma resposta. Eu tentei:

joined = over(A, B)

e overlay, mas não conseguiu o que eu queria.

Existe uma maneira de fazer essa junção diretamente ou seria necessária uma transformação intermediária de A para outro formato?

Conceitualmente, quero selecionar os pontos de A que se enquadram nas codeáreas de B (semelhante a "unir com base na localização espacial no ArcGIS").

Alguém teve esse problema e resolveu?

ben_aaron
fonte
Você já viu no point.in.polygon()pacote sp?
@ arvi1000 Tenho e tentarei novamente. Meu pensamento point.in.polygonera se isso preservaria as variáveis monthe crime_type. Você sabe sobre aquilo?
ben_aaron
Eu tentei um pouco mais point.in.polye finalmente selecionei os pontos que se enquadram nos polígonos relevantes. Obrigado.
ben_aaron
Então talvez você deva responder sua própria pergunta com sua solução. Lembre-se de que são boas respostas sobre este site.
SlowLearner

Respostas:

8

A função point.in.poly no pacote spatialEco retorna um objeto SpatialPointsDataFrame dos pontos que cruzam um objeto de polígono sp e, opcionalmente, adiciona os atributos de polígono.

Primeiro vamos adicionar os pacotes requeridos e criar alguns dados de exemplo.

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

Agora, vamos dar uma olhada rápida nos dados e plotá-los.

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

Finalmente, podemos cruzar os pontos com os polígonos. Os resultados serão um objeto SpatialPointsDataFrame com, nesse caso, dois atributos extras (PIDS, y) que estavam contidos nos dados do polígono srdf.

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

Se não houver uma coluna de identificação exclusiva nos dados do polígono, você poderá facilmente adicionar uma.

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

Depois que os pontos e polígonos estiverem interceptados, podemos agregar os pontos usando os IDs de polígonos exclusivos que eram um atributo nos dados do 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
fonte
@ arvi1000, sim, mas sp :: point.in.polygon produz uma lógica. O spatialEco: point.in.poly é um wrapper para over, mas retorna um sp SpatialPointsDataFrame e atalhos para algumas etapas na relação dos atributos de polígono, como raster: intersect faz para rgeos :: gIntersect.
Jeffrey Evans
sp::point.in.polygonna verdade, retorna um valor numérico (0 = ponto está fora, 1 = dentro, 2 = na aresta, 3 = no vértice). Pode ser a coisa certa para algumas circunstâncias. Pensei que era útil notar aqui, uma vez que este é um resultado de topo no Google por termos relacionados
arvi1000
27

over()do pacote sppode ser um pouco confuso, mas funciona bem. Suponho que você já tenha feito o "A" espacial com coordinates(A) <- ~longitude+latitude:

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

Em vez de um objeto espacial pontual, isso simplesmente fornece um quadro de dados, com o mesmo não. linhas como A e uma única variável "código" de cada polígono que se cruza de B.

# Add that data back to A:
A$bcode <- a.data$code
Simbamangu
fonte
Descobri over()que há problemas com pontos nos vértices dos polígonos, embora eu ache que essa seja a solução mais fácil que encontrei até agora.
JMT2080AD
Que problemas você já teve?
Simbamangu
Exclusão. Eu preciso explorar mais. Ainda hoje, apresentamos alguns dados e podemos analisá-los juntos, se você estiver interessado. Posso estar errado, mas tenho certeza de que existem algumas degenerescências no algoritmo que precisam ser resolvidas, pelo menos para os meus dados.
JMT2080AD
Deixa pra lá. Deve ser algo com meus dados. Este conjunto experimental funciona bem. r-fiddle.org/#/fiddle?id=m5sTjE4N&version=1
JMT2080AD
1
Essa é uma abordagem muito mais direta do que a resposta aceita e não requer a instalação de pacotes adicionais além do rgdal.
Bryce Frank
0

Aqui está uma solução semelhante ao 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"))

Os dados da população são provenientes de: https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/parl Parliamentaryconstituencymidyearpopulationestimates

Eu tive que converter os arquivos de forma baixados de para geoJson: https://geoportal.statistics.gov.uk/datasets/westminster-parliamentary-constituencies-december-2018-uk-bgc/data?page=1

Você pode fazer isso:

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
fonte