Gabungkan data titik spasial ke poligon di R

21

Saya mencoba melakukan penggabungan spasial antara data titik dan data poligon.

Saya memiliki data yang menunjukkan koordinat spasial dari suatu peristiwa dalam file csv saya A dan memiliki file lain, shapefile B, yang berisi batas-batas suatu area sebagai poligon.

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>

Saya ingin bergabung dengan data kejahatan A ke shapefile B saya untuk memetakan peristiwa kejahatan yang terjadi di daerah saya A. Sayangnya saya tidak dapat melakukan penggabungan atribut berdasarkan codekode dalam A mengacu pada unit yang berbeda dari kode di B.

Saya telah membaca sejumlah tutorial dan posting tetapi tidak dapat menemukan jawaban. Saya mencoba:

joined = over(A, B)

dan overlay, tetapi tidak mencapai apa yang saya inginkan.

Apakah ada cara untuk melakukan ini secara langsung atau akankah transformasi antara dari A ke format lain diperlukan?

Secara konseptual saya ingin memilih titik-titik A yang termasuk dalam codearea B (mirip dengan "bergabung berdasarkan lokasi spasial di ArcGIS").

Apakah seseorang memiliki masalah ini dan menyelesaikannya?

ben_aaron
sumber
Sudahkah Anda melihat point.in.polygon()dalam paket sp?
@ arvi1000 Saya sudah dan akan mencoba ini lagi. Pikiranku tentang point.in.polygonapakah ini akan melestarikan variabel monthdan crime_type. Apakah kamu tahu tentang itu?
ben_aaron
Saya sudah mencoba sedikit lebih banyak point.in.polydan akhirnya memilih titik-titik yang termasuk dalam poligon yang relevan. Terima kasih.
ben_aaron
Maka mungkin Anda harus menjawab pertanyaan Anda sendiri dengan solusi Anda. Ingat, jawaban yang baik adalah semua tentang situs ini.
SlowLearner

Jawaban:

8

Fungsi point.in.poly dalam paket spatialEco mengembalikan objek SpatialPointsDataFrame dari titik-titik yang memotong objek poligon sp dan secara opsional menambahkan atribut poligon.

Pertama mari kita tambahkan paket yang dibutuhkan dan buat beberapa contoh data.

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

Sekarang, mari kita melihat data dan memplotnya.

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

Akhirnya, kita bisa memotong titik dengan poligon. Hasilnya akan menjadi objek SpatialPointsDataFrame dengan, dalam hal ini, dua atribut tambahan (PIDS, y) yang terkandung dalam data poligon srdf.

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

Jika tidak ada kolom identifikasi unik dalam data poligon, Anda dapat dengan mudah menambahkannya.

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

Setelah titik dan poligon berpotongan, kami dapat mengumpulkan titik menggunakan ID poligon unik yang merupakan atribut dalam data poligon.

# 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
sumber
@ arvi1000, ya tapi sp :: point.in.polygon menghasilkan yang logis. SpatialEco: point.in.poly adalah pembungkus untuk over tetapi mengembalikan sp SpatialPointsDataFrame dan cara pintas beberapa langkah dalam menghubungkan atribut poligon, seperti raster: intersect tidak untuk rgeos :: gIntersect.
Jeffrey Evans
sp::point.in.polygonsebenarnya mengembalikan nilai numerik (0 = titik di luar, 1 = di dalam, 2 = di tepi, 3 = di titik). Bisa jadi hal yang tepat untuk beberapa keadaan. Berpikir itu membantu untuk dicatat di sini, karena ini adalah hasil google teratas untuk istilah terkait
arvi1000
27

over()dari paket spbisa sedikit membingungkan tetapi berfungsi dengan baik. Saya berasumsi Anda telah membuat spasial "A" dengan coordinates(A) <- ~longitude+latitude:

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

Alih-alih objek spasial titik, ini hanya memberi Anda bingkai data, dengan no yang sama. baris sebagai A, dan "kode" variabel tunggal dari setiap poligon berpotongan dari B.

# Add that data back to A:
A$bcode <- a.data$code
Simbamangu
sumber
Saya menemukan over()ada masalah dengan poin pada simpul poligon, meskipun saya pikir ini adalah solusi termudah yang saya temukan sejauh ini.
JMT2080AD
Masalah apa yang Anda miliki?
Simbamangu
Pengecualian. Saya perlu menjelajah lebih jauh. Saya akan memberi Anda beberapa data nanti hari ini dan kami dapat melihatnya bersama-sama jika Anda tertarik. Saya mungkin salah, tetapi saya cukup yakin ada beberapa kemunduran dalam algoritma yang perlu dijaga, setidaknya untuk data saya.
JMT2080AD
Sudahlah. Pasti ada sesuatu dengan data saya. Set eksperimental ini berfungsi dengan baik. r-fiddle.org/#/fiddle?id=m5sTjE4N&version=1
JMT2080AD
1
Ini adalah pendekatan yang jauh lebih mudah daripada jawaban yang diterima, dan tidak perlu menginstal paket tambahan selain rgdal.
Bryce Frank
0

Berikut ini adalah solusi seperti 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"))

Data populasi berasal dari: https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/parlemenaryconstituencymidyearpopulationestimates

Saya harus mengonversi file bentuk yang diunduh dari ke geoJson: https://geoportal.statistics.gov.uk/datasets/westminster-parliamentary-constituencies-december-2018-uk-bgc/data?page=1

Anda dapat melakukannya dengan:

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
sumber