2011-04-20 95 views
8

我想根據高程將多邊形圖層切割成兩部分(上部和下部)。多邊形可能是凸起或凹進的,並且切割的位置可能彼此不同。輪廓線的間隔爲5米,這意味着我可能需要生成輪廓線較濃的輪廓線,例如1米間隔。關於如何做,ArcGIS或R更好的想法? 下面是用於Q運行例如:提前在多邊形圖層下使用等高線切割多邊形

library(sp) 
library(raster) 
r<-raster(ncol=100,nrow=100) 
values(r)<-rep(1:100,100) 
plot(r) ### I have no idea why half of the value is negative... 
p1<-cbind(c(-100,-90,-50,-100),c(60,70,30,30,60)) 
p2<-cbind(c(0,50,100,0),c(0,-25,10,0)) 
p1p<-Polygons(list(Polygon(p1,hole=T)),"p1") 
p2p<-Polygons(list(Polygon(p2,hole=T)),"p2") 
p<-SpatialPolygons(list(p1p,p2p),1:2) 
plot(p,add=T) 
segments(-90,80,-90,20) ##where the polygon could be devided 
segments(50,20,50,-30) ## 

謝謝〜

馬爾科

+2

看看http://cran.r-project.org/web/views/Spatial.html – 2011-04-20 07:17:47

+0

Thanks @ Roman〜我會研究它 – Marco 2011-04-20 07:42:52

+0

你能給我們一些玩具的例子嗎?如果你能告訴我們你正在使用哪些對象/軟件包,這也很好。有多種可能性。 – 2011-04-20 11:27:36

回答

7

如果我理解正確的話,你可以使用rgeos包和相關的空間工具R.

我採取了訣竅來緩衝相交線,然後從該網站生成「差異」多邊形:

http://www.chopshopgeo.com/blog/?p=89

生成示例柵格和覆蓋的多邊形。

vdata <- list(x = 1:nrow(volcano), y = 1:ncol(volcano), z = volcano) 

## raw polygon data created using image(vdata); xy <- locator() 

xy <- structure(list(x = c(43.4965355534823, 41.7658494766076, 36.2591210501883, 
25.560334393145, 13.7602020508178, 18.7949251835441, 29.179041644792, 
40.6645037913237, 44.2832110429707, 47.272577903027, 47.5872480988224 
), y = c(30.0641086410103, 34.1278207016757, 37.6989616034726, 
40.900674136118, 32.7732500147872, 27.4781100569505, 22.5523984682652, 
22.7986840476995, 24.5226831037393, 29.3252519027075, 33.8815351222414 
)), .Names = c("x", "y")) 

## close the polygon 
coords <- cbind(xy$x, xy$y) 
coords <- rbind(coords, coords[1,]) 

library(sp) 

## create a Spatial polygons object 
poly <- SpatialPolygons(list(Polygons(list(Polygon(coords, hole = FALSE)), "1"))) 


## create a contour line that cuts the polygon at height 171 
cl <- contourLines(vdata, levels = 171) 

## for ContourLines2SLDF 
library(maptools) 

clines <- ContourLines2SLDF(cl) 

現在,相交線的多邊形,然後稍微緩衝行和差再次與多邊形給予多聚。

library(rgeos) 
lpi <- gIntersection(poly, clines) 

blpi <- gBuffer(lpi, width = 0.000001) 

dpi <- gDifference(poly, blpi) 

繪製原始數據,並從空間對象手動提取多邊形的一半。

par(mfrow = c(2,1)) 

image(vdata) 
plot(poly, add = TRUE) 

plot(SpatialPolygons(list(Polygons(list([email protected][[1]]@Polygons[[1]]), "1"))), 
    add = TRUE, col = "lightblue") 

image(vdata) 
plot(poly, add = TRUE) 
cl <- contourLines(vdata, levels = 171) 

plot(SpatialPolygons(list(Polygons(list([email protected][[1]]@Polygons[[2]]), "2"))), 
    add = TRUE, col = "lightgreen") 

這適用於這個相當簡單的情況,它可能對您的情況有用。

+0

Thanks @ mdsumner ~~我會在我的案例上測試它:) – Marco 2011-05-05 10:55:49