Определите линейную функцию на растровой карте и верните объект линейной формы, используя R

Ответ 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)

}

Вместо буферизации более ранней версии я вычисляю расстояние от каждой середины до линии многоугольника и беру только точки, которые являются а) внутри, и б) дальше от края, чем 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){[email protected][[1]]@Polygons[[1]]@coords}

Захват пространственных линий с функцией "захвата":

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

Ответ 2

Вы можете получить границу этого многоугольника как SpatialLines прямым принуждением:

rLines <- as(rPoly, "SpatialLinesDataFrame")

Подведение итогов координат до одной "центральной линии" было бы возможно, но ничего непосредственного, о котором я не знаю. Я думаю, что этот процесс обычно называют "скелетонизацией":

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

Ответ 3

Благодаря @klewis на gis.stackexchange.com для ссылки на этот элегантный алгоритм для нахождения центральной линии (в ответ на соответствующий вопрос, который я задал там).

Процесс требует нахождения координат на краю многоугольника, описывающего линейную особенность, и выполнения тероселяции Вороного этих точек. Координаты плиток Вороного, которые попадают в многоугольник линейной функции, попадают на осевую линию. Поверните эти точки в строку.

Теоретизация Вороного выполняется эффективно в 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:

Ответ 4

Я думаю, что идеальным решением было бы создать такой отрицательный буфер, который динамически достигает минимальной ширины и не ломается, когда значение слишком велико; сохраняет продолжение объекта и в конечном итоге рисует линию, если значение достигнуто. Но, к сожалению, это может быть очень требовательным к вычислениям, потому что это будет сделано, вероятно, поэтапно и проверяет, достаточно ли значения для конкретной точки, чтобы иметь точку (нашей средней линии). Возможно, ему не нужно иметь бесконечное число шагов или, по крайней мере, какое-то параметризованное значение.

Я не знаю, как реализовать это сейчас.