Obliczanie gęstości drogi w R za pomocą gęstości jądra? [Zamknięte]

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:

1. Jak mogę zinterpretować wyjście density.pspfunkcji? Jakie są jednostki?
2. Jak mogę wybrać sigmaparametr, density.pspaby wyniki były zgodne z bardziej bezpośrednim, intuicyjnym podejściem powyżej?
3. 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.

Zadałem to pytanie na listserv R-sig-Geo i otrzymałem pomocną odpowiedź od Adriana Baddeleya, jednego z autorów spatstatów . Tutaj przedstawię moją interpretację jego odpowiedzi na potomstwo.

Adrian zauważa, że ​​ta funkcja spatstat::pixellate.psp()lepiej pasuje do mojego zadania. Ta funkcja konwertuje wzór segmentu linii (lub SpatialLinesobiekt z konwersją) na obraz pikselowy (lub RasterLayerz konwersją), gdzie wartość w każdej komórce jest długością segmentów linii przechodzących przez tę komórkę. Właśnie tego szukam!

Rozdzielczość wynikowego obrazu można zdefiniować za pomocą epsparametru lub dimyxparametru, który określa wymiary (liczbę wierszy i kolumn).

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

Wyniki są dokładnie takie, jak pożądane.

Adrian również odpowiedział na moje pytania dotyczące spatstat::density.psp(). Wyjaśnia, że ​​ta funkcja:

oblicza splot jądra Gaussa z liniami. Intuicyjnie oznacza to, że density.psp„rozmazuje” linie w dwuwymiarową przestrzeń. To density(L)jest jak rozmyta wersja pixellate(L). W rzeczywistości density(L)jest bardzo podobny do tego, blur(pixellate(L))gdzie blurjest inna spatstatfunkcja, która rozmywa obraz. [Parametr] sigmato przepustowość jądra Gaussa. Wartość dla density.psp(L)danego piksela u jest czymś w rodzaju całkowitej długości linii w okręgu sigma promienia wokół piksela u, z tym wyjątkiem, że jest to naprawdę ważona średnia takich udziałów z różnych promieni okręgu. Jednostki mają długość ^ (- 1), tj. Długość linii na jednostkę powierzchni.

Nie jest dla mnie jasne, kiedy metoda jądra Gaussa density.psp()byłaby lepsza niż bardziej intuicyjne podejście do bezpośredniego obliczania długości linii pixellate(). Chyba będę musiał zostawić to ekspertom.