Матрунич Консалтинг

Анализ данных, визуализация, маркетинговые исследования, язык R

Tag: OpenStreetMap

Как определить длину пути между адресами

Проблема и цель проекта

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

Подбор инструментов

Преимущества свободных данных и софта

Для реализации проекта предлагается использовать данные и программное обеспечение, которые распространяются под свободными лицензиями. Преимущества такого подхода перед использованием картографических сервисов от Google или Яндекса следующие:

  • Проект подразумевает обработку данных для внутренних целей, а лицензионные договоры Google и Яндекса требуют использования их сервисов только для создания общедоступных публикаций в виде карт. Т.е. использование услуг этих компаний в нужных нам целях будет незаконным. Конечно, вероятность того, что Яндекс обратит внимание на незначительную нагрузку, которую мы генерируем, очень мала. Но с точки зрения развития бизнеса или его продажи в будущем, лучше сразу отказаться от включения в работу предприятия заведомо нелегальных алгоритмов.
  • В процессе работы будут постепенно «всплывать» ошибки в используемой базе данных: например, какие-то адреса будут неправильно преобразовываться в координаты. Если мы работаем на основе картографического сервиса Яндекса, то в случае обнаружения ошибок (в нашем опыте такие ошибки у Яндекса были), каждый ошибочный адрес придётся обрабатывать как исключение, поскольку сам Яндекс по нашему запросу исправит эти ошибки только через несколько месяцев. При использовании данных из базы OpenStreetMap, мы можем сами внести исправление в базу, и оно практически сразу станет доступно для нашего приложения.
  • При необходимости внести какие-либо изменения в работу алгоритмов по геокодированию или прикладке маршрутов, мы можем отказаться от использования открытого API, а запустить эти сервисы на своём оборудовании, поскольку и база данных дорожной сети, и программное обеспечение доступны на условиях свободных лицензий.

Схема работы

Алгоритм вычисления расстояния между двумя адресами состоит из двух этапов:

  1. Геокодирование, т.е. преобразование текстовой строки, содержащей адрес, в географические координаты.
  2. Прокладка маршрута и вычисление его длины.
Open Source Routing Machine

Open Source Routing Machine

Геокодирование российских адресов

Международный сайт OpenStreetMap.org для геокодирования использует программу Nominatim. У нас есть опыт работы с API Nominatim через статистический пакет R. Исходные коды R-скриптов доступны в репозитории на github.

Для национального проекта OpenStreetMap.ru разработано собственное приложение для поиска адресов и точек интереса. Приложение «заточено» под специфику адресов в России, »понимает« нумерацию домов в духе »14 К2 С1« и может искать дома с двойным адресом. У поисковика доступно API, документация опубликована в форуме. Например, для определения координат дома по адресу «Псков, Маркса, 23» следует отправить запрос вида http://openstreetmap.ru/api/search?q=Псков, Маркса, 23. Ответ вернётся в формате JSON, координаты хранятся в полях lat и lon.

Для снижения количества ошибок геокодирования и уменьшения объёма лишней информации в ответе в строку запроса рекомендуется добавить параметры точки начала поиска (координаты центра города; аргументы запроса lat и lon), поиск только по адресам (stype=addr), количество ответов установить равным 1 (cnt=1).

Прокладка маршрута

Корифей OpenStreetMap и ведущий новостного блога ШТОСМ Илья Зверев порекомендовал использовать для маршрутизации Open Source Routing Machine. У разработчиков запущен демо-сервер с доступным API (Илья просил не более 20 запросов в минуту, также существует ряд требований).

Предположим, мы хотим узнать длину маршрута между адресами Маркса, 23 и Рижский, 69 в городе Пскове. Через API OpenStreetMap.ru мы узнаём координаты:

Маркса, 23
57.818023681640625, 28.342592239379883
Рижский, 69
57.812679290771484, 28.28318214416504

Полученные координаты мы отправляем на API Open Source Routing Machine в виде запроса http://router.project-osrm.org/viaroute?loc=57.818023681640625,28.342592239379883&loc=57.812679290771484,28.28318214416504. Ответ приходит в формате JSON, в поле total_time содержится предполагаемое время на дорогу (в нашем случае 467 секунд), в поле total_distance — длина маршрута (4582 метра).

OSRM позволяет прокладывать маршрут сразу между последовательностью точек, но не более 25.

Подготовлено для кафе «Моя Италия»

Рисуем картограмму в R средствами ggplot2

Что такое картограмма

Плотность населения районов Псковской области. Картограмма сделана в ggplot2

Плотность населения районов Псковской области. Картограмма сделана в ggplot2


Картограмма – это график, демонстрирующий распределение какого-либо статистического показателя по географическим единицам. На картограмме границы географических единиц (континенты, страны, районы, кварталы) обозначены линиями, а площадь этих единиц залита разными цветами. Цвет показывает категорию или интенсивность показателя. Мы построим картограмму плотности населения районов Псковской области. Чем выше плотность населения района, тем более интесивным цветом он закрашен.

Общая схема создания картограммы с помощью пакета ggplot2 в среде R

Для построения картограммы необходимо осуществить такие шаги:

  1. Получить массив с географическими (пространственными) данными о границах районов.
  2. Получить массив с теми статистическими данными, распределение которых мы хотим отобразить на картограмме (плотность населения в нашем случае).
  3. Связать пространственные и статистические данные, чтобы программа знала, к какому куску плоскости относится данная статистика.
  4. Построить картограмму.

Получаем данные об административных границах районов Псковской области

Что такое [векторные] пространственные данные об административных границах? Это массив данных, в котором хранится информация о координатах точек, составляющих линию границы каждого района. Самый законно-бесплатный способ получить такой массив – это взять его из базы данных OpenStreetMap. В ней хранятся данные с административными границами нашего дорого Отечества (и всех других стран), федеральных округов, субъектов Федерации. По Псковской области также есть информация о границах районов и о границах некоторых волостей.

Спасибо активистам сообщества gis-lab.info, которые на регулярной основе “выдергивают” из всего громадного массива данных OpenStreetMap актуальные куски для стран, ранее входивших в СССР, и преобразуют их в разные форматы. Нам пригодятся данные по административным границам муниципальных районов России, доступные на сайте gis-lab.info в формате ESRI Shapefile. К сожалению, самая последняя версия этого файла не содержит информацию о границе Струго-Красненского района. Видимо, надо что-то исправлять в OpenStreetMap. Поэтому я в архиве нашёл версию данных от 1 января 2012 года, в которых граница ещё присутствует 🙂 Воспользуемся этим файлом.

# Скачиваем файл
download.file('http://gis-lab.info/data/yav/adm/20120101/adm6_district.7z', 'adm6_district.7z')

# Разархивируем. Системная команда для Linux.
# В Windows, наверно, следует прописать путь к архиватору 7z.
system('7z x adm6_district.7z')

# Читаем shape-файл.
library(rgdal)
# Первый параметр у readOGR - имя файла. Второй - название слоя.
shapes <- readOGR('adm6_district.shp', 'adm6_district')

Полученный массив содержит информацию по границам муниципальных районов всей России. Оставим только Псковскую область.

shapes <- shapes[shapes$ADM4_NAME=='Псковская область', ]

Как выяснилось при создании картограммы в ggplot2, функция geom_map не рисует дырки в полигонах, по крайней мере, я не нашёл способ это сделать. Поэтому выходит, что Псков и Великие Луки лежат на территории Псковского и Великолукского районов. И ещё один нюанс, что по умолчанию полигоны заливаются в той последовательности, в которой идут их идентификаторы. Скорее всего, есть более элегантный способ решить этот вопрос, а я пока иду "грязным" путём, заменив идентификаторы указанных городов в shape-массиве:

# При экспорте из OSM объекты получают идентификаторы со знаком "минус".
# Поэтому объект с любым положительным идентификатором будет отрисовываться после всех остальных.
shapes$OSM_ID[shapes$NAME=="Великие Луки"] <- 2
shapes$OSM_ID[shapes$NAME=="Псков"] <- 1

Функции из пакета ggplot2 любят (требуют) работать с данными в формате data.frame. Поэтому преобразуем массив класса 'sp' в массив класса data.frame:

# Преобразуем sp-объект в обычный data.frame, подходящий для ggplot2.
# Параметр region указывает, по какому полю группировать точки.
# В нашем случае каждый район обладает уникальным идентификатором OSM_ID из OpenStreetMap
shapesm <- fortify(shapes, region='OSM_ID')

Если вы захотите каждый район отметить своим названием, то функция geom_text из ggplot2 захочет, чтобы вы указали ей координаты тех мест, где ей нужно будет разместить текстовые строки с названиями. В общем случае, нам сгодятся и координаты центров районов, т.е. центроидов.

# получаем координаты центроида каждого района,
# чтобы разместить названия районов
# пользуемся функцией centroid.polygon из пакета maps
# вызываем её явно из пакета, т.е. с использованием :::, т.к.
# она внутренняя и не доступна в обычном пространстве имён
lp <- ddply(shapesm, .(id), maps:::centroid.polygon)
names(lp) <- c('OSM_ID', 'lpx', 'lpy')

Далее нам придётся бороться с ограничениями этого метода: метки Пскова и Лук накладываются на метки Псковского и Великолукского районов.

# Раздвигаем метки по вертикали.
# Координаты для сдвига получены опытным путём.
lp$lpy[lp$OSM_ID==-956303] <- 58
lp$lpy[lp$OSM_ID==1] <- 57.7
lp$lpy[lp$OSM_ID==-957758] <- 56.52
lp$lpy[lp$OSM_ID==2] <- 56.25

Итак, исходные пространственные данные для построения картограммы готовы. Осталось получить фактическую информацию о плотности населения.

Получаем данные о плотности населения районов Псковской области

Данные о населении районов и об их площади я взял из Википедии и сохранил в csv-файл.

# Читаем массив с данными по площади и населению районов
data <- read.table('/home/sas/1irr/data/regions_square_pop.csv', header=T, sep=';', dec='.')

# Считаем плотность
data$density <- data$population / data$square

# объединяем данные по плотности с координатами центроидов
data <- join(data, lp)

Дело остаётся за малым: построить график.

Строим график в ggplot2

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

# Функция по "прилизываю" графика.
polish_map <- function(x) {
  # Добавляем названия районов
  x <- x + geom_text(family='Courier', size=4, fontface='bold')
  # Удаляем линии сетки (промежуточные и основные)
  x <- x + theme(panel.grid.minor=element_blank(), panel.grid.major=element_blank())
  # Удаляем подписи значений на осях
  x <- x + scale_y_continuous(breaks=NULL) + scale_x_continuous(breaks=NULL)
  # Удаляем названия осей
  x <- x + ylab('') + xlab('')
  # Передвигаем легенду в верхний левый угол.
  # У Псковской области там как раз пустое пространство.
  x <- x + theme(legend.position=c(1, 1), legend.justification=c(1,1)
  # Можно залить фон графика белым: добавьте след.строку в параметры theme выше.
 # panel.background=element_rect(fill='white', colour='white')
               )
  return(x)
}

Рисуем график.

# Создаём переменную класса ggplot
# В ней используется массив data, 
# параметры x и y - это имена колонок из data, которые используются для определения
# горизонтальной и вертикальной позиции,
# map_id - по какой колонке происходит сопоставление данных и геообъектов.
# label - откуда брать текст.
m <- ggplot(data, aes(x=lpx, y=lpy, map_id = OSM_ID, label=shortname))

# Указываем, что в нашем графике данные будут визуализироваться в виде карты.
# При этом цвет заливки регионов карты определяется содержанием колонки density.
# Параметр colour указывает на то, что обрамление регионов должно присутствовать, но не 
# сопоставлено какой-либо переменной.
# Параметр map указывает, где хранятся данные о карте.
m <- m + geom_map(aes(fill=density), colour=1, map=shapesm) 

# Функция expand_limits() задаёт диапазон данных по осям. 
# Без этой функции края области будут обрезаны.
m <- m + expand_limits(x = shapesm$long, y = shapesm$lat)

# В параметре name задаётся название легенды. В нашем случае
# мы дополнительно используем функцию expression, чтобы добавить
# верхний индекс.
# В параметре trans задаём способ трансформации исходных данных - 
# поскольку у нас два города с высокой плотностью, а остальные - с 
# очень низкой, то применяем логарифмичекое преобразование, чтобы 
# сгладить разрыв.
# В labels указываем функцию round(), которая укругляет значения,
# отображаемые в легенде.
m <- m + scale_fill_gradient2(name = expression(paste("Человек/км"^2)), trans='log', labels=round)
m <- polish_map(m)
m <- m + ggtitle('Плотность населения районов\nПсковской области')
ggsave('density.png')

Powered by WordPress & Theme by Anders Norén