Verbinden Sie räumliche Punktdaten mit Polygonen in R

21

Ich versuche, eine räumliche Verknüpfung zwischen Punktdaten und Polygondaten durchzuführen.

Ich habe Daten, die die räumlichen Koordinaten eines Ereignisses in meiner CSV-Datei A angeben, und eine andere Datei, Shapefile B, die die Grenzen eines Bereichs als Polygone enthält.

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>

Ich möchte die Verbrechensdaten A mit meinem Shapefile B verknüpfen, um die in meinem Bereich A aufgetretenen Verbrechensereignisse abzubilden. Leider kann ich keine Attributverknüpfung durchführen, codeda sich der Code in A auf andere Einheiten als der Code in B bezieht.

Ich habe eine Reihe von Tutorials und Beiträgen gelesen, konnte aber keine Antwort finden. Ich habe es versucht:

joined = over(A, B)

und overlay, aber habe nicht erreicht, was ich wollte.

Gibt es eine Möglichkeit, diesen Join direkt auszuführen, oder wäre eine Zwischentransformation von A in ein anderes Format erforderlich?

Konzeptionell möchte ich die Punkte von A auswählen, die in die codeBereiche von B fallen (ähnlich wie "Verbinden basierend auf der räumlichen Position in ArcGIS").

Hatte jemand dieses Problem und löste es?

ben_aaron
quelle
Hast du point.in.polygon()im Paket angeschaut sp?
@ arvi1000 Ich habe und werde das nochmal versuchen. Ich dachte darüber point.in.polygonnach, ob dies die Variablen monthund erhalten würde crime_type. Wissen Sie darüber Bescheid?
ben_aaron
Ich habe ein bisschen mehr mit versucht point.in.polyund schließlich diejenigen Punkte ausgewählt, die in die relevanten Polygone fallen. Vielen Dank.
ben_aaron
Dann sollten Sie vielleicht Ihre eigene Frage mit Ihrer Lösung beantworten. Denken Sie daran, auf dieser Website dreht sich alles um gute Antworten.
SlowLearner

Antworten:

8

Die point.in.poly-Funktion im spaciousEco-Paket gibt ein SpatialPointsDataFrame-Objekt der Punkte zurück, die ein sp-Polygonobjekt schneiden, und fügt optional die Polygonattribute hinzu.

Fügen Sie zunächst die erforderlichen Pakete hinzu und erstellen Sie einige Beispieldaten.

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

Lassen Sie uns nun einen kurzen Blick auf die Daten werfen und sie plotten.

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

Schließlich können wir die Punkte mit den Polygonen schneiden. Das Ergebnis ist ein SpatialPointsDataFrame-Objekt mit in diesem Fall zwei zusätzlichen Attributen (PIDS, y), die in den srdf-Polygondaten enthalten waren.

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

Wenn die Polygondaten keine eindeutige Identifizierungsspalte enthalten, können Sie problemlos eine hinzufügen.

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

Sobald wir die Schnittpunkte und Polygone haben, können wir die Punkte unter Verwendung der eindeutigen Polygon-IDs aggregieren, die ein Attribut in den Polygondaten waren.

# 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
quelle
@ arvi1000, ja, aber sp :: point.in.polygon erzeugt eine logische. Die Datei "spatialEco: point.in.poly" ist ein Wrapper für over, gibt jedoch einen sp-Wert für "SpatialPointsDataFrame" zurück und führt einige Schritte zum Verknüpfen der Polygonattribute aus, ähnlich wie bei "raster: intersect" für "rgeos :: gIntersect".
Jeffrey Evans
sp::point.in.polygonGibt tatsächlich einen numerischen Wert zurück (0 = Punkt liegt außerhalb, 1 = innerhalb, 2 = an der Kante, 3 = am Scheitelpunkt). Könnte unter bestimmten Umständen das Richtige sein. Dachte, es war hilfreich, hier zu beachten, da dies ein Top-Google-Ergebnis für verwandte Begriffe ist
arvi1000
27

over()from package spkann etwas verwirrend sein, funktioniert aber gut. Ich gehe davon aus, dass Sie "A" bereits räumlich gemacht haben mit coordinates(A) <- ~longitude+latitude:

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

Anstelle eines Punkt-Raumobjekts erhalten Sie einfach einen Datenrahmen mit derselben Nummer. Zeilen als A und eine einzelne Variable "Code" von jedem sich überschneidenden Polygon von B.

# Add that data back to A:
A$bcode <- a.data$code
Simbamangu
quelle
Ich habe over()Probleme mit Punkten an den Eckpunkten der Polygone festgestellt, obwohl ich denke, dass dies die einfachste Lösung ist, die ich bisher gefunden habe.
JMT2080AD
Welche Probleme hatten Sie?
Simbamangu
Ausschluss. Ich muss es weiter erforschen. Ich werde Ihnen heute noch einige Daten zusenden und wir können sie bei Interesse gemeinsam prüfen. Ich könnte mich irren, aber ich bin mir ziemlich sicher, dass es einige Entartungen im Algorithmus gibt, auf die zumindest für meine Daten geachtet werden muss.
JMT2080AD
Keine Ursache. Es muss etwas mit meinen Daten sein. Dieser experimentelle Satz funktioniert gut. r-fiddle.org/#/fiddle?id=m5sTjE4N&version=1
JMT2080AD
1
Dies ist eine viel einfachere Vorgehensweise als die akzeptierte Antwort und erfordert keine Installation zusätzlicher Pakete außer rgdal.
Bryce Frank
0

Hier ist eine Dplyr-ähnliche Lösung:

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

Die Bevölkerungsdaten stammen von: https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/parliamentaryconstituencymidyearpopulationestimates

Ich musste die Formdateien, die von heruntergeladen wurden, in geoJson konvertieren: https://geoportal.statistics.gov.uk/datasets/westminster-parliamentary-constituencies-december-2018-uk-bgc/data?page=1

Sie können dies tun, indem Sie:

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
quelle