У меня есть большой (~ 70 МБ) шейп-файл дорог, и я хочу преобразовать его в растр с плотностью дорог в каждой ячейке. В идеале я хотел бы сделать это в R вместе с инструментами командной строки GDAL, если это необходимо.
Мой первоначальный подход состоял в том, чтобы напрямую рассчитать длину отрезков в каждой ячейке в соответствии с этим потоком . Это дает желаемые результаты, но довольно медленно даже для шейп-файлов, намного меньших, чем у меня. Вот очень упрощенный пример, для которого правильные значения ячеек очевидны:
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
Выглядит хорошо, но не масштабируемо! В паре других вопросов spatstat::density.psp()
функция была рекомендована для этой задачи. Эта функция использует подход плотности ядра. Я могу реализовать это, и это кажется быстрее, чем вышеупомянутый подход, но мне неясно, как выбрать параметры или интерпретировать результаты. Вот приведенный выше пример с использованием 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
Я подумал, что в этом случае подход ядра вычисляет плотность, а не длину на ячейку, поэтому я преобразовал:
# 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
Но ни в том, ни в другом случае подход к ядру не приближается к более прямому подходу, описанному выше.
Итак, мои вопросы:
- Как я могу интерпретировать вывод
density.psp
функции? Какие единицы? - Как выбрать
sigma
параметрdensity.psp
таким образом, чтобы результаты соответствовали более прямому, интуитивному подходу, описанному выше? - Бонус: что на самом деле делает плотность линий ядра? У меня есть некоторый смысл для того, как эти подходы работают для точек, но я не вижу, как это распространяется на линии.