Calcul de la densité de la route dans R en utilisant la densité du noyau? [fermé]

J'ai un grand fichier de formes (~ 70 Mo) de routes et je souhaite le convertir en un raster avec une densité de route dans chaque cellule. Idéalement, je voudrais le faire en R avec les outils de ligne de commande GDAL si nécessaire.

Mon approche initiale était de calculer directement les longueurs des segments de ligne dans chaque cellule selon ce fil . Cela produit les résultats souhaités, mais est assez lent, même pour des fichiers de formes beaucoup plus petits que le mien. Voici un exemple très simplifié pour lequel les valeurs de cellule correctes sont évidentes:

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

Semble bon, mais pas évolutif! Dans quelques autres questions, la spatstat::density.psp()fonction a été recommandée pour cette tâche. Cette fonction utilise une approche de densité du noyau. Je suis capable de l'implémenter et cela semble plus rapide que l'approche ci-dessus, mais je ne sais pas comment choisir les paramètres ou interpréter les résultats. Voici l'exemple ci-dessus en utilisant 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

J'ai pensé que cela pourrait être le cas si l'approche du noyau calcule la densité plutôt que la longueur par cellule, j'ai donc 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

Mais, dans aucun des cas, l'approche du noyau ne se rapproche de l'approche plus directe ci-dessus.

Donc, mes questions sont:

1. Comment puis-je interpréter la sortie de la density.pspfonction? Quelles sont les unités?
2. Comment puis-je choisir le sigmaparamètre pour density.pspque les résultats s'alignent sur l'approche plus directe et intuitive ci-dessus?
3. Bonus: que fait réellement la densité de la ligne du noyau? J'ai une idée de la façon dont ces approches fonctionnent pour les points, mais je ne vois pas comment cela s'étend aux lignes.

J'ai posté cette question sur la liste de diffusion R-sig-Geo et j'ai reçu une réponse utile d'Adrian Baddeley, l'un des auteurs de spatstats . Je posterai ici mon interprétation de sa réponse pour la postérité.

Adrian note que la fonction spatstat::pixellate.psp()correspond mieux à ma tâche. Cette fonction convertit un motif de segment de ligne (ou un SpatialLinesobjet avec conversion) en une image pixel (ou RasterLayeravec conversion), où la valeur dans chaque cellule est la longueur des segments de ligne passant par cette cellule. Exactement ce que je recherche!

La résolution de l'image résultante peut être définie avec le epsparamètre ou le dimyxparamètre, qui définit les dimensions (nombre de lignes et de colonnes).

require(sp)
require(raster)
require(maptools)
require(spatstat)

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

# Convert SpatialLines to psp object using maptools library
pspSl <- as.psp(sl)
# Pixellate with resolution of 0.5, i.e. 2x2 pixels
px <- pixellate(pspSl, eps=0.5)
# This can be converted to raster as desired
rLength <- raster(px)
# Values:
round(as.matrix(rLength),3)

     [,1] [,2]
[1,]  0.5  0.0
[2,]  1.0  0.5

Les résultats sont exactement comme souhaité.

Adrian a également répondu à mes questions spatstat::density.psp(). Il explique que cette fonction:

calcule la convolution du noyau gaussien avec les lignes. Intuitivement, cela signifie que density.psp«étale» les lignes dans un espace à deux dimensions. C'est density(L)comme une version floue de pixellate(L). En fait, density(L)c'est très similaire à blur(pixellate(L))où blurest une autre spatstatfonction qui brouille une image. [Le paramètre] sigmaest la bande passante du noyau gaussien. La valeur de density.psp(L)à un pixel u donné est quelque chose comme la quantité totale de longueur de ligne dans un cercle de rayon sigma autour du pixel u, sauf que c'est vraiment une moyenne pondérée de ces contributions de rayons de cercle différents. Les unités sont la longueur ^ (- 1), c'est-à-dire la longueur de ligne par unité de surface.

Je ne sais pas density.psp()trop quand l'approche du noyau gaussien de serait préférée à l'approche plus intuitive du calcul direct des longueurs de ligne dans pixellate(). Je suppose que je devrai laisser cela aux experts.