Eu tenho um grande arquivo de forma (~ 70 MB) de estradas e quero convertê-lo em um raster com densidade de estrada em cada célula. Idealmente, eu gostaria de fazer isso no R junto com as ferramentas de linha de comando GDAL, se necessário.
Minha abordagem inicial foi calcular diretamente os comprimentos dos segmentos de linha em cada célula conforme esse segmento . Isso produz os resultados desejados, mas é bastante lento, mesmo para arquivos de forma muito menores que os meus. Aqui está um exemplo muito simplificado para o qual os valores corretos das células são óbvios:
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
Parece bom, mas não escalável! Em algumas outras questões, a spatstat::density.psp()
função foi recomendada para esta tarefa. Esta função usa uma abordagem de densidade do kernel. Eu sou capaz de implementá-lo e parece mais rápido que a abordagem acima, mas não sei como escolher os parâmetros ou interpretar os resultados. Aqui está o exemplo acima usando 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
Eu pensei que poderia ser o caso de a abordagem do kernel calcular a densidade em oposição ao comprimento por célula, então eu converti:
# 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
Mas, em nenhum dos casos, a abordagem do kernel chega perto de se alinhar com a abordagem mais direta acima.
Então, minhas perguntas são:
- Como posso interpretar a saída da
density.psp
função? Quais são as unidades? - Como posso escolher o
sigma
parâmetrodensity.psp
para que os resultados se alinhem à abordagem mais direta e intuitiva acima? - Bônus: o que a densidade de linha do kernel está realmente fazendo? Eu tenho algum senso de como essas abordagens funcionam para obter pontos, mas não vejo como isso se estende às linhas.