Считаем плотность дорог с QGIS на основе данных OpenStreetMap

Как-то раз мне понадобилась карта плотности дорожного покрытия одного из регионов России. Конечно же, как и любой нормальный человек, я сразу же открыл Google. К своему удивлению, я не обнаружил ничего касательного своего региона, да и вообще России. Сложность добавило одно из моих требований — карта должна была быть под свободной лицензией. Я достаточно быстро понял, что продолжать поиски бессмысленно. Но что же делать? Делать карту самому!
Disclaimer: Автор данной статьи человек достаточно далёкий от геодезии и ГИС. А способ, описанный в данной статье, не претендует на абсолютную точность и правильность. Автор допускает существование более простого, правильного и лаконичного варианта решения данной задачи.
Так передо мной оформился следующий порядок действий:
- Скачать данные с OSM
- Отфильтровать данные, оставив только дороги федерального и регионального значения, а так же полигон региона
- Импортировать дороги и полигон региона в одну из ГИС
- ???
- PROFIT!!!
Скачиваем данные
Нужная мне область конечно же оказалась “слишком велика для экспорта в качестве данных” для стандартного средства экспорта OSM. Поэтому у меня осталось два варианта — воспользоваться Overpass API или скачать уже готовую выгрузку данных. Первый вариант позволяет сразу скачать только нужные данные, что очень хорошо. Но мои запросы на требуемой области уходили в таймаут, что совершенно не удивительно для такого размера. Ну что же, не будем его мучить. Придётся качать готовую выгрузку. Мой выбор пал на Geofabrik, где я легко отыскал требуемую область. Скачивал в формате *.pbf. После чего я сразу распаковал его в более привычный и читабельный xml, т.е. *.osm.
Распаковка. pbf to osm
Эту процедуру я произвёл простенькой утилитой Osmconvert. Команда для распаковки в моём случае более чем тривиальна: osmconvert north-caucasus-fed-district-latest.osm.pbf -o=north-caucasus-fed-district-latest.osm На распаковку ушло некоторое время, а готовый файл оказался внушительного(по сравнению с исходным) размера:

Фильтрация
Меня интересовали автомагистрали, важные дороги и дороги регионального и областного значения. Для фильтрации данных применил Osmfilter: osmfilter north-caucasus-fed-district-latest.osm — keep=”highway=motorway =trunk =primary =secondary” -o=north-caucasus-fed-district-filtered.osm

Полигон региона

Существует онлайн конвертер osm отношения в полигон с единой геометрией — polygons.openstreetmap.fr Я выбрал GeoJSON в качестве выходного формата данных. Так же там можно сразу увидеть наглядный результат:

Выбор ГИС
Я перебрал несколько вариантов, в том числе и знаменитый ArcGIS. В конце концов я остановился на QGIS. Почему? Функциональность, относительно простой интерфейс, наличие большого сообщества пользователей как на русском, так и на английском языке, лицензия GNU GPL.
Импорт
Я использовал последнюю версию QGIS — 2.18, где уже сразу встроена поддержка OSM. В версиях < 2 QGIS OSM Plugin для работы с OSM необходимо устанавливать отдельно.
Перед началом импорта нужно сделать небольшое отступление и сказать пару слов о выборе системы координат. По умолчанию в качестве системы координат для нового проекта выбирается EPSG:4326 — WGS 84, которая является одной из самых известных и распространённых. Всё бы хорошо, только WGS 84 географическая система координат, где единицами измерения являются градусы. В моём же случае необходимы линейные единицы измерения.
Я выбрал систему координат UTM со своей зоной — WGS 84 / UTM zone 38N. Для более больших областей, не влезающих в зоны, необходимо искать другую подходящую систему координат. Сами зоны можно посмотреть тут: http://gis-lab.info/qa/proj-sk-faq.html#19. Там же можно узнать и другую полезную информацию о системах координат.

Так же включил преобразование координат ‘на лету’ в свойствах проекта:

Теперь непосредственно сам импорт. Идём в меню: Вектор > OpenStreetMap > Импортировать топологию… и выбираем свой файл с отфильтрованными данными.

После этого: Вектор > OpenStreetMap > Экспортировать топологию OSM SpatiaLite… и выбираем получившийся выходной файл SpatiaLite. В тип объектов: Полилинии. Загружаем теги из БД и выбираем нужные теги. Я выбрал только два. Также необходимо, чтобы галочка “Добавить результат в проект” была активирована.

После чего дороги должны появится:

Теперь загружаем полигон региона. Я воспользовался онлайн конвертером mapshaper.org GeoJSON в Shapefile. Последний формат прекрасно импортировался в QGIS.

И выбираем полученный файл с расширением *.shp:

На этом импорт данных закончен. Вот наши два слоя, с которыми далее предстоит работать:

Для наглядности можно отключить заливку слоя региона. Делается в свойствах слоя:

Так уже лучше:

Тут стоит снова вспомнить о системах координат, чтобы не получить ошибки на следующих этапах. Дело в том, что наши два слоя импортировались в WGS 84. А работать я собираюсь в WGS 84 UTM. Поэтому необходимо сохранить наши слои в последней. Нажимаем правой кнопкой по слою и выбираем Сохранить как… В открывшимся окне необходимо выбрать формат, имя файла, систему координат и активировать галочку изменения размера. Также необходимо, чтобы “Добавить результат в проект” была активна.

Тоже самое проделываем и со вторым слоем.
Можно заметить, что некоторые дороги выходят за контур. В моём случае эти хвосты были ненужны, поэтому я решил их обрезать. Делается это инструментом Пересечение. Найти его можно в меню: Вектор > Геообработка. В качестве исходного тот слой, который нужно обрезать, т.е. наши дороги. В качестве слоя пересечения — полигон региона. Выполнение данной операции занимает некоторое время.

Получится ещё один слой — Пересечение. Старый слой с дорогами можно удалять.

Теперь необходимо проверить единицы измерения. Должно быть так:

Теперь можно разбивать полигон региона на квадраты, внутри которых после подсчитаем длину дорог. Для этого нужен инструмент Векторная сетка. Находится тут: Вектор > Выборка. Границу сетки нужно задать по слою региона — North Caucasus. Отступ и шаг измеряется в метрах. Я выбрал сетку 20 на 20 километров. Вид сетки — полигоны.

Теперь у нас есть сетка. Изначально она создаётся с заливкой, которую я уже отключил.

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

Теперь необходимо подсчитать длину дорог в каждой из ячеек. Инструмент — Сумма расстояний. Вектор > Анализ. Линии — наш обрезанный слой с дорогами. Полигоны — обрезанная сетка. Lines length field name — имя столбца с длиной линий. Lines count field name — имя столбца с количеством линий. Можно назвать как угодно, главное, чтобы понятно было.

Теперь у нас есть новый слой Длина линии. Старую сетку можно удалять.

Можно заглянуть в таблицу нашего нового слоя(правой кнопкой по слою > Открыть таблицу атрибутов) и увидеть что получилось:

Теперь нужна наглядная визуализация полученных данных. Для этого идём в свойства слоя и выбираем стиль Градуированный знак вместо Обычного.

После необходимо выбрать столбец с длиной линии, понравившийся градиент, количество градаций и нажать “Классифицировать”. С настройкой градиента и его режима можно долго разбираться, но это отдельная тема. Легенда так же в метрах.

Получается так:

Поставленная задача решена, карту можно уже экспортировать. (Да, знаю, что тут нет легенды и масштаба)
Экспорт
Тут множество вариантов. Это и экспорт в обычную картинку, и в различные другие более специфичные форматы.
Самое простое картинка. Проекты >Сохранить как изображение. На выходе получается png:

Кроме того я хотел импортировать результат в Ai. Для этого необходимо сохранить слой в DXF. С такими настройками, чтобы сохранились цвета.

А DXF формат прекрасно открывается Ai:

Что дальше?
Добавить легенду, масштаб и другие украшательства. Но это уже в другой раз.
