Разрезать полигоны, используя контурную линию под слоями полигонов - PullRequest
Новая
       37

Разрезать полигоны, используя контурную линию под слоями полигонов

9 голосов
/ 20 апреля 2011

Я хотел бы разрезать слой многоугольника, в зависимости от высоты, на две части (верхняя и нижняя часть). Многоугольник может выпуклым или вогнутым, и положение для резки может отличаться друг от друга. Контурная линия имеет интервал 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)  ##

Заранее спасибо ~

Marco

1 Ответ

8 голосов
/ 03 мая 2011

Если я правильно понимаю, вы можете использовать пакет rgeos и связанные с ним инструменты Spatial в 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(dpi@polygons[[1]]@Polygons[[1]]), "1"))), 
     add = TRUE, col = "lightblue")

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

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

Это работает для этого довольно простого случая, это может быть полезно для вашего сценария.

...