Połącz dane punktu przestrzennego z wielokątami w R.

21

Próbuję wykonać połączenie przestrzenne między danymi punktów i danymi wielokąta.

Mam dane wskazujące współrzędne przestrzenne zdarzenia w moim pliku csv A i mam inny plik, plik kształtu B, który zawiera granice obszaru jako wielokąty.

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>

Chcę dołączyć dane przestępstwa A do mojego pliku kształtu B, aby zmapować zdarzenia przestępcze, które mają miejsce w mojej okolicy A. Niestety nie mogę wykonać łączenia atrybutów, codeponieważ kod w A odnosi się do innych jednostek niż kod w B.

Przeczytałem wiele samouczków i postów, ale nie mogłem znaleźć odpowiedzi. Próbowałem:

joined = over(A, B)

i overlay, ale nie osiągnąłem tego, co chciałem.

Czy istnieje sposób na wykonanie tego połączenia bezpośrednio, czy też potrzebna byłaby przejściowa transformacja z A do innego formatu?

Koncepcyjnie chcę wybrać te punkty A, które mieszczą się w codeobszarach B (podobne do „łączenia na podstawie położenia przestrzennego w ArcGIS”).

Czy ktoś miał ten problem i go rozwiązał?

ben_aaron
źródło
Czy oglądałeś point.in.polygon()w paczce sp?
@ arvi1000 Mam i spróbuję ponownie. Myślałem o point.in.polygontym, czy to pozwoli zachować zmienne monthi crime_type. Czy wiesz o tym?
ben_aaron
Próbowałem jeszcze trochę point.in.polyi ostatecznie wybrałem te punkty, które pasują do odpowiednich wielokątów. Dzięki.
ben_aaron
Być może powinieneś odpowiedzieć na swoje pytanie z rozwiązaniem. Pamiętaj, że o to właśnie chodzi w tej witrynie.
SlowLearner

Odpowiedzi:

8

Funkcja point.in.poly w pakiecie spatialEco zwraca obiekt SpatialPointsDataFrame punktów, które przecinają obiekt wielokąta sp i opcjonalnie dodaje atrybuty wielokąta.

Najpierw dodajmy wymagane pakiety i utwórz przykładowe dane.

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

Teraz rzućmy okiem na dane i wykreślmy je.

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

Na koniec możemy przecinać punkty wielokątami. Wynikiem będzie obiekt SpatialPointsDataFrame z dwoma dodatkowymi atrybutami (PIDS, y) zawartymi w danych wielokąta srdf.

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

Jeśli w danych wielokąta nie ma unikalnej kolumny identyfikacyjnej, można ją łatwo dodać.

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

Po przecięciu punktów i wielokątów możemy agregować punkty przy użyciu unikalnych identyfikatorów wielokątów, które były atrybutem w danych wielokąta.

# 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
źródło
@ arvi1000, tak, ale sp :: point.in.polygon tworzy logikę. SpatialEco: point.in.poly jest opakowaniem dla over, ale zwraca spatial SpektrPointsDataFrame i skraca niektóre kroki w powiązaniu atrybutów wielokąta, podobnie jak raster: intersect robi dla rgeos :: gIntersect.
Jeffrey Evans,
sp::point.in.polygonfaktycznie zwraca wartość liczbową (0 = punkt jest na zewnątrz, 1 = wewnątrz, 2 = na krawędzi, 3 = na wierzchołku). Może być właściwa w niektórych okolicznościach. Pomyślałem, że warto tutaj odnotować, ponieważ jest to najlepszy wynik w Google dla pokrewnych terminów
arvi1000 14.04.16
27

over()z pakietu spmoże być trochę mylące, ale działa dobrze. Zakładam, że masz już przestrzenne „A” z coordinates(A) <- ~longitude+latitude:

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

Zamiast punktowego obiektu przestrzennego daje to po prostu ramkę danych o tym samym numerze. wiersze jako A i pojedynczy zmienny „kod” z każdego przecinającego się wielokąta z B.

# Add that data back to A:
A$bcode <- a.data$code
Simbamangu
źródło
Odkryłem, że over()mam problemy z punktami na wierzchołkach wielokątów, chociaż myślę, że jest to najłatwiejsze rozwiązanie, jakie do tej pory znalazłem.
JMT2080AD
Jakie miałeś problemy?
Simbamangu,
Wykluczenie. Muszę to zbadać dalej. Prześlemy ci dzisiaj trochę danych i możemy je obejrzeć razem, jeśli jesteś zainteresowany. Mogę się mylić, ale jestem prawie pewien, że w algorytmie są pewne zwyrodnienia, które należy załatwić, przynajmniej w przypadku moich danych.
JMT2080AD
Nieważne. To musi być coś z moimi danymi. Ten zestaw eksperymentalny działa dobrze. r-fiddle.org/#/fiddle?id=m5sTjE4N&version=1
JMT2080AD
1
Jest to znacznie prostsze podejście niż zaakceptowana odpowiedź i nie wymaga instalowania dodatkowych pakietów innych niż rgdal.
Bryce Frank
0

Oto rozwiązanie podobne do 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"))

Dane dotyczące populacji pochodzą z: https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/parliamentaryconstituencymidyearpopulationestimates

Musiałem przekonwertować pliki kształtów pobrane z geoJson: https://geoportal.statistics.gov.uk/datasets/westminster-parliamentary-constituencies- grudnia 2018-uk-bgc/data?page=1

Możesz to zrobić przez:

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
źródło