Ich habe ein großes (~ 70 MB) Shapefile von Straßen und möchte dieses in ein Raster mit Straßendichte in jeder Zelle konvertieren. Idealerweise möchte ich dies in R zusammen mit GDAL-Befehlszeilentools tun, falls erforderlich.
Mein anfänglicher Ansatz bestand darin, die Länge der Liniensegmente in jeder Zelle anhand dieses Threads direkt zu berechnen . Dies führt zu den gewünschten Ergebnissen, ist jedoch selbst für Shapefiles, die viel kleiner als meine sind, recht langsam. Hier ist ein sehr vereinfachtes Beispiel, für das die korrekten Zellenwerte offensichtlich sind:
require(sp)
require(raster)
require(rgeos)
require(RColorBrewer)
# Create some sample lines
l1 <- Lines(Line(cbind(c(0,1),c(.25,0.25))), ID="a")
l2 <- Lines(Line(cbind(c(0.25,0.25),c(0,1))), ID="b")
sl <- SpatialLines(list(l1,l2))
# Function to calculate lengths of lines in given raster cell
lengthInCell <- function(i, r, l) {
r[i] <- 1
rpoly <- rasterToPolygons(r, na.rm=T)
lc <- crop(l, rpoly)
if (!is.null(lc)) {
return(gLength(lc))
} else {
return(0)
}
}
# Make template
rLength <- raster(extent(sl), res=0.5)
# Calculate lengths
lengths <- sapply(1:ncell(rLength), lengthInCell, rLength, sl)
rLength[] <- lengths
# Plot results
spplot(rLength, scales = list(draw=TRUE), xlab="x", ylab="y",
col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")),
sp.layout=list("sp.lines", sl),
par.settings=list(fontsize=list(text=15)))
round(as.matrix(rLength),3)
#### Results
[,1] [,2]
[1,] 0.5 0.0
[2,] 1.0 0.5
Sieht gut aus, ist aber nicht skalierbar! In einigen anderen Fragen wurde die spatstat::density.psp()
Funktion für diese Aufgabe empfohlen. Diese Funktion verwendet einen Kerneldichteansatz. Ich kann es implementieren und es scheint schneller zu sein als der obige Ansatz, aber ich weiß nicht, wie ich die Parameter auswählen oder die Ergebnisse interpretieren soll. Hier ist das obige Beispiel mit density.psp()
:
require(spatstat)
require(maptools)
# Convert SpatialLines to psp object using maptools library
pspSl <- as.psp(sl)
# Kernel density, sigma chosen more or less arbitrarily
d <- density(pspSl, sigma=0.01, eps=0.5)
# Convert to raster
rKernDensity <- raster(d)
# Values:
round(as.matrix(rKernDensity),3)
#### Results
[,1] [,2]
[1,] 0.100 0.0
[2,] 0.201 0.1
Ich dachte, es könnte der Fall sein, dass der Kernel-Ansatz die Dichte im Gegensatz zur Länge pro Zelle berechnet, also habe ich konvertiert:
# Convert from density to length per cell for comparison
rKernLength <- rKernDensity * res(rKernDensity)[1] * res(rKernDensity)[2]
round(as.matrix(rKernLength),3)
#### Results
[,1] [,2]
[1,] 0.025 0.000
[2,] 0.050 0.025
In keinem Fall kommt der Kernel-Ansatz jedoch einer Ausrichtung auf den oben beschriebenen direkteren Ansatz nahe.
Meine Fragen lauten also:
- Wie kann ich die Ausgabe der
density.psp
Funktion interpretieren ? Was sind die Einheiten? - Wie kann ich den
sigma
Parameterdensity.psp
so auswählen, dass die Ergebnisse mit dem direkteren, intuitiveren Ansatz oben übereinstimmen? - Bonus: Was macht eigentlich die Kernel-Liniendichte? Ich habe ein Gespür dafür, wie diese Ansätze für Punkte funktionieren, sehe aber nicht, wie sich das auf Linien erstreckt.