2015-02-08 4 views
0

Я хотел бы выяснить, какие точки, которые определяют правильную решетку, находятся в пределах многоугольника. Приведенный ниже код делает это, но ОЧЕНЬ ОЧЕНЬ медленно:Определение того, какие точки в регулярной решетке находятся внутри границ многоугольника

#the polygon that I want to check each point against 
glasgow_single <- readShapePoly(
    fn="data/clipped/glasgow_single" 
) 

#interpolated contains the coordinates of the regular grid 
points_to_check <- expand.grid(
    x=interpolated$x, 
    y=interpolated$y 
    ) 

#function to be called by plyr  
fn <- function(X){ 
    this_coord <- data.frame(lon=X["x"], lat=X["y"]) 
    this_point <- SpatialPoints(this_coord) 
    out <- gContains(glasgow_single, this_point) 
    out <- data.frame(x=X["x"], y=X["y"], val=out) 
    return(out) 
} 

#plyr call 
vals <- adply(points_to_check, 1, fn, .progress="text") 
vals$val <- as.numeric(vals$val) 

Принимая во внимание как время на обдумывание и времени вычисления, есть гораздо более быстрый способ сделать это?

ответ

1

Да, есть намного лучший подход. Для этого и многих других топологических операций пакет rgeos вы хорошо накрыли. Здесь Вы желаете rgeos::gWithin():

## Required packages 
library(rgdal) 
library(raster) ## For example polygon & functions used to make example points 
library(rgeos) 

## Reproducible example 
poly <- readOGR(system.file("external", package="raster"), "lux")[1,] 
points <- as(raster(extent(poly)), "SpatialPoints") 
proj4string(points) <- proj4string(poly) 

## Test which points fall within polygon 
win <- gWithin(points, poly, byid=TRUE) 

## Check that it works 
plot(poly) 
points(points, col=1+win) 

enter image description here

+0

Отлично. Большое спасибо. – JonMinton