Базовая оценка транспортной доступности средствами GRASS GIS и QGIS

Эта страница является черновиком статьи.

Построение карт транспортной доступности на основе данных OpenStreetMap средствами открытых геоинформационных систем GRASS GIS и QGIS

Картографирование транспортной доступности на основе данных о дорожной сети - одна из классических задач ГИС. Наиболее распространенным способом моделирования транспортной доступности является построение изохрон - линий равных затрат времени на преодоление пространства относительно заданных точек. В представленной статье обсуждается алгоритм построения изохрон по данным OpenStreetMap с использованием открытых ГИС GRASS GIS и QGIS. В QGIS будет осуществляться подготовка данных и картографическое представление результатов, а в GRASS собственно моделирование. Всю работу можно выполнить целиком в GRASS, но, по мнению автора, общие манипуляции геоданными и представление картографических материалов удачнее и удобнее реализованы в QGIS.

Описываемые в статье действия выполнялись в средах Ubuntu 14.04 LTS и Windows 8.1 (x64), GRASS 7.2, QGIS 2.14.

Получение и подготовка данных

Для осуществления расчётов нам потребуется набор векторных линейных геоданных, содержащий информацию о дорожно-транспортной сети исследуемой территории. Заполучить подобную информацию можно различными способами: приобрести у специализированных поставщиков, найти в одном из источников открытых данных, оцифровать атлас автомобильных дорог, нарисовать по космическому снимку, и так далее (не забывайте про лицензии данных!). В данном случае мы воспользуемся данными OSM как достаточно качественными и подробными, и, что важно, доступными бесплатно и легально. Начать знакомство с OSM вы можете с этой страницы. Этап подготовки данных являются ключевым: помимо технических этапов, здесь нам предстоит определить транспортные характеристики, которые будут использованы для моделирования.

Для демонстрации принципов базового транспортного моделирования рассмотрим Тосненский район Ленинградской области. Транспорт: автомобиль.

Загрузка данных OSM

Для загрузки данных OSM существует множество возможностей. Одним из наиболее простых способов является загрузка уже подготовленных наборов данных по слоям в формате ESRI Shapefile. Популярные сервисы в сети: Geofabrik, где вы можете найти комплекты данных на весь мир, и NextGIS, где вы найдёте вместе с данными оформленные QGIS-проекты по субъектам РФ и странам СНГ. Также удобный способ быстро получить данные OSM - воспользоваться одним из соответствующих плагинов для QGIS, например, OSMDownloader. При его установке в интерфейс QGIS добавится кнопка выделения прямоугольной области, для которой и будут загружены данные в указанное вами расположение. Для демонстрации в этой статье используется набор данных OSM от NextGIS на территорию Ленинградской области от 13 ноября 2016 года, доступный бесплатно (загрузить).

Для нашей задачи понадобятся только данные о дорожной сети (в случае работы с выгрузками NextGIS - слой highway-line из директории data; в случае работы с OSMDownloader - все загруженные объекты линейного типа, у которых в атрибуте "Highway" содержится какое-либо значение). Для оформительских задач также может быть использован набор геоданных с границами административно-территориального деления (например, для обрезки изохрон по границам района).

Первичная подготовка данных

Для начала подготовим данные по дорожной сети: обрежем их по границам района интереса, спроецируем куда нужно, отфильтруем лишнее. Распаковываем архив RU-LEN.7z в удобное расположение, лучше без кириллицы и пробелов в путях. Затем запускаем QGIS и добавляем два векторных слоя из папки .../RU-LEN/data/:

• highway-line
• boundary-polygon

Нужно проявить осторожность с кодировками, и добавлять слои с кодировкой UTF-8 (кодировку можно выбрать в меню выбора векторного слоя при его добавлении). Выглядит это как-то так (цвета могут быть произвольными):

В первую очередь, избавимся от всего лишнего. Найдём полигон Тосненского района и сохраним его в отдельный слой. Для этого в таблице атрибутов слоя boundary-polygon (доступно через контекстное меню слоя) найдём объект с соответствующим содержанием поля NAME. Для этого можно воспользоваться поиском по запросу, либо просто отсортировать колонку и выделить нужную строку. Затем через контекстное меню слоя выбираем "Сохранить как", определяя новое расположение для файла с границами Тосненского района (в моём случае Tosnensky_boundary.shp), и не забываем активировать флаг "Сохранить только выделенные объекты".

Затем необходимо обрезать набор геоданных по дорожной сети по этим границам. Для этого воспользуемся инструментом обрезки (Вектор - Инструменты геообработки - Обрезка). В качестве исходного слоя выбираем highway-line, в качестве слоя обрезки Tosnensky_bondary, результат сохраняем в удобное новое расположение (в моём случае Tosnensky_roads.shp).

Теперь, чтобы избежать возможных проблем с расчётом длин участков дорог (и некоторых других) в будущем, перепроецируем результирующий слой Tosnensky_roads в подходящую систему координат. Для расчётных задач рекомендуется использовать поперечную цилиндрическую проекцию (UTM или Гаусса-Крюгера) для вашей зоны. В данном случае используем WGS 84 / UTM zone 36N (EPSG:32636), для этого через меню "Сохранить как" создаем копию слоя Tosnensky_roads, в выпадающем списке "Система координат", выбрав UTM zone 36N и задав удобное расположение для нового файла (Tosnensku_roads_utm).

Аналогичным образом перепроецируем слой с границами Тосненского района, на всякий случай (получаем дополнительно Tosnensky_boundary_utm). Данные почти готовы! Теперь нужно отфильтровать геоданные по дорогам с учётом требований к моделированию. Это первый тематический шаг подготовки: отбор только тех объектов дорожной сети, которые могут быть использованы для перемещения предполагаемыми транспортными средствами. Здесь вы должны решить, по объектам каких типов будет разрешено перемещаться, например, для моделирования пешеходного движения можно оставить все типы дорог; для легковых автомобилей исключить пешеходные дорожки, тропы, просеки; для грузовых автомобилей исключить дороги с соответствующими ограничениями на массу, и так далее. Шаг этот очень ответственный. При этом подход к фильтрации сильно зависит от качества исходных данных и их детализации. Вспоминаем, что имеем дело с OSM, где различия между объектами существуют на уровне тегов, которые записаны в таблице атрибутов в поле HIGHWAY. Посмотрим на список уникальных значений для нашего набора данных:

path steps footway residential primary service unclassified secondary road track construction tertiary raceway tertiary_link trunk secondary_link proposed pedestrian bridleway primary_link living_street trunk_link

Все эти теги описаны на ресурсах OSM, даже с картинками. Нам ещё предстоит вернуться к этим описаниям, а сейчас важно определить, по объектам каких типов может нормально передвигаться легковой автомобиль, а по каким нет. Очевидно, что нужно исключить объекты следующих типов (пригодны только для пешеходов, проектируются, строятся и т.д.):

path, steps, footway, construction, proposed, pedestrian, bridleway

Также много объектов с типом unclassified. При осуществлении серьезной работы, конечно, стоит прояснить природу каждого из таких объектов по вспомогательным данным (и внести информацию в OSM, конечно же). В данном случае оставим их как подходящие для работы. Объекты всех отобранных типов нужно удалить. Для этого переходим в таблицу атрибутов слоя Tosnensky_roads_utm, выбираем инструмент "Выбрать по выражению" и формируем запрос вида:

"HIGHWAY" IN (  'path' ,  'steps' ,  'footway' ,  'construction' ,  'proposed' ,  'pedestrian' ,  'bridleway' )

отмечая в списке, соответственно, выбранные вами типы дорог, которые должны быть проигнорированы. Нажимаем кнопку "выбрать" и видим на карте лишние дороги.

Переходим в режим редактирования и удаляем выбранные объекты, сохраняем изменения. Если что, все исходные данные остаются в других слоях. Общая подготовка закончена.

Определение транспортных характеристик данных о дорожной сети

Второй этап подготовки является самым важным и трудоёмким в реальных условиях: необходимо каждому участку дороги назначить некоторую среднюю скорость перемещения по нему. Для примера мы используем простую логику: каждому семейству объектов (по тегу) назначим общую ожидаемую скорость, не вникая в каждый отдельный объект. При серьезном исследовании, опять же, стоит более внимательно отнестись к этому этапу, по возможности используя дополнительные материалы. Итак, рассмотрим те категории, которые у нас остались после фильтрации:

• trunk - важнейшие и крупнейшие дороги, например, для нашей территории, Московское шоссе. Ожидаемая скорость 90 км/ч
• primary - крупные шоссе, следующий уровень после trunk. Ожидаемая скорость 90 км/ч
• secondary - относительно крупные дороги, следующий уровень после primary. Ожидаемая скорость 60 км/ч
• tertiary - обычные автомобильные дороги между небольшими населенными пунктами. Ожидаемая скорость 60 км/ч
• living_street - жилые зоны, где у пешеходов явное преимущество в праве передвижения. Ожидаемая скорость 15 км/ч
• residential - автомобильные дороги в жилых кварталах. Ожидаемая скорость около 40 км/ч
• service - сервисные подъезды, въезды и проч. Ожидаемая скорость 30 км/ч
• road - автомобильная дорога неизвестного типа. Примем скорость 60 км/ч
• track - грунтовые дороги, обычно для сельхоз-техники. Примем скорость 30 км/ч
• raceway - дороги для автомобильных видов спорта. Примем скорость 90 км/ч
• tertiary_link - участки, соединяющие tertiary с другими tertiary или дорогами других типов. Примем скорость 40 км/ч
• secondary_link - участки, соединяющие secondary с другими secondary или дорогами других типов. Примем скорость 40 км/ч
• primary_link - участки, соединяющие primary с другими primary или дорогами других типов. Примем скорость 40 км/ч
• trunk_link - участки, соединяющие trunk с другими trunk или дорогами других типов. Примем скорость 40 км/ч
• unclassified - дороги без тега. Примем скорость 40 км/ч

Принятые значения достаточно условны, но позволят продемонстрировать принципы дальнейшей работы. Назначаемые скорости довольно существенно зависят от конкретной территории, местного законодательства и так далее. Важно заметить, что у некоторых участков (обычно их немного) в свойствах заполнено значение MAXSPEED, оно может быть как численным (в км/ч), так и вида RU:urban. Расшифровку этих обозначений можно найти, к примеру, здесь и использовать при определении скоростей. В целом этот этап является ключевым, и в боевых условиях здесь важно внимательно и тщательно определить характеристики дорожной сети.

Для того, чтобы применить принятые характеристики, воспользуемся калькулятором полей, запустив его из таблицы атрибутов слоя Tosnensky_roads_utm. Создадим новый атрибут SPEED как целочисленный и заполним его следующим выражением:

CASE
WHEN  "HIGHWAY" = 'trunk' THEN 90
WHEN  "HIGHWAY" = 'primary' THEN 90
WHEN  "HIGHWAY" = 'secondary' THEN 60
WHEN  "HIGHWAY" = 'tertiary' THEN 60
WHEN  "HIGHWAY" = 'living_street' THEN 15
WHEN  "HIGHWAY" = 'residential' THEN 40
WHEN  "HIGHWAY" = 'service' THEN 30
WHEN  "HIGHWAY" = 'road' THEN 60
WHEN  "HIGHWAY" = 'track' THEN 30
WHEN  "HIGHWAY" = 'raceway' THEN 90
WHEN  "HIGHWAY" = 'tertiary_link' THEN 40
WHEN  "HIGHWAY" = 'secondary_link' THEN 40
WHEN  "HIGHWAY" = 'primary_link' THEN 40
WHEN  "HIGHWAY" = 'trunk_link' THEN 40
WHEN  "HIGHWAY" = 'unclassified' THEN 40
END

Здесь мы проверяем принадлежность каждого объекта слоя к одной из категорий и назначаем скорость в соответствии с определенным ранее значениями. Очень удобно использовать условный оператор CASE.

После этой операции у каждого объекта слоя есть характеристика скорости. Теперь, при необходимости, можно поправить её для отдельных объектов в таблице атрибутов с учётом изучения реальных характеристик территории. Мы этот шаг пропустим и приступим к следующему - расчёту времени, - которое потребуется для преодоления каждого отдельного участка дорожной сети. Для этого нужно сначала рассчитать длины всех объектов слоя. В калькуляторе атрибутов создаем новое поле LENGTH как десятичное число по формуле

$length/1000

По умолчанию расчёт длин и площадей производится в единицах проекции, т.е. в нашем случае в метрах (квадратных метрах для площади). Делением на 1000 мы приводим расстояния к километрам. Следующий шаг - расчёт времени (предпочтительно в минутах), которое потребуется для преодоления каждого отдельного участка дороги с учётом его длины и ожидаемой скорости. В калькуляторе атрибутов создаем новое поле TIME как десятичное число по формуле

"LENGTH"/"SPEED" * 60

Теперь у нас есть для каждого участка дорожной сети оценка времени, требуемая для его преодоления. На этом подготовка данных в QGIS завершена, можно приступать к моделированию и переходить в GRASS. Не забудьте сохранить изменения в слое и выйти из режима редактирования.

Моделирование транспортной доступности

Подготовка проекта и данных в GRASS

Запускаем графический интерфейс GRASS GIS 7. В первую очередь, необходимо создать так называемую "локацию" или "область". Если вы запускаете GRASS впервые, укажите папку для хранения данных (database directory). На скриншоте ниже вы видите, что в качестве такой папки выбрана E:\transport_article

Теперь нужно создать новую локацию (location). Для этого используем единственную активную кнопку New и попадаем в меню, где указываем название директории для хранения данных области внутри database directory, и заголовок для неё. Для простоты в примере использованы те же названия: "transport_article". Таким образом, в папке E:\transport_article будет создана папка "transport_article", доступная для выбора по такому же заголовку.

Далее необходимо выбрать систему координат (мы уже знаем, что нам нужна UTM 36N, которую можно найти по номеру или по названию

Трансформацию датума на следующем шаге, при использовании UTM, можно выбрать произвольно, т.к. эффект будет одинаков (обратите внимание, что в опции с преобразованием все смещения нулевые).

Далее проверяем, что все параметры правильные, и завершаем создание области.

По завершению создания области будет предложено установить охват и разрешение региона. Это мы сделаем позднее, а сейчас откажемся.

На предложение же создать новый набор ("mapset") ответим утвердительно, и создадим его с именем tosno, что будет отсылкой к конкретному району исследования.

В итоге в стартовом меню появляются созданные локация (tranport_article) и набор (tosno), выбрав которые можно запускать собственно GRASS (кнопка Start GRASS session)

Наконец-то, начинается что-то интересное! Для начала нужно установить плагин, который и будет ответственным за построение изохрон: v.isochrones. Для этого обратимся к меню Settings - Addons extensions - Install extensions from addons.

В открывшемся меню в подгруппе vector находим модуль v.isochrones и устанавливаем его.

В консоли эта операция тоже выполняется просто:

g.extension extension=v.isochrones operation=add

Если у вас возникают какие-либо проблемы с доступом к репозиторию с модулями GRASS (плагины не доступны для загрузки через меню "extensions"), это не критично: можно загрузить исходный Python-скрипт по ссылке, и запускать модуль, обращаясь напрямую к скрипту через File - Launch script. При таком сценарии у вас откроется тот же графический интерфейс, который будет рассмотрен в дальнейшем.

Теперь давайте добавим в GRASS те данные, которые мы подготовили в QGIS. Для этого используем функцию импорта векторных данных (v.in.ogr) или v.import. В графическом интерфейсе искомое меню открывается через File - Import vector data - common import formats. Выбираем наш набор геоданных "Tosnensky_roads_utm.shp":

v.import input=E:\qgis_grass_transport\data\Tosnensky_roads_utm.shp layer=Tosnensky_roads_utm output=Tosnensky_roads_utm

В консоли мы видим результат импорта и базовые сведения о данных, а в окне карты видим простую визуализацию.

Теперь создадим точечный слой, в котором будут содержаться местоположения, относительно которых требуется рассчитать изохроны. Это удобнее делать в интерактивном режиме в окне карты. Перейдем из режима 2D карты в режим оцифровки:

Затем в крайнем слева выпадающем списке выберем опцию создания нового векторного слоя и зададим для него какое-нибудь подходящее название, например, start_points. При этом атрибуты нам не нужны? и соответствующий флаг в меню создания слоя можно снять.

Теперь выбираем инструмент "Создать новую точку" и добавляем интересующие нас локации. Для демонстрации поставим точки отсчёта в центры населенных пунктов г. Тосно и г. Любань. После этого выходим из режима оцифровки, соглашаясь на сохранение изменений.

Обратим внимание, что теперь на вкладке Слои (Layers) имеем два набора: импортированную дорожную сеть и созданные начальные точки.

Продолжаем подготовку данных. Поскольку в планах у нас занятия сетевым анализом, нужно набор линейных данных преобразовать в набор сетевых данных с помощью инструмента v.net, который доступен в меню Vector - Network analysis - Network preparation. Всё, что нам нужно, это применить операцию nodes к нашему слою "Tosnensky_roads_utm". На вкладке "Необходимо" выбираем метод nodes, на вкладке Arcs выбираем слой "Tosnensky_roads_utm@tosno", на вкладке Опционный задаем имя для результирующего набора. Например, "Tosnensky_roads_utm_net".

v.net -c input=Tosnensky_roads_utm operation=nodes output=Tosnensky_roads_utm_net

Следующий (и последний) шаг подготовки данных - это задание региона и разрешения для моделирования с помощью модуля g.region, который доступен в меню "Settings - Регион - Установить регион". На первой вкладке в выпадающем списке Set region to match vector maps выбираем слой "Tosnensky_roads_utm_net", созданный на предыдущем шаге, на вкладке "Границы" устанавливаем флаг "Подогнать регион под разрешение" и на вкладке "Разрешение" назначаем для "2D grid resolution" подходящее значение. Этот последний параметр очень важен - это размер ячеек сетки, до которых и будет производиться расчёт времени перемещения. Он должен быть достаточно небольшим, чтобы изохроны были плавными, но и не совсем маленьким, потому что чем подробнее сетка, тем более ресурсоёмкими получаются вычисления. Для территории подобной демонстрационной подходит размер 20 метров. Другие параметры не изменяются.

g.region vector=Tosnensky_roads_utm_net res=20 -a

Теперь всё готово для создания самих изохрон и некоторых других выходных данных.

Создание изохрон

Модуль v.isochrones позволяет работать в двух режимах: создавать непрерывные изохроны, то есть поля затрат на перемещение (с использованием модуля r.cost), и просто откладывать преодолеваемое расстояние вдоль дорог (что хорошо подходит для случая, когда перемещение за пределами транспортной сети невозможно совсем - с использованиям модуля v.net.iso). Нас интересуют непрерывные изохроны. Запускаем модуль v.isochrones, набрав такую команду в консоли:

v.isochrones

Можно также запустить загруженный скрипт через File - Launch script. На первой вкладке в качестве дорожной сети выбираем сетевой набор данных "Tosnensky_roads_utm_net@tosno", название атрибута с ценой - SPEED, стартовые точки - "start_points@tosno", выбираем удобное название для слоя с изохронами (например, "isochrones_rcost"), выбираем метод r.cost и (важное!) собственно значения изохрон в минутах через запятую. Здесь значения могут быть любыми в зависимости от ваших задач. Зададим следующие отсечки по времени: 15, 30, 60, 90, 120 и 150 минут. Далее, на вкладке r.cost задаются очень важные параметры: название для поверхности времени перемещений (опциональный продукт, но очень интересный, сохраним результат в слой timemap_rcost и посмотрим, что в нем будет), скорость для преодоления пространства вне дорожной сети (вышли на обочине и пошли пешком через поле, установим 5 км/ч) и доступная для расчётов память. Обычно хватает 1000 МБ.

v.isochrones map=Tosnensky_roads_utm_net@tosno roads_layer=1 cost_column=SPEED start_points=start_points@tosno isochrones=isochrones_rcost time_steps=15,30,60,90,120,150 timemap=timemap_rcost memory=1000 method=r.cost

После выполнения обработки в дереве слоёв появляются два новых: изохроны и поверхность затрат времени. Это и есть искомый результат. Оформим его несколько позднее, пока просто взглянем:

Выглядит неплохо! Изучим данные немного детальнее в QGIS, для этого экспортируем данные в стандартные ГИС-форматы с помощью команд v.out.ogr и r.out.gdal, они доступы в меню File - экспорт векторного слоя и File - экспорт растрового слоя соответственно.

v.out.ogr input=isochrones_rcost@tosno output=E:\qgis_grass_transport\data\isochrones.shp format=ESRI_Shapefile

r.out.gdal input=timemap_rcost@tosno output=E:\qgis_grass_transport\data\timemap.tif format=GTiff

Теперь можно использовать данные в других ГИС и делиться ими!

Представление результатов в QGIS

Посмотрим на данные в QGIS, и оформим простую карту. Сначала добавим растровый слой timemap.tif, что он из себя представляет? Оказывается, это очень интересный набор данных, как бы предшествующий изохронам - растр с тем самым разрешением 20х20 метров, указанным на этапе задания области, в каждой ячейке которого хранится время в минутах, сколько туда добираться из указанных точек. Такой набор данных может быть очень полезен при разнообразном растровом анализе как один из входных слоёв.

Слой же с изохронами содержит полигоны, охватывающие территории, соответствующие диапазонам времени перемещения: от 0 до 15, от 15 до 30, и так далее. Диапазоны мы ранее задавали самостоятельно. Важно, что в атрибутах хранятся данные о временном диапазоне, а внешние полигоны (более дальние) не содержат внутренних (ближних), такие данные открывают широкие возможности для разнообразного анализа.

Применив простую заливку по уникальным значениям поля traveltime и подложив какую-нибудь базовую карту, можно быстро получить достаточно симпатичную визуализацию изохрон.

Хотя ценнее в данном случае, конечно, сами данные.

Итак, рассмотрен простой подход к транспортному моделированию способом построения изохрон в GRASS GIS. Возможности такого подхода весьма ограничены, особенно в аспектах, касающихся учёта разнородных препятствий, и так далее. Тем не менее, при внимательной подготовке исходных данных о дорожной сети, можно получать интересные результаты, которые применимы при различных сценариях анализа транспортной доступности.