Соединить данные пространственной точки с полигонами в R

21

Я пытаюсь выполнить пространственное соединение между данными точек и данными многоугольников.

У меня есть данные, которые указывают пространственные координаты события в моем CSV-файле A, и у меня есть другой файл, шейп-файл B, который содержит границы области в виде полигонов.

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>

Я хочу соединить данные преступления A с моим шейп-файлом B, чтобы отобразить события преступления, которые происходят в моей области A. К сожалению, я не могу выполнить соединение атрибутов на основе, codeпоскольку код в A ссылается на другие единицы, чем код в B.

Я прочитал несколько уроков и постов, но не смог найти ответ. Я старался:

joined = over(A, B)

и overlay, но не выполнить то, что я хотел.

Есть ли способ сделать это соединение напрямую или потребуется промежуточное преобразование из А в другой формат?

Концептуально я хочу выбрать те точки A, которые попадают в codeобласти B (аналогично «объединению на основе пространственного расположения в ArcGIS»).

Был ли у кого-то этот вопрос и решил его?

ben_aaron
источник
Вы смотрели point.in.polygon()в упаковке sp?
@ arvi1000 У меня есть и попробую это снова. Я думал о point.in.polygonтом, сохранит ли это переменные monthи crime_type. Ты знаешь об этом?
ben_aaron
Я попробовал немного больше point.in.polyи наконец выбрал те точки, которые попадают в соответствующие многоугольники. Благодарю.
ben_aaron
Тогда, возможно, вы должны ответить на свой вопрос с вашим решением. Помните, что хорошие ответы - это то, о чем этот сайт.
SlowLearner

Ответы:

8

Функция point.in.poly в пакетеatialEco возвращает объект SpatialPointsDataFrame точек, которые пересекают объект sp polygon, и дополнительно добавляет атрибуты polygon.

Сначала давайте добавим требуемые пакеты и создадим пример данных.

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

Теперь давайте взглянем на данные и нанесем их на график.

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

Наконец, мы можем пересекать точки с полигонами. Результатом будет объект SpatialPointsDataFrame, в данном случае два дополнительных атрибута (PIDS, y), которые содержались в данных многоугольника srdf.

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

Если в данных многоугольника нет уникальной идентификационной колонки, ее можно легко добавить.

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

После пересечения точек и многоугольников мы можем объединить точки, используя уникальные идентификаторы многоугольников, которые были атрибутом в данных многоугольника.

# 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)
Джеффри Эванс
источник
@ arvi1000, да, но sp :: point.in.polygon производит логическое. ПространственноеEco: point.in.poly является оболочкой для over, но возвращает sp SpatialPointsDataFrame и сокращает некоторые шаги в связывании атрибутов полигона, так же как и растр: intersect для rgeos :: gIntersect.
Джеффри Эванс
sp::point.in.polygonфактически возвращает числовое значение (0 = точка снаружи, 1 = внутри, 2 = по краю, 3 = по вершине). Может быть правильным для некоторых обстоятельств. Думаю, что это было полезно отметить здесь, так как это лучший результат Google для связанных терминов
arvi1000
27

over()из пакета spможет быть немного запутанным, но работает хорошо. Я предполагаю, что вы уже сделали "A" пространственным с coordinates(A) <- ~longitude+latitude:

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

Вместо точечного пространственного объекта это просто дает вам фрейм данных с тем же самым номером. строки как A, и одна переменная "код" из каждого пересекающегося многоугольника из B.

# Add that data back to A:
A$bcode <- a.data$code
Simbamangu
источник
У меня over()возникли проблемы с точками в вершинах многоугольников, хотя я думаю, что это самое простое решение, которое я нашел до сих пор.
JMT2080AD
Какие проблемы у вас были?
Симбамангу
Исключение. Мне нужно изучить это дальше. Я отправлю вам некоторые данные позже сегодня, и мы можем посмотреть их вместе, если вам интересно. Возможно, я ошибаюсь, но я уверен, что в алгоритме есть некоторые вырождения, о которых нужно позаботиться, по крайней мере, для моих данных.
JMT2080AD
Ничего. Это должно быть что-то с моими данными. Этот экспериментальный набор работает отлично. r-fiddle.org/#/fiddle?id=m5sTjE4N&version=1
JMT2080AD
1
Это гораздо более простой подход, чем принятый ответ, и не требует установки дополнительных пакетов, кроме rgdal.
Брайс Франк
0

Вот решение, похожее на 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"))

Данные о населении получены по адресу : https://www.ons.gov.uk/peoplepopulationandcommunity/populationandmigration/populationestimates/datasets/par Parliamentaryconstituencymidyearpopulationestimates.

Мне пришлось преобразовать файлы форм, загруженные из, в geoJson: https://geoportal.statistics.gov.uk/datasets/westminster-par Parliamentary-constituencies-de December-2018-uk-bgc/data?page=1

Вы можете сделать это:

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")
Шери
источник