大規模なデータセットに対してポイントをポリゴンに入れて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は、各地区のデータ・セットで、イベント用のデータセットです。

小規模なデータセットの場合、このアルゴリズムは問題なく動作しますが、データセットでは10時間以上実行されます（現在のサイズの1/4のデータセットの場合、このアルゴリズムは1〜2分）。どんな所見がどの地区に属するのかを知る良い方法があれば、私は疑問に思っていますか？私の問題は、point.in.polygon関数が時間がかかり過ぎて、これを行う他の関数があるかどうか疑問に思っていますか？

PS：現在のデータは実際に処理する必要がある実際のデータの1/10に過ぎないため、これを行うにはもっと速い方法が本当に必要です。

Answer 1

それで、しばらく前に、光線の概念を使用するW. Randolph Franklinによってポリゴンアルゴリズムのポイントを移植しました。私。点がポリゴン内にある場合は、奇数回通過する必要があります。さもなければ、それが偶数であるとき、それはポリゴンの外側にあるべきです。

Rcppを使用して記述されているため、コードはかなり高速です。これは2つの部分に分割されています。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; 
} 

Answer 2

あなたのコードは、かなり単純です、あなたのつまずきの石ではなく、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) 

このコードでは、各ポイントが1つのポリゴンしか含まないと仮定しています。内部ループとtapplyは、おそらくdplyrライブラリを使用することで改善されました。 PIPアルゴリズムを使用したもう1つのソリューションでは、組み込みのメソッドを上回る効果が得られます。

Answer 3

それに対応するパッケージ、つまりptinpolyがあります。私はrbindを使用する理由あなたが行列を必要とする単一のものをテストする場合は、（下記参照）のいくつかのポイントをテストすることができますが、

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 

注意、それはです。ポイントは多角形の内部にある場合

あなたが0取得、-1そうでない場合：

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

あなたは、同時に複数の点をテストすることができます前に、私が言ったように：

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

パッケージはまたのためにテストすることができます

3次元多面体の点封じ込め