识别栅格地图上的线性特征并使用R返回线性对象 [英] Identify a linear feature on a raster map and return a linear shape object using R

问题描述

我想在栅格地图上识别诸如道路和河流之类的线性特征，并使用R将其转换为线性空间对象(SpatialLines类).

I would like to identify linear features, such as roads and rivers, on raster maps and convert them to a linear spatial object (SpatialLines class) using R.

raster和sp包可用于将要素从栅格转换为多边形矢量对象(SpatialPolygons类). rasterToPolygons()将从栅格中提取特定值的像元并返回一个多边形对象.可以使用dissolve=TRUE选项简化该产品，该选项调用rgeos程序包中的例程来完成此操作.

The raster and sp packages can be used to convert features from rasters to polygon vector objects (SpatialPolygons class). rasterToPolygons() will extract cells of a certain value from a raster and return a polygon object. The product can be simplified using the dissolve=TRUE option, which calls routines in the rgeos package to do this.

这一切都很好，但是我希望它是一个SpatialLines对象.我该怎么办?

This all works just fine, but I would prefer it to be a SpatialLines object. How can I do this?

请考虑以下示例:

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

最好的方法是找到穿过多边形的中心线吗?
还是有可用的现有代码来做到这一点?

Would the best approach be to find a centre line through the polygon?
Or is there existing code available to do this?

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

Thanks to @mdsumner for pointing out that this is called skeletonization.

推荐答案

这是我的努力.该计划是:

Here's my effort. The plan is:

• 压实线条
• 计算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的点最远离边缘的点.

Instead of the buffering of the earlier version I compute the distance from each midpoint to the line of the polygon, and only take points that are a) inside, and b) further from the edge than 1.5 of the distance of the inside point that is furthest from the edge.

如果多边形以较长的节距并没有致密化地折回自身，则可能会出现问题.在这种情况下，图形是一棵树，代码会报告它.

Problems can arise if the polygon kinks back on itself, with long segments, and no densification. In this case the graph is a tree and the code reports it.

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

As a test, I digitized a line (s, SpatialLines object), buffered it (p), then computed the centreline and superimposed them:

 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事物中获取一个环的坐标，而该空间只能是一个环:

The 'onering' function just gets the coordinates of one ring from a SpatialPolygons thing that should only be one ring:

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

使用捕获"功能捕获空间线特征:

Capture spatial lines features with the 'capture' function:

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

这篇关于识别栅格地图上的线性特征并使用R返回线性对象的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！