2012-02-22 16 views
7

Mam wiele zestawów punktów (dla różnych lat ~ 20)Jak tworzyć wieloboki thiessen z punktów za pomocą pakietów R?

Chcę generować wieloboki Thiessen dla każdego zestawu punktów za pomocą pakietów przestrzennych r.

wiem, można to zrobić za pomocą GIS, ale jak chcesz coś proces wsadowy w R byłoby

pomocne.

+2

'install.packages ("SOS"); biblioteka ("sos"); findFn ("thiessen") ' –

+0

Używam ArcGIS do tego samego teraz ... – user2760

Odpowiedz

16

Nie dałeś nam dostępu do twoich danych, ale oto przykład dla punktów reprezentujących miasta świata, używając podejścia opisanego przez Carson Farmer pod numerem his blog. Mam nadzieję, że będzie to dobry początek ...

# Carson's Voronoi polygons function 
voronoipolygons <- function(x) { 
    require(deldir) 
    require(sp) 
    if (.hasSlot(x, 'coords')) { 
    crds <- [email protected] 
    } else crds <- x 
    z <- deldir(crds[,1], crds[,2]) 
    w <- tile.list(z) 
    polys <- vector(mode='list', length=length(w)) 
    for (i in seq(along=polys)) { 
    pcrds <- cbind(w[[i]]$x, w[[i]]$y) 
    pcrds <- rbind(pcrds, pcrds[1,]) 
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i)) 
    } 
    SP <- SpatialPolygons(polys) 
    voronoi <- SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1], 
    y=crds[,2], row.names=sapply(slot(SP, 'polygons'), 
    function(x) slot(x, 'ID')))) 
} 

Przykład 1: Wejście jest SpatialPointsDataFrame:

# Read in a point shapefile to be converted to a Voronoi diagram 
library(rgdal) 
dsn <- system.file("vectors", package = "rgdal")[1] 
cities <- readOGR(dsn=dsn, layer="cities") 

v <- voronoipolygons(cities) 

plot(v) 

Voronoi diagram of cities

Przykład 2: Wejście to wektory x, y współrzędne:

dat <- data.frame(x=runif(100), y=runif(100)) 
v2 <- voronoipolygons(dat) 
plot(v2) 

Another voronoi diagram

+0

Zmodyfikowałem funkcję tak, aby akceptowała wektory współrzędnych (oczekiwane w porządku x, y) oraz' SpatialPointsDataFrame'. – jbaums

+0

To działało..TY! – user2760

+0

Miło to słyszeć. – jbaums

1

samej zasadzie, jak pokazano przez jbaums, ale prostszy kod:

library(dismo) 
library(rgdal) 
cities <- shapefile(file.path(system.file("vectors", package = "rgdal")[1], "cities")) 

v <- voronoi(cities) 
plot(v) 
Powiązane problemy