17

我想在栅格地图上识别线性特征，例如道路和河流，并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)

线性特征的 SpatialPolygons 表示

最好的方法是找到穿过多边形的中心线吗？
或者是否有可用的现有代码来执行此操作？

编辑：感谢@mdsumner 指出这称为骨架化。

4

4 回答 4

19

这是我的努力。计划是：

  • 使线条变密
  • 计算 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)))}
于 2012-03-10T01:04:06.763 回答
5

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

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

Voronoi tessellation 在 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 对象： 描述栅格线状特征的 SpatialLines 对象

于 2012-03-09T04:09:56.350 回答
4

您可以通过直接强制将该多边形的边界作为 SpatialLines：

rLines <- as(rPoly, "SpatialLinesDataFrame")

将坐标总结为一条“中心线”是可能的，但我不知道什么是即时的。我认为这个过程通常被称为“骨架化”：

http://en.wikipedia.org/wiki/Topological_skeleton

于 2012-03-07T03:05:51.327 回答
0

我认为理想的解决方案是构建这样的负缓冲区，它动态地达到最小宽度并且在值太大时不会中断；保持持续的对象，并最终在达到该值时画一条线。但不幸的是，这可能对计算要求很高，因为这可能会分步完成并检查特定点的值是否足以拥有一个点（我们的中线）。可能需要无限多的步骤，或者至少需要一些参数化的值。

我现在不知道如何实现这一点。

于 2019-07-19T08:59:15.193 回答