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

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.

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

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

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

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

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 

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

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.