r - 使用 R 的多个多边形中的点

Question

目前我有两个data.frame，一个是多边形(poly.x, poly.y, enum)，一个(pt.x, pt.y)是多边形enum的id。我试图确定哪些点属于哪些多边形，所以我得到了一个 data.frame (pt.x, pt.y, enum)。

我的第一次尝试使用point.in.polygon包sp和lapply函数来查找该点所属的多边形。虽然我的代码有效，但在大型数据集上需要很长时间。

我的第二次尝试over也从sp包中使用，由关于 gis stackexchange 的问题拼凑而成。虽然它要快得多，但我似乎无法从中获得正确的输出，over因为它是1s 和NAs 的数据帧。

下面我包含了一个简化的工作示例（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 

标签： rgispoint-in-polygon

Answer 1

over返回几何内点的索引。也许是这样的：

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

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