Плотность дорожной сети

[rhot]

Есть вектор дорог.
Требуется вычислить плотность дорожной сети (в км на км^2).

В ArcGIS всё решается одной специальной командой - Плотность линий.

1) Есть ли аналоги этой команды в свободных ГИС?
2) Не могу разобраться с опцией Численность в инструменте Плотность линий. Объясните на примере, что она означает.

--
Из того, что нашёл по теме:
а) http://gis.stackexchange.com/questions/ ... 3587#63587 - не подходит, так как территория очень большая - лучше в растре;
б) r.resamp.filter - не понятно как работает;
в) в R spatstat::density - тоже не понятно как работает. Сомневаюсь, что вообще подходит.

[trir]

1. Разрезать линни по областям
2. Получить суммарную длину линий для каждой области
3. Разделить сумарную длину на площадь области

[rhot]

1. Разрезать линни по областям
2. Получить суммарную длину линий для каждой области
3. Разделить сумарную длину на площадь области

Это пункт (а) из того, что я нашёл по теме. Не подходит - очень ресурсозатратно.

[Александр Мурый]

<r.resamp.filter> здесь точно ни при чём.

[Александр Мурый]

Может, этот пример как-то поможет?

[rhot]

Делал всё по инструкции, остановился на

Тут, скорее всего, можно на чистом SQL, но я сделал в виде shell-скрипта

.

Пробывал сделать UPDATE в одну строчку вместе с WITH, но у меня GRASS собран с SQLite 3.7.13 и без поддержки PostgreSQL. Может можно как то ещё 2 запроса в один уместить. В шелле не хочу делать, потому что этот рассчёт плотности сети является частью R-скрипта.

[gamm]

не мудрите. Наставьте погуще точек вдоль линий (скриптом в R), в атрибуты запишите среднюю длину, приходящуюся на точку (длина линии/(n-1)), а потом просуммируйте точки по полигонам (в том же R, там есть point_in_poly()) - присвойте точке ID полигона, и аггрегируйте

[Александр Мурый]

Ещё вариант в R, подсмотренный здесь:

Код: Выделить всё

require(sp)
require(raster)
require(spatstat)
require(maptools)
#require(RColorBrewer)

# Load data
lines <- shapefile('lines.shp')

# Convert to a line segment pattern object with maptools
linesPSP <- as.psp(as(lines, 'SpatialLines'))

# Calculate lengths per cell (resolution ~ 1000m)
linesLengthIM <- pixellate.psp(linesPSP, eps=1000)

# Convert pixel image to raster
linesLength <- raster(linesLengthIM, crs=projection(lines))

# Convert raster to polygons (if needed)
poly <- rasterToPolygons(linesLength)

# spplot(linesLength, scales = list(draw=TRUE), xlab="x", ylab="y", col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")), sp.layout=list("sp.lines", lines), par.settings=list(fontsize=list(text=15)))

В итоге получается растр с разрешением около 1км в СК исходного вектора и полигоны со значениями плотности (если они вообще нужны). Через spplot можно всё это изобразить.

Снимок.png (57.77 КБ) 8922 просмотра

При необходимости скрипт трансформируется для использования в QGIS:

Код: Выделить всё

##[Own Scripts] = group
##Vect = vector
##lines_density_rast = output raster
##lines_density_vect = output vector

library(sp)
library(raster)
library(spatstat)
library(maptools)

roadsPSP <- as.psp(as(Vect, 'SpatialLines'))
roadLengthIM <- pixellate.psp(roadsPSP, eps=1000)
roadLength <- raster(roadLengthIM, crs=projection(Vect))
lines_density_rast <- roadLength
lines_density_vect <- rasterToPolygons(lines_density_rast)

[rhot]

Пробовал вышеописанный метод, но не разобрался по-хорошему и бросил

Код: Выделить всё

library(maptools)
library(statspat)
win <- square(1000)
spatstat.options(checksegments=F) # some segments do not lie entirely inside the window, this option allows to omit the error
pspSl <- as.psp.SpatialLinesDataFrame(roads.bel) # coerce to psp class
# Pixellate with resolution of 1000, i.e. polygon 1000x1000 pixels
px <- pixellate(pspSl, W=win)
rLength <- raster(px)
plot(rLength)

Проблема в том, что я выбирал опцию 'W', а не 'ops'. И растр у меня получался размером 1000х1000.

По мотивам [1] сосчитал плотность дорог через 'kernel density' в GRASS. Радиус функции выбрал 500 м, т.к. этот рассчёт плотности дорог нужен для оценки лесохозяйственного потенциала, а форвардер обычно дальше полкилометра в лес не ездит - не экономично.

Код: Выделить всё

execGRASS('v.to.points', input='roads', type='line', output='p_roads', dmax=200, flags='overwrite')
execGRASS('g.region', vector='bounds', res='200', flags='a')
execGRASS('g.remove', type='raster', name='roads.ker500', flags='f')
execGRASS('v.kernel', input='p_roads', output='roads.ker500', radius=500)
execGRASS('r.mapcalc', expression='roads.ker500=20000*roads.ker500', flags='overwrite')
execGRASS('g.region', res='1000', flags='a')
execGRASS('r.resamp.stats', input='roads.ker500', output='roads.ker500.1km', method='sum', flags=c('w','n','overwrite'))

--
1.) Cai, X., Wu, Z. & Cheng, J. (2013). Using kernel density estimation to assess the spatial pattern of road density and its impact on landscape fragmentation. International Journal of Geographical Information Science, 27(2), pp 222–230.