r - 如何使用R对大型数据集有效地检查点是否在多边形中？

Question

我是 R 的新手，对于我当前的项目，我必须绘制与特定事件相关的热图。大约有 200 万次此类事件的观测，并且在每个观测中都有一个长坐标和纬度坐标。另外，我将 map 数据转换为数据框，数据框包含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 移植了多边形算法中的一个点使用射线的概念。 IE。如果一个点在多边形中，它应该通过奇数次。否则，当它有偶数时，它应该位于多边形的外部。

代码相当快，因为​​它是使用 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;
}