r - 识别栅格 map 上的线性特征并使用 R 返回线性形状对象

Question

我想在栅格 map 上识别线性特征，例如道路和河流，并使用 R 将它们转换为线性空间对象(SpatialLines 类)。
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

这是我的努力。计划是:

增密线条

计算德劳内三角剖分

取中点，并取多边形中的那些点

构建距离加权最小生成树

找到它的图直径路径

初学者的致密代码:

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)))}