Mam duży (~ 70 MB) plik kształtów dróg i chcę przekonwertować go na raster z gęstością dróg w każdej komórce. Idealnie chciałbym to zrobić w R wraz z narzędziami wiersza poleceń GDAL, jeśli to konieczne.
Moje początkowe podejście polegało na bezpośrednim obliczeniu długości odcinków linii w każdej komórce zgodnie z tym wątkiem . Daje to pożądane rezultaty, ale jest dość powolne, nawet dla plików kształtów znacznie mniejszych niż moje. Oto bardzo uproszczony przykład, dla którego prawidłowe wartości komórek są oczywiste:
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
Wygląda dobrze, ale nie można go skalować! W kilku innych pytaniach spatstat::density.psp()
funkcja została zalecona do tego zadania. Ta funkcja wykorzystuje podejście gęstości jądra. Jestem w stanie go wdrożyć i wydaje się to szybsze niż powyższe podejście, ale nie jestem pewien, jak wybrać parametry lub zinterpretować wyniki. Oto powyższy przykład, używając 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
Pomyślałem, że może być tak, że podejście jądra oblicza gęstość w przeciwieństwie do długości na komórkę, więc przekonwertowałem:
# 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
Ale w żadnym przypadku podejście do jądra nie jest zbliżone do bardziej bezpośredniego podejścia opisanego powyżej.
Tak więc moje pytania to:
- Jak mogę zinterpretować wyjście
density.psp
funkcji? Jakie są jednostki? - Jak mogę wybrać
sigma
parametr,density.psp
aby wyniki były zgodne z bardziej bezpośrednim, intuicyjnym podejściem powyżej? - Bonus: co tak naprawdę robi gęstość linii jądra? Mam pewien sens, jak te podejścia działają w przypadku punktów, ale nie widzę, jak to się rozciąga na linie.