Точки в нескольких полигонах, используя R

Question

В настоящее время у меня есть два data.frames, один из полигонов (poly.x, poly.y, enum) и одна из точек (pt.x, pt.y), где enum - это идентификатор многоугольника. Я пытаюсь определить, какие точки принадлежат каким полигонам, поэтому я получаю data.frame (pt.x, pt.y, enum).

Моя первая попытка использует point.in.polygon из пакета sp и функции lapply, чтобы найти, к какому полигону относится точка. Хотя мой код работает, на больших наборах данных требуется long .

Моя вторая попытка использует over также из пакета sp, собранного воедино из вопросов об обмене стеками ГИС. Хотя он намного быстрее, я не могу получить правильный вывод из over, поскольку это кадр данных 1 с и NA с.

Ниже я включил упрощенный рабочий пример (npoly можно изменить для проверки скорости различных методов), а также мою рабочую попытку с использованием sp::point.in.polygon и бессмысленный вывод из моей попытки sp::over. Я не смущаюсь, какой метод я использую, пока он быстрый.

Любая помощь будет принята с благодарностью!

#-------------------------------------------
# Libraries
library(ggplot2)  # sample plots
library(dplyr)    # bind_rows(), etc
library(sp)       # spatial data
# Sample data
npoly = 100
#  polygons
localpolydf <- data.frame(
  x = rep(c(0, 1, 1, 0), npoly) + rep(0:(npoly-1), each = 4),
  y = rep(c(0, 0, 1, 1), npoly),
  enum = rep(1:npoly, each = 4))
#  points
offsetdf <- data.frame(
  x = seq(min(localpolydf$x) - 0.5, max(localpolydf$x) + 0.5, by = 0.5),
  y = runif(npoly*2 + 3, 0, 1))
# Sample plot
ggplot() + 
  geom_polygon(aes(x, y, group = enum), 
               localpolydf, fill = NA, colour = "black") +
  geom_point(aes(x, y), offsetdf)

#-------------------------------------------
# Dplyr and lapply solution for point.in.polygon
ptm <- proc.time() # Start timer
#  create lists
offsetlist <- split(offsetdf, rownames(offsetdf))
polygonlist <- split(localpolydf, localpolydf$enum)
# lapply over each pt in offsetlist
pts <- lapply(offsetlist, function(pt) {
  # lapply over each polygon in polygonlist
  ptpoly <- lapply(polygonlist, function(poly) {
    data.frame(
      enum = poly$enum[1],
      ptin = point.in.polygon(pt[1,1], pt[1,2], poly$x, poly$y))
  })
  ptpoly <- bind_rows(ptpoly) %>% filter(ptin != 0)
  if (nrow(ptpoly) == 0) return(data.frame(x = pt$x, y = pt$y, enum = NA, ptin = NA))
  ptpoly$x = pt$x
  ptpoly$y = pt$y
  return(ptpoly[c("x", "y", "enum", "ptin")])
})
pts_apply <- bind_rows(pts)
proc.time() - ptm # end timer

#-------------------------------------------
# Attempted sp solution for over
ptm <- proc.time() # Start timer
# Split the dataframe into a list based on enum and then remove enum from df in the list
polygonlist <- split(localpolydf, localpolydf$enum)
polygonlist <- lapply(polygonlist, function(x) x[,c("x", "y")])
# Convert the list to Polygon, then create a Polygons object
polygonsp <- sapply(polygonlist, Polygon)
polygonsp <- Polygons(polygonsp, ID = 1)
polygonsp <- SpatialPolygons(list(polygonsp))
plot(polygonsp)
# Convert points to coordinates
offsetps <- offsetdf
coordinates(offsetps) <- ~x+y
points(offsetps$x, offsetps$y)
# Determine polygons points are in
pts_sp <- over(offsetps, polygonsp)
proc.time() - ptm # end timer

#===========================================
# Output
#  Apply: point.in.polygon
> head(pts_apply)
   x   y         enum ptin
1 -0.5 0.2218138   NA   NA
2  4.0 0.9785541    4    2
3  4.0 0.9785541    5    2
4 49.0 0.3971479   49    2
5 49.0 0.3971479   50    2
6 49.5 0.1177206   50    1
user  system elapsed 
4.434   0.002   4.435 
# SP: over
> head(pts_sp)
1  2  3  4  5  6 
NA  1  1 NA  1 NA 
user  system elapsed 
0.048   0.000   0.047 

Answer 1

После еще одного взгляда я понял, что Роман сделал pts_sp == 1, потому что у меня был только 1 идентификатор для всех моих квадратов, т.е. когда я сделал ID = 1.

Как только я исправил это, я смогстолбец с ID = enum.Для обработки точек в нескольких многоугольниках я могу использовать returnList = TRUE и добавить дополнительные строки для преобразования списка в data.frame, но здесь это не обязательно.

# Attempted sp solution
ptm <- proc.time() # Start timer
# Split the dataframe into a list based on enum and then remove enum from df in the list
polygonlist <- split(localpolydf, localpolydf$enum)
# Convert the list to Polygon, then create a Polygons object
polygonsp <- sapply(polygonlist, function(poly){
  Polygons(list(Polygon(poly[, c("x", "y")])), ID = poly[1, "enum"])
})
# polygonsp <- Polygons(polygonsp, ID = 1)
polygonsp <- SpatialPolygons(polygonsp)
plot(polygonsp)
# Convert points to coordinates
offsetps <- offsetdf
coordinates(offsetps) <- ~x+y
points(offsetps$x, offsetps$y)
# Determine polygons points are in
pts_sp <- over(offsetps, polygonsp)
pts_sp <- data.frame(
  x = offsetps$x, y = offsetps$y,
  enum = unique(localpolydf$enum)[pts_sp])
proc.time() - ptm # end timer

Answer 2

Альтернативой использованию over является использование sf::intersection, поскольку пакет sf становится все более и более популярным.

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

Вот как это сделать:

library(tidyverse)
library(sf)

# convert into st_polygon friendly format (all polygons must be closed)
# must be a nicer way to do this!
localpoly <- localpolydf %>% split(localpolydf$enum) %>% 
lapply(function(x) rbind(x,x[1,])) %>%
lapply(function(x) x[,1:2]) %>%
lapply(function(x) list(as.matrix(x))) %>%
lapply(function(x) st_polygon(x))

# convert points into sf object
points <- st_as_sf(offsetdf,coords=c('x','y'),remove = F)

#convert polygons to sf object and add id column
polys <- localpoly %>% st_sfc() %>% st_sf(geom=.) %>% 
mutate(id=factor(1:100)) 

#find intersection
joined <- polys  %>% st_intersection(points) 


# Sample plot
ggplot() + geom_sf(data=polys) +
geom_sf(data=joined %>% filter(id %in% c(1:10)),aes(col=id)) +
lims(x=c(0,10))

Обратите внимание, что использовать geom_sf ввремя написания вам нужно будет установить версию ggplot для разработки.

вывод графика:

Answer 3

over возвращает индекс точек внутри геометрии. Возможно, что-то вроде этого:

xy <- offsetps[names(na.omit(pts_sp == 1)), ]

plot(polygonsp, axes = 1, xlim = c(0, 10))
points(offsetps)
points(xy, col = "red")