识别栅格地图上的线性要素并使用R返回线性形状对象

Question

我想在栅格地图上识别线性要素,例如道路和河流,并SpatialLines使用R 将它们转换为线性空间对象(类).

的raster和sp封装可以用于从特征栅格到多边形矢量对象(转换SpatialPolygons类别). rasterToPolygons()将从栅格中提取特定值的单元格并返回多边形对象.可以使用dissolve=TRUE选项简化产品,该选项调用rgeos包中的例程来执行此操作.

这一切都很好,但我更喜欢它是一个SpatialLines对象.我怎样才能做到这一点？

考虑这个例子:

## Produce a sinuous linear feature on a raster as an example
library(raster)
r <- raster(nrow=400, ncol=400, xmn=0, ymn=0, xmx=400, ymx=400)
r[] <- NA
x <-seq(1, 100, by=0.01)
r[cellFromRowCol(r, round((sin(0.2*x) + cos(0.06*x)+2)*100), round(x*4))] <- 1

## Quick trick to make it three cells wide
r[edge(r, type="outer")] <- 1

## Plot
plot(r, legend=FALSE, axes=FALSE)

## Convert linear feature to a SpatialPolygons object
library(rgeos)
rPoly <- rasterToPolygons(r, fun=function(x) x==1, dissolve=TRUE)
plot(rPoly)

最好的方法是找到穿过多边形的中心线吗？
或者是否有可用的代码？

编辑:感谢@mdsumner指出这称为骨架化.

Answer 1

这是我的努力.计划是:

• 密集线
• 计算一个delaunay三角剖分
• 取中点,并取出多边形中的那些点
• 构建距离加权最小生成树
• 找到它的图形直径路径

初学者的密集代码:

densify <- function(xy,n=5){
  ## densify a 2-col matrix
  cbind(dens(xy[,1],n=n),dens(xy[,2],n=n))
}

dens <- function(x,n=5){
  ## densify a vector
  out = rep(NA,1+(length(x)-1)*(n+1))
  ss = seq(1,length(out),by=(n+1))
  out[ss]=x
  for(s in 1:(length(x)-1)){
    out[(1+ss[s]):(ss[s+1]-1)]=seq(x[s],x[s+1],len=(n+2))[-c(1,n+2)]
  }
  out
}

现在主要课程:

simplecentre <- function(xyP,dense){
require(deldir)
require(splancs)
require(igraph)
require(rgeos)

### optionally add extra points
if(!missing(dense)){
  xy = densify(xyP,dense)
} else {
  xy = xyP
}

### compute triangulation
d=deldir(xy[,1],xy[,2])

### find midpoints of triangle sides
mids=cbind((d$delsgs[,'x1']+d$delsgs[,'x2'])/2,
  (d$delsgs[,'y1']+d$delsgs[,'y2'])/2)

### get points that are inside the polygon 
sr = SpatialPolygons(list(Polygons(list(Polygon(xyP)),ID=1)))
ins = over(SpatialPoints(mids),sr)

### select the points
pts = mids[!is.na(ins),]

dPoly = gDistance(as(sr,"SpatialLines"),SpatialPoints(pts),byid=TRUE)
pts = pts[dPoly > max(dPoly/1.5),]

### now build a minimum spanning tree weighted on the distance
G = graph.adjacency(as.matrix(dist(pts)),weighted=TRUE,mode="upper")
T = minimum.spanning.tree(G,weighted=TRUE)

### get a diameter
path = get.diameter(T)

if(length(path)!=vcount(T)){
  stop("Path not linear - try increasing dens parameter")
}

### path should be the sequence of points in order
list(pts=pts[path+1,],tree=T)

}

而不是缓冲早期版本,我计算从每个中点到多边形线的距离,并且只取a)内部的点,和b)距边缘的距离比内部点的距离的1.5更远.离边缘最远.

如果多边形自身扭结,长段并且没有致密化,则可能出现问题.在这种情况下,图形是树,代码报告它.

作为测试,我将一行(s,SpatialLines对象)数字化,缓冲它(p),然后计算中心线并叠加它们:

 s = capture()
 p = gBuffer(s,width=0.2)
 plot(p,col="#cdeaff")
 plot(s,add=TRUE,lwd=3,col="red")
 scp = simplecentre(onering(p))
 lines(scp$pts,col="white")

'onering'函数只从一个SpatialPolygons事物中得到一个环的坐标,它应该只是一个环:

onering=function(p){p@polygons[[1]]@Polygons[[1]]@coords}

使用"捕获"功能捕获空间线要素:

capture = function(){p=locator(type="l")
            SpatialLines(list(Lines(list(Line(cbind(p$x,p$y))),ID=1)))}

Answer 2

感谢gis.stackexchange.com上的 @klewis链接到这个用于查找中心线的优雅算法（回应我在那里提出的相关问题）。

该过程需要找到描述线性特征的多边形边缘上的坐标，并对这些点执行 Voronoi 细分。落在线性要素多边形内的 Voronoi 切片的坐标落在中心线上。将这些点变成一条线。

Voronoi 曲面细分在 R 中使用该包非常高效地完成deldir，并且使用该rgeos包进行多边形和点的交集。

## Find points on boundary of rPoly (see question)
rPolyPts <-  coordinates(as(as(rPoly, "SpatialLinesDataFrame"),
                "SpatialPointsDataFrame"))

## Perform Voronoi tessellation of those points and extract coordinates of tiles
library(deldir)
rVoronoi <- tile.list(deldir(rPolyPts[, 1], rPolyPts[,2]))
rVoronoiPts <- SpatialPoints(do.call(rbind, 
                 lapply(rVoronoi, function(x) cbind(x$x, x$y))))

## Find the points on the Voronoi tiles that fall inside 
## the linear feature polygon
## N.B. That the width parameter may need to be adjusted if coordinate
## system is fractional (i.e. if longlat), but must be negative, and less
## than the dimension of a cell on the original raster.
library(rgeos)
rLinePts <- gIntersection(gBuffer(rPoly, width=-1), rVoronoiPts)

## Create SpatialLines object
rLine <- SpatialLines(list(Lines(Line(rLinePts), ID="1")))

生成的 SpatialLines 对象：