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

Кроме QGIS
Ответить
Аватара пользователя
rhot
Гуру
Сообщения: 1727
Зарегистрирован: 25 янв 2011, 17:50
Репутация: 194
Ваше звание: доктор
Откуда: Архангельск

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

Сообщение rhot » 21 июл 2015, 13:06

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

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

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

--
Из того, что нашёл по теме:
а) http://gis.stackexchange.com/questions/ ... 3587#63587 - не подходит, так как территория очень большая - лучше в растре;
б) r.resamp.filter - не понятно как работает;
в) в R spatstat::density - тоже не понятно как работает. Сомневаюсь, что вообще подходит.
___________(¯`·.¸(¯`·.¸ Scientia potentia est _/ {SILVA}:::{FOSS}:::{GIS} \_ Знание сила ¸.·´¯)¸.·´¯)___________

trir
Гуру
Сообщения: 5278
Зарегистрирован: 09 апр 2010, 19:30
Репутация: 1014
Ваше звание: просто мимо прохожу
Откуда: Ё-бург

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

Сообщение trir » 21 июл 2015, 13:13

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

Аватара пользователя
rhot
Гуру
Сообщения: 1727
Зарегистрирован: 25 янв 2011, 17:50
Репутация: 194
Ваше звание: доктор
Откуда: Архангельск

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

Сообщение rhot » 21 июл 2015, 13:18

trir писал(а):1. Разрезать линни по областям
2. Получить суммарную длину линий для каждой области
3. Разделить сумарную длину на площадь области
Это пункт (а) из того, что я нашёл по теме. Не подходит - очень ресурсозатратно.
___________(¯`·.¸(¯`·.¸ Scientia potentia est _/ {SILVA}:::{FOSS}:::{GIS} \_ Знание сила ¸.·´¯)¸.·´¯)___________

Александр Мурый
Гуру
Сообщения: 5173
Зарегистрирован: 26 сен 2009, 16:26
Репутация: 793
Ваше звание: званий не имею
Откуда: Москва

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

Сообщение Александр Мурый » 21 июл 2015, 13:24

<r.resamp.filter> здесь точно ни при чём.
Редактор материалов, модератор форума

Александр Мурый
Гуру
Сообщения: 5173
Зарегистрирован: 26 сен 2009, 16:26
Репутация: 793
Ваше звание: званий не имею
Откуда: Москва

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

Сообщение Александр Мурый » 21 июл 2015, 13:34

Может, этот пример как-то поможет?
Редактор материалов, модератор форума

Аватара пользователя
rhot
Гуру
Сообщения: 1727
Зарегистрирован: 25 янв 2011, 17:50
Репутация: 194
Ваше звание: доктор
Откуда: Архангельск

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

Сообщение rhot » 21 июл 2015, 18:23

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

Пробывал сделать UPDATE в одну строчку вместе с WITH, но у меня GRASS собран с SQLite 3.7.13 и без поддержки PostgreSQL. :( Может можно как то ещё 2 запроса в один уместить. В шелле не хочу делать, потому что этот рассчёт плотности сети является частью R-скрипта.
___________(¯`·.¸(¯`·.¸ Scientia potentia est _/ {SILVA}:::{FOSS}:::{GIS} \_ Знание сила ¸.·´¯)¸.·´¯)___________

gamm
Гуру
Сообщения: 4046
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 1050
Ваше звание: программист
Откуда: Казань

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

Сообщение gamm » 21 июл 2015, 18:41

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

Александр Мурый
Гуру
Сообщения: 5173
Зарегистрирован: 26 сен 2009, 16:26
Репутация: 793
Ваше звание: званий не имею
Откуда: Москва

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

Сообщение Александр Мурый » 22 июл 2015, 14:34

Ещё вариант в 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
Снимок.png (57.77 КБ) 7910 просмотров

При необходимости скрипт трансформируется для использования в 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
Гуру
Сообщения: 1727
Зарегистрирован: 25 янв 2011, 17:50
Репутация: 194
Ваше звание: доктор
Откуда: Архангельск

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

Сообщение rhot » 23 июл 2015, 09:26

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

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

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. :D

По мотивам [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.
Последний раз редактировалось rhot 23 июл 2015, 09:28, всего редактировалось 1 раз.
___________(¯`·.¸(¯`·.¸ Scientia potentia est _/ {SILVA}:::{FOSS}:::{GIS} \_ Знание сила ¸.·´¯)¸.·´¯)___________

Ответить

Вернуться в «Свободные, бесплатные, открытые ГИС»

Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и 4 гостя