Tengo un gran archivo de formas (~ 70 MB) de carreteras y quiero convertir esto en una trama con densidad de carreteras en cada celda. Idealmente, me gustaría hacer esto en R junto con las herramientas de línea de comandos GDAL si es necesario.
Mi enfoque inicial fue calcular directamente las longitudes de los segmentos de línea en cada celda según este hilo . Esto produce los resultados deseados, pero es bastante lento incluso para archivos de forma mucho más pequeños que los míos. Aquí hay un ejemplo muy simplificado para el cual los valores de celda correctos son obvios:
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
¡Se ve bien, pero no escalable! En algunas otras preguntas, la spatstat::density.psp()
función se ha recomendado para esta tarea. Esta función utiliza un enfoque de densidad del núcleo. Puedo implementarlo y parece más rápido que el enfoque anterior, pero no tengo claro cómo elegir los parámetros o interpretar los resultados. Aquí está el ejemplo anterior 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
Pensé que podría ser el caso de que el enfoque del núcleo calcula la densidad en lugar de la longitud por celda, por lo que convertí:
# 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
Pero, en ninguno de los casos, el enfoque del núcleo se acerca a alinearse con el enfoque más directo anterior.
Entonces, mis preguntas son:
- ¿Cómo puedo interpretar la salida de la
density.psp
función? ¿Cuáles son las unidades? - ¿Cómo puedo elegir el
sigma
parámetrodensity.psp
para que los resultados se alineen con el enfoque más directo e intuitivo anterior? - Bonificación: ¿qué está haciendo realmente la densidad de línea del núcleo? Tengo algún sentido sobre cómo funcionan estos enfoques para los puntos, pero no veo cómo se extiende a las líneas.