Форумы GIS-Lab.info

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

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

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

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

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

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

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

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

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

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

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

.

Пробывал сделать UPDATE в одну строчку вместе с WITH, но у меня GRASS собран с SQLite 3.7.13 и без поддержки PostgreSQL.

Может можно как то ещё 2 запроса в один уместить. В шелле не хочу делать, потому что этот рассчёт плотности сети является частью R-скрипта.

не мудрите. Наставьте погуще точек вдоль линий (скриптом в 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 КБ) 8775 просмотров

При необходимости скрипт трансформируется для использования в 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)

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

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

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.

Форумы GIS-Lab.info

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

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

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

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

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

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

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

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

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

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