如何使用 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 函数花费太多时间，我想知道是否有其他函数可以做到这一点？

PS：当前的数据实际上只是我必须处理的真实数据的1/10，所以我真的需要一种更快的方法来做到这一点。

Answer 1

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

该代码相当快，因为​​它是使用Rcpp编写的。它分为两部分：1. PIP 算法和 2. 用于分类的包装函数。

画中画算法

#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;
}