Wie überprüft man, ob ein Punkt in einem Polygon liegt und R für große Datenmengen verwendet?

Question

Ich bin neu in R und für mein momentanes Projekt muss ich eine Heatmap zeichnen, die sich auf ein bestimmtes Ereignis bezieht. Es gibt ungefähr 2 Millionen Beobachtungen eines solchen Ereignisses und bei jeder Beobachtung gibt es eine Längs- und eine Breitenkoordinate. Auch habe ich die Kartendaten in einen Datenrahmen umgewandelt und der Datenrahmen enthält 71 Bezirk, jeder Bezirk wird mit einer Menge von Koordinaten definiert. Ich muss entscheiden, welche Beobachtung der Veranstaltung zu welchem ​​Bezirk gehört. Ich bin mit dem folgenden Code:

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 die für das Ereignis gesetzt Daten vorhanden sind, ist polygonOfdis die für jeden Stadtteil Daten.

Für kleine Datenmenge funktioniert dieser Algorithmus in Ordnung, aber für meinen Datensatz wird es definitiv mehr als zehn Stunden oder mehr laufen (Für einen Datensatz nur 1/400 der aktuellen Größe läuft dieser Algorithmus für 1 bis 2 Protokoll). Ich frage mich, ob es irgendeinen besseren Weg gibt herauszufinden, welche Beobachtung zu welchem ​​Bezirk gehört? Mein Problem ist, dass die point.in.polygon Funktion zu viel Zeit braucht und ich frage mich, ob es irgendeine andere Funktion kann das tun?

PS: Die aktuellen Daten sind tatsächlich nur 1/10 der realen Daten, die ich verarbeiten muss, also brauche ich wirklich einen schnelleren Weg, dies zu tun.

Answer 1

Also, vor einer Weile portierte ich über einen Punkt in einem Polygon-Algorithmus von W. Randolph Franklin, die den Begriff der Strahlen verwendet. I.e. Wenn sich ein Punkt im Polygon befindet, sollte er ungerade Male durchlaufen werden. Andernfalls, wenn es eine gerade Zahl hat, sollte es auf der Außenseite des Polygons liegen.

Der Code ist sehr schnell, da er mit Rcpp geschrieben wird. Es ist in zwei Teile aufgeteilt: 1. Der PIP-Algorithmus und 2. Eine Wrapper-Funktion für die Klassifizierung.

PIP-Algorithmus

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

Klassifizierungsalgorithmus

//' 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

Ihr Code ist ziemlich geradlinig, Ihr Stein des Anstoßes ist die Verwendung von Schleifen statt der Vektorisierung Leistung des R. Dieser Code sollte funktionieren, aber ohne Daten Ich kann es nicht überprüfen kann:

# 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) 

Dieser Code geht auch davon aus, dass jeder Punkt in einer und nur 1 Polygon ist. Die innere Schleife und das Tapply könnten höchstwahrscheinlich durch Verwendung der dplyr-Bibliothek verbessert werden. Die andere aufgelistete Lösung mit dem PIP-Algorithmus könnte einen Schub gegenüber der integrierten Methode bieten.

Answer 3

Es gibt ein Paket dafür, nämlich 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 

Beachten Sie, dass mehrere Punkte (siehe unten) testen können, aber wenn Sie eine einzelne Sie eine Matrix müssen testen, das ist, warum ich rbind verwenden.

Sie erhalten 0, wenn der Punkt innerhalb des Polygons ist, -1 anders:

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

Wie gesagt, bevor Sie gleichzeitig mehrere Punkte testen:

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

Das Paket auch für testen können Punkteindämmung in einem 3D-Polyeder.