Linienformdatei in Raster konvertieren, Wert = Gesamtlänge der Linien in der Zelle

14

Ich habe ein Linien-Shapefile, das ein Straßennetz darstellt. Ich möchte diese Daten rastern, wobei die resultierenden Werte im Raster die Gesamtlänge der Linien anzeigen, die in die Rasterzelle fallen.

Die Daten befinden sich in der British National Grid-Projektion, sodass die Einheiten Meter sind.

Idealerweise möchte ich diese Operation mit ausführen R, und ich vermute, dass die rasterizeFunktion aus dem rasterPaket dazu beitragen würde. Ich kann einfach nicht herausfinden, welche Funktion angewendet werden soll.

JPD
quelle
Vielleicht vignette('over', package = 'sp')könnte helfen.

Antworten:

12

Nach einer kürzlich gestellten Frage möchten Sie möglicherweise die Funktionen des rgeos- Pakets nutzen, um Ihr Problem zu lösen. Aus Gründen der Reproduzierbarkeit habe ich ein Shapefile der tansanischen Straßen von DIVA-GIS heruntergeladen und in mein aktuelles Arbeitsverzeichnis gestellt. Für die anstehenden Aufgaben benötigen Sie drei Pakete:

  • rgdal für den allgemeinen Umgang mit Geodaten
  • Raster zur Rasterung der Shapefile-Daten
  • rgeos , um die Kreuzung von Straßen mit einer Rastervorlage zu überprüfen und die Straßenlängen zu berechnen

Folglich sollten Ihre ersten Zeilen von could so aussehen:

library(rgdal)
library(raster)
library(rgeos)

Danach müssen Sie die Shapefile-Daten importieren. Beachten Sie, dass DIVA-GIS-Shapefiles in EPSG: 4326 verteilt sind. Daher projiziere ich das Shapefile in EPSG: 21037 (UTM 37S), um Meter statt Grad zu verarbeiten.

roads <- readOGR(dsn = ".", layer = "TZA_roads")
roads_utm <- spTransform(roads, CRS("+init=epsg:21037"))

Für die anschließende Rasterung benötigen Sie eine Rastervorlage, die die räumliche Ausdehnung Ihres Shapefiles abdeckt. Die Rastervorlage besteht standardmäßig aus 10 Zeilen und 10 Spalten, wodurch zu lange Rechenzeiten vermieden werden.

roads_utm_rst <- raster(extent(roads_utm), crs = projection(roads_utm))

Nachdem Sie die Vorlage eingerichtet haben, durchlaufen Sie alle Zellen des Rasters (das derzeit nur aus NA-Werten besteht). Durch Zuweisen eines Werts von '1' zu der aktuellen Zelle und anschließende Ausführung rasterToPolygonswird im resultierenden Shapefile 'tmp_shp' automatisch die Ausdehnung des aktuell verarbeiteten Pixels gespeichert. gIntersectsErkennt, ob sich diese Ausdehnung mit Straßen überschneidet. Wenn nicht, gibt die Funktion den Wert '0' zurück. Andernfalls wird das Straßen-Shapefile von der aktuellen Zelle abgeschnitten und die Gesamtlänge von "SpatialLines" in dieser Zelle wird mit berechnet gLength.

lengths <- sapply(1:ncell(roads_utm_rst), function(i) {
  tmp_rst <- roads_utm_rst
  tmp_rst[i] <- 1
  tmp_shp <- rasterToPolygons(tmp_rst)

  if (gIntersects(roads_utm, tmp_shp)) {
    roads_utm_crp <- crop(roads_utm, tmp_shp)
    roads_utm_crp_length <- gLength(roads_utm_crp)
    return(roads_utm_crp_length)
  } else {
    return(0)
  }
})

Zuletzt können Sie die berechneten Längen (die in Kilometer umgerechnet werden) in die Rastervorlage einfügen und Ihre Ergebnisse visuell überprüfen.

roads_utm_rst[] <- lengths / 1000

library(RColorBrewer)
spplot(roads_utm_rst, scales = list(draw = TRUE), xlab = "x", ylab = "y", 
       col.regions = colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout = list("sp.lines", roads_utm), 
       par.settings = list(fontsize = list(text = 15)), at = seq(0, 1800, 200))

lines2raster

fdetsch
quelle
Das ist fantastisch! Ich habe sapply()auf pbsapply()das Cluster-Argument gewechselt und es verwendet cl = detectCores()-1. Jetzt kann ich dieses Beispiel parallel ausführen!
philiporlando
11

Das Folgende ist aus der Lösung von Jeffrey Evans modifiziert. Diese Lösung ist viel schneller, da keine Rasterung verwendet wird

library(raster)
library(rgdal)
library(rgeos)

roads <- shapefile("TZA_roads.shp")
roads <- spTransform(roads, CRS("+proj=utm +zone=37 +south +datum=WGS84"))
rs <- raster(extent(roads), crs=projection(roads))
rs[] <- 1:ncell(rs)

# Intersect lines with raster "polygons" and add length to new lines segments
rsp <- rasterToPolygons(rs)
rp <- intersect(roads, rsp)
rp$length <- gLength(rp, byid=TRUE) / 1000
x <- tapply(rp$length, rp$layer, sum)
r <- raster(rs)
r[as.integer(names(x))] <- x
Robert Hijmans
quelle
Scheint die effizienteste und eleganteste der aufgelisteten Methoden zu sein. Außerdem, hatte ich raster::intersect() vorher nicht gesehen , gefällt mir, dass es im Gegensatz zu den Attributen der geschnittenen Features kombiniert rgeos::gIntersection().
Matt SM
+1 immer schön, effizientere Lösungen zu sehen!
Jeffrey Evans
@RobertH, ich habe versucht, diese Lösung für ein anderes Problem zu verwenden, bei dem ich das Gleiche tun möchte, wie in diesem Thread, aber mit einem sehr großen Shapefile von Straßen (für einen gesamten Kontinent). Es scheint geklappt zu haben, aber wenn ich versuche, die Figur wie von @ fdetsch zu machen, erhalte ich kein zusammenhängendes Gitter, sondern nur ein paar farbige Quadrate im Bild (siehe tinypic.com/r/20hu87k/9 ).
Doon_Bogan
Und am effizientesten ... mit meinem Beispieldatensatz: Laufzeit 0,6 Sekunden für diese Lösung im Vergleich zu 8,25 Sekunden für die Lösung mit den meisten positiven Stimmen.
user3386170
1

Sie benötigen keine for-Schleife. Schneiden Sie einfach alles auf einmal und fügen Sie dann mit der Funktion "SpatialLinesLengths" in sp den neuen Liniensegmenten Linienlängen hinzu. Mithilfe der Raster-Package-Rasterize-Funktion mit dem Argument fun = sum können Sie dann ein Raster mit der Summe der Zeilenlängen erstellen, die jede Zelle schneiden. Unter Verwendung der obigen Antwort und der zugehörigen Daten handelt es sich hier um Code, der die gleichen Ergebnisse generiert.

require(rgdal)
require(raster)
require(sp)
require(rgeos)

setwd("D:/TEST/RDSUM")
roads <- readOGR(getwd(), "TZA_roads")
  roads <- spTransform(roads, CRS("+init=epsg:21037"))
    rrst <- raster(extent(roads), crs=projection(roads))

# Intersect lines with raster "polygons" and add length to new lines segments
rrst.poly <- rasterToPolygons(rrst)
  rp <- gIntersection(roads, rrst.poly, byid=TRUE)
    rp <- SpatialLinesDataFrame(rp, data.frame(row.names=sapply(slot(rp, "lines"), 
                               function(x) slot(x, "ID")), ID=1:length(rp), 
                               length=SpatialLinesLengths(rp)/1000) ) 

# Rasterize using sum of intersected lines                            
rd.rst <- rasterize(rp, rrst, field="length", fun="sum")

# Plot results
require(RColorBrewer)
spplot(rd.rst, scales = list(draw=TRUE), xlab="x", ylab="y", 
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout=list("sp.lines", rp), 
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))
Jeffrey Evans
quelle
Zum ersten Mal sehe ich SpatialLinesLengths. Schätze, es ist nie zu spät zu lernen, danke (: rasterizedauert allerdings ziemlich lange (7-mal länger als der obere Ansatz auf meinem Computer).
fdetsch
Mir ist aufgefallen, dass das Rastern langsam ist. Bei großen Problemen verlangsamt eine for-Schleife jedoch die Geschwindigkeit und ich denke, dass Sie mit dem Rasterisieren eine viel optimiertere Lösung finden werden. Außerdem ist der Entwickler des Raster-Pakets bestrebt, jedes Release optimierter und schneller zu machen.
Jeffrey Evans
Ein mögliches Problem, das ich bei dieser Technik festgestellt habe, ist, dass die rasterize()Funktion alle Linien enthält, die eine bestimmte Zelle berühren . Dies führt in einigen Fällen dazu, dass Liniensegmentlängen zweimal gezählt werden: einmal in der Zelle, die sie haben sollen, und einmal in der angrenzenden Zelle, die der Linienendpunkt gerade berührt.
Matt SM
Ja, aber "rp" ist das Objekt, das gerastert wird und der Schnittpunkt von Polygonen und Punkten ist.
Jeffrey Evans
1

Hier ist noch ein anderer Ansatz. Es weicht von den bereits mit dem spatstatPaket gegebenen ab . Soweit ich das beurteilen kann, verfügt dieses Paket über eine eigene Version von räumlichen Objekten (z. B. im imVergleich zu rasterObjekten), aber das maptoolsPaket ermöglicht die Konvertierung zwischen spatstatObjekten und räumlichen Standardobjekten.

Dieser Ansatz stammt aus diesem R-sig-Geo-Beitrag .

require(sp)
require(raster)
require(rgdal)
require(spatstat)
require(maptools)
require(RColorBrewer)

# Load data and transform to UTM
roads <- shapefile('data/TZA_roads.shp')
roadsUTM <- spTransform(roads, CRS("+init=epsg:21037"))

# Need to convert to a line segment pattern object with maptools
roadsPSP <- as.psp(as(roadsUTM, 'SpatialLines'))

# Calculate lengths per cell
roadLengthIM <- pixellate.psp(roadsUTM, dimyx=10)

# Convert pixel image to raster in km
roadLength <- raster(dtanz / 1000, crs=projection(roadsUTM))

# Plot
spplot(rtanz, scales = list(draw=TRUE), xlab="x", ylab="y", 
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout=list("sp.lines", roadsUTM), 
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))

Imgur

Das langsamste Bit konvertiert die Straßen von SpatialLinesin ein Liniensegmentmuster (dh spatstat::psp). Sobald dies erledigt ist, sind die tatsächlichen Längenberechnungsteile ziemlich schnell, sogar für viel höhere Auflösungen. Zum Beispiel auf meinem alten 2009 MacBook:

system.time(pixellate(tanzpsp, dimyx=10))
#    user  system elapsed 
#   0.007   0.001   0.007

system.time(pixellate(tanzpsp, dimyx=1000))
#    user  system elapsed 
#   0.146   0.032   0.178 
Matt SM
quelle
Hmmmm ... ich wünschte, diese Äxte wären nicht in wissenschaftlicher Notation. Hat jemand eine Idee, wie man das behebt?
Matt SM
Sie können die globalen R-Einstellungen ändern und die Notation mit: options (scipen = 999) deaktivieren. Ich weiß jedoch nicht, ob das Gitter die globalen Umgebungseinstellungen berücksichtigt.
Jeffrey Evans
0

Lassen Sie mich Ihnen die Paket- Ader mit mehreren Funktionen zum Bearbeiten von räumlichen Linien und zum Importieren von sf und data.table vorstellen

library(vein)
library(sf)
library(cptcity)
data(net)
netsf <- st_as_sf(net) #Convert Spatial to sf
netsf <- st_transform(netsf, 31983) # Project data
netsf$length_m  <- st_length(netsf)
netsf <- netsf[, "length_m"]
g <- make_grid(netsf, width = 1000) #Creat grid of 1000m spacing with columns id for each feature
# Number of lon points: 12
# Number of lat points: 11

gnet <- emis_grid(netsf, g)
plot(gnet["length_m"])

sf_to_raster <- function(x, column, ncol, nrow){
  x <- sf::as_Spatial(x)
  r <- raster::raster(ncol = ncol, nrow = nrow)
  raster::extent(r) <- raster::extent(x)
  r <- raster::rasterize(x, r, column)
  return(r)
}

rr <- sf_to_raster(gnet, "length_m", 12, 11)
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = cpt(), scales = list(draw = T))
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = cpt(pal = 5176), scales = list(draw = T))
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = lucky(), scales = list(draw = T))
# Colour gradient: neota_flor_apple_green, number: 6165

Bildbeschreibung hier eingeben

Bildbeschreibung hier eingeben

Bildbeschreibung hier eingeben

Sergio
quelle
-1

Dies mag ein bisschen naiv klingen, aber wenn es sich um ein Straßensystem handelt, wählen Sie die Straßen aus und speichern Sie sie in einer Zwischenablage. Suchen Sie dann ein Tool, mit dem Sie der Zwischenablage einen Puffer hinzufügen können. Stellen Sie ihn auf die zulässige Straßenbreite ein, dh 3 Meter +/- Denken Sie daran, dass sich der Puffer für jede Seite von der Mittellinie bis zur Kante * 2 i befindet, sodass ein 3-Meter-Puffer tatsächlich eine 6-Meter-Straße von einer Seite zur anderen ist.

derekB
quelle
Was hat die Straßenbreite mit der Straßenlänge zu tun? Diese Antwort geht nicht auf die Frage ein.
Alphabetasoup