Identificare una caratteristica lineare su una mappa raster e restituire un oggetto forma lineare utilizzando R

Question

Vorrei identificare le caratteristiche lineari, come strade e fiumi, su mappe raster e convertirle in un oggetto territoriale lineare (classe) utilizzando R

I pacchetti raster e sp possono essere utilizzati per convertire le funzioni dai raster agli oggetti vettoriali poligono (classe SpatialPolygons). rasterToPolygons() estrarrà celle di un certo valore da un raster e restituirà un oggetto poligono. Il prodotto può essere semplificato utilizzando l'opzione dissolve=TRUE, che chiama le routine nel pacchetto rgeos per farlo.

Questo funziona tutto bene, ma preferirei che fosse un oggetto . Come posso fare questo?

Considerate questo esempio:

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

Sarebbe il miglior approccio essere quello di trovare una linea mediana del poligono?
Oppure esiste un codice esistente per farlo?

EDIT: Grazie a @mdsumner per aver sottolineato che questo è chiamato scheletrizzazione.

Answer 1

Ecco il mio sforzo.Il piano è:

• densificare le linee
• calcolare una triangolazione delaunay
• prendere i punti medi, e prendere quei punti che sono poligono
• costruire un minimo albero ricoprente distanza ponderata
• trovare il suo percorso di diametro grafico

Il codice ridensificante per cominciare:

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 
} 

E ora il corso principale:

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) 

} 

Invece del buffer della versione precedente ho calcolare la distanza da ogni punto medio alla linea del poligono, e solo prendere i punti che sono a) all'interno e b) più lontano dal bordo di 1,5 della distanza del punto interno più lontano dal bordo.

I problemi possono sorgere se il poligono si piega su se stesso, con segmenti lunghi e senza densificazione. In questo caso il grafico è un albero e il codice lo segnala.

Come prova, ho digitalizzato una linea (s, SpatialLines oggetto), tamponata (p), poi calcolato mezzeria e sovrapposti:

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

La funzione 'onering' ottiene solo le coordinate di un anello da una cosa SpatialPolygons che dovrebbe essere un solo anello:

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

cattura linee spaziali caratteristiche con la funzione di 'cattura' :

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

Answer 2

È possibile ottenere il confine di quel poligono come SpatialLines di coercizione diretta:

rLines <- as(rPoly, "SpatialLinesDataFrame") 

Riassumendo le coordinate fino a un singolo "linea centrale" sarebbe possibile, ma nulla immediato, che io sappia. Credo che questo processo è generalmente chiamato "scheletrizzazione":

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

Answer 3

Grazie a @klewis a gis.stackexchange.com per il collegamento a this elegant algorithm per trovare la linea centrale (in risposta ad una domanda correlata ho chiesto lì).

Il processo richiede la ricerca delle coordinate sul bordo di un poligono che descriva la caratteristica lineare e l'esecuzione di una tassellazione Voronoi di quei punti. Le coordinate delle tessere Voronoi che rientrano nel poligono della feature lineare cadono sulla linea centrale. Trasforma questi punti in una linea.

La tessellation di Voronoi viene eseguita in modo veramente efficiente in R utilizzando il pacchetto deldir e le intersezioni di poligoni e punti con il pacchetto 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"))) 

I SpatialLines risultanti oggetto: