如何使用R有效地检查大数据集中的点是否在多边形内？

Question

我是R语言的新手，在我的项目中，我需要绘制与特定事件相关的热图。这个事件有大约200万个观测值，每个观测值都有一个经纬度坐标。此外，我已经将地图数据转换为数据框，该数据框包含71个区域，每个区域都用一组坐标定义。我需要确定事件的哪个观测值属于哪个区域。我正在使用以下代码：

for (row in 1:nrow(data2015)){
  point.x=data2015[row,"Latitude"]
  point.y=data2015[row,"Longitude"]
  for (name in names(polygonOfdis)){
    if (point.in.polygon(point.x, point.y, polygonOfdis[[name]]$lat,   polygonOfdis[[name]]$long, mode.checked=FALSE)){
    count[[name]]<-count[[name]]+1
    break
}
}
}

data2015是事件的数据集，polygonOfdis是每个区域的数据集。

对于小型数据集，此算法效果良好，但对于我的数据集，它肯定会运行超过十个小时甚至更长时间（对于当前大小的数据集的1/400，此算法运行时间为1到2分钟）。我想知道有没有更好的方法来查找哪个观察值属于哪个区域？我的问题是point.in.polygon函数需要太多时间，我想知道是否有其他函数可以完成此操作？

注：当前数据实际上只有我需要处理的真实数据的1/10，因此我真的需要一种更快的方法来处理这个问题。

- Optimus Prime

Answer 1

不久前，我移植了一个使用射线概念的多边形中的点算法，该算法由W. Randolph Franklin编写。即，如果一个点在多边形内部，它应该通过奇数次。否则，当它有偶数次时，它应该位于多边形外部。

代码非常快，因为它使用Rcpp编写。它分为两个部分：1. PIP算法和2.用于分类的包装函数。

PIP算法

#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]

//' @param points A \code{rowvec} with x,y  coordinate structure.
//' @param bp     A \code{matrix} containing the boundary points of the polygon. 
//' @return A \code{bool} indicating whether the point is in the polygon (TRUE) or not (FALSE)
// [[Rcpp::export]]
bool pnpoly(const arma::rowvec& point, const arma::mat& bp) {
    // Implementation of the ray-casting algorithm is based on
    // 
    unsigned int i, j;
    
    double x = point(0), y = point(1);
    
    bool inside = false;
    for (i = 0, j = bp.n_rows - 1; i < bp.n_rows; j = i++) {
      double xi = bp(i,0), yi = bp(i,1);
      double xj = bp(j,0), yj = bp(j,1);
      
      // See if point is inside polygon
      inside ^= (((yi >= y) != (yj >= y)) && (x <= (xj - xi) * (y - yi) / (yj - yi) + xi));
    }
    
    // Is the cat alive or dead?
    return inside;
}

分类算法

//' PIP Classifier
//' @param points A \code{matrix} with x,y  coordinate structure.
//' @param names  A \code{vector} of type \code{string} that contains the location name.
//' @param bps    A \code{field} of type {matrix} that contains the polygon coordinates to test against.
//' @return A \code{vector} of type \code{string} with location information.
// [[Rcpp::export]]
std::vector<std::string> classify_points(const arma::mat& points, 
                                         std::vector<std::string> names,
                                         const arma::field<arma::mat>& bps){
  unsigned int i, j;
  
  unsigned int num_points = points.n_rows;
  
  std::vector<std::string> classified(num_points);
  
  for(i = 0; i < num_points; i++){
    
    arma::rowvec active_row = points.row(i);
    
    // One of the coordinate lacks a value
    if( !arma::is_finite(active_row(0)) || !arma::is_finite(active_row(1)) ){
      classified[i] = "Missing";
      continue; // skip trying to find a location
    }

    // Try to classify coordinate based on supplied boundary points for area j
    for(j = 0; j < names.size(); j++){
      if( pnpoly(active_row, bps(j)) ){
        classified[i] = names[j];
        break; // Break loop
      }
    }
    
  }
  
  return classified;
}

- coatless

Answer 2

这个函数来自SMDTools软件包，表现良好。

- Conner M.

---

链接现在已经失效* - Brennan

“SDMTools”软件包已从CRAN存储库中移除。 - Jeff Bezos

Answer 3

我刚发现这个对我很有效：

library(secr)
## 100 random points in unit square
xy <- matrix(runif(200, -0.5, 1.5), ncol = 2)
## triangle centred on (0.5, 0.5)
poly <- data.frame(x = c(1, 1, 0, 0, 1), y = c(1,0,0, 1, 1))
plot(xy, pch = 1 + pointsInPolygon(xy, poly))
lines(poly)

- JAQuent

Answer 4

有一个专门的包可以解决这个问题，它叫做ptinpoly。

library(ptinpoly)
# define a square 
square <- rbind(
  c(0,0),
  c(0,1),
  c(1,0),
  c(1,1)
)

pinside <- rbind(c(0.5,0.5)) # point inside the square
poutside <- rbind(c(2,1)) # point outside the square

请注意，您可以测试多个点（见下文），但如果您只测试一个点，则需要矩阵，这就是为什么我使用rbind的原因。

如果该点在多边形内部，则返回0，否则返回-1：

> pip2d(square, pinside)
[1] 0
> pip2d(square, poutside)
[1] -1

就像我之前所说的，你可以同时测试多个点：

> pip2d(square, rbind(pinside, poutside))
[1]  0 -1

该软件包还允许在3D多面体中测试点的包含性。

- Stéphane Laurent

---

我发现pip2d函数不太可靠，会导致RStudio会话崩溃并出现致命错误。 - Conner M.

@ConnerM。我最近进行了深入的测试，并确认了您的评论。 - Stéphane Laurent

我尝试了很多次，但在以下示例中它无法正常工作： square <- rbind( c(-26.606963 ,16.02765), c(-26.606963, 46.02765, c(-6.606963, 46.02765), c(-6.606963 ,16.02765) ) pip2d(square,rbind(c(-24.18359 36.11769))) #应该是内部，给出-1 pip2d(square,rbind(c(50.38219 11.98086))) #应该是外部，也给出-1 - Rosanne

Pip2d仍然会导致R会话崩溃 - 是否已经找到了任何修复或解决方案？ - jpd527

pip2d 是无用的，会导致 R 会话崩溃 :/ - Shrimp

是的，这段代码似乎不起作用。一般来说，似乎正方形对此无效？如果我模拟你的正方形并绘制它，我得到三个三角形而不是一个正方形。 - JAQuent

Answer 5

你的代码非常直观，但你使用了循环而不是 R 的向量化功能，这是你遇到的障碍。这段代码应该可以工作，但没有任何数据我无法验证：

# create a column onto the dataframe to store the results 
data2015$poly<-"blank"
point.x=data2015$Latitude
point.y=data2015$Longitude
for (name in names(polygonOfdis)){
    #point.in.polygon returns a arrary of 0 to 3 for point location
    inpoly<-point.in.polygon(point.x, point.y, polygonOfdis[[name]]$lat,
                         polygonOfdis[[name]]$long, mode.checked=FALSE)
    #if the element in >0 in poly assign poly name to poly column 
    data2015$poly[inpoly>0]<-name
  }
  #additional processing (returns count per polygon)
  tapply(data2015$poly, INDEX = data2015$poly, FUN=length)

这段代码还假设每个点只属于一个多边形。使用 dplyr 库可能可以改进内部循环和 tapply 的效率。另一个使用 PIP 算法的解决方案可能比内置方法更有效。

- Dave2e

Answer 6

您可以使用cgalPolygons包（尚未在CRAN上发布）。

library(cgalPolygons)
# define a square 
square <- rbind(
  c(0, 0),
  c(0, 1),
  c(1, 0),
  c(1, 1)
)

pinside   <- c(0.5, 0.5) # point inside the square
poutside  <- c(2, 1)     # point outside the square
ponsquare <- c(1, 0.5)   # point on the boundary of the square

请注意，您可以使用单个命令测试多个点（参见下文）。

如果点位于多边形内部，则返回1，如果在外部则返回-1，如果在边界上则返回0：

> plg <- cgalPolygon$new(square)
> plg$whereIs(pinside)
[1] 1
> plg$whereIs(poutside)
[1] -1
> plg$whereIs(ponsquare)
[1] 0

就像我之前说的那样，您可以同时测试多个点：

> plg$whereIs(rbind(pinside, poutside, ponsquare))
[1]  1 -1  0

- Stéphane Laurent

Answer 7

根据@conner-m的建议：

library(tidyverse)
library(furrr)
library(SMDTools)

plan(multiprocess)
future_map2_dfr(
  polygonOfdis,
  names(polygonOfdis),
  ~tibble(
    district = .y,
    pip = 
      pnt.in.poly(
        data2015[, c('Latitude', 'Longitude')], 
        .x
      )$pip
  )
) %>% 
  group_by(district) %>% 
  summarise(count = sum(pip))

- Augusto Fadel

Answer 8

我更多地涉及空间数据。我会将它们转换为空间对象以执行操作（对我来说需要更少的时间，因为可能不是非常高效）

xyDf <- data.frame(X = MyYPtsCoordshere), Y = MyYPtsCoordshere) # points coords
coordMat <- data.frame(X = MyYPolygonCoordshere, Y = MyYPolygonCoordshere) # polygon coords

## Filter points by bounding box (easy=
posCoord <- which(
     xyDf$X <= max(coordMat$X) & # west
     xyDf$X >= min(coordMat$X) & # east
     xyDf$Y <= max(coordMat$Y) & # north
     xyDf$Y >= min(coordMat$Y) )# south

#check how many: good for debug      
str(posCoord)
plot(coordMat[, c('X', 'Y')], type = 'b')
points(xyDf[, c('X', 'Y')], col = 2, pch = 20)
points(xyDf[posCoord, c('X', 'Y')], col = 4, pch = 2)
      
# Filter for real using the bbox
xySel <- xyDf[posCoord, ]

#Make the polygon spatial
spDf <<- sp:SpatialPolygonsDataFrame(
        SpatialPolygons(list(Polygons(list(Polygon(coordMat)), 1) # polgons
          )), data = data.frame(ID = 1), match.ID = FALSE)

#Make póints spatial and make the query
posSel <- sp::over(sp::SpatialPoints(xySel[, c('X', 'Y')]), spDf)

posSel 是一个数据框，指示每个点是否在内部，值为1或NA。另一个选项（速度较慢）是使用raster包+extract函数。

system.time(ov_ap_mat <- sp::over(coord_pts, ap)) # 0.45
#system.time(ov_ap_matR <- raster::extract( ap_eco_sp, mat[, c('X_a1', 'Y_a1')])) # 1.25

- gonzalez.ivan90