R: Формирование матрицы расстояний

Непосредственное измерение расстояний между населенными пунктами для заполнения матрицы расстояний — трудоёмкая задача. В случае развитой дорожной сети:

• сложно определить кратчайший маршрут движения между пунктами;
• велик объём работы (для \(n\) населенных пунктов число измерений составит \(n^2 - n \over 2 \));
• сложно объединить матрицы расстояний нескольких районов.

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

Сформировать матрицу кратчайших расстояний между всеми парами вершин произвольного связного графа можно при помощи функции floyd.warshall.all.pairs.sp из пакета RBGL. В качестве единственного аргумента эта функция принимает объект класса graph, поддерживаемый одноименным пакетом. Оба этих пакета являются частью проекта Bioconductor, с особенностями их установки и использования можно ознакомиться на соответствующих страницах: RBGL и graph.

Виртуальный класс graph объединяет несколько подклассов, отличающихся способом описания графа. Класс graphAM хранит граф в виде матрицы связности его вершин, класс graphBAM — в виде таблицы с описанием его дуг.

Применим перечисленные инструменты на примере представления в виде графа дорожной сети Можгинского района Удмуртской республики — сформируем и сохраним в формате csv таблицу, первый и второй столбец которой соответствуют узлам дорожной сети (населенным пунктам и перекресткам), а третий столбец — протяженности соединяющих их участков дороги с твёрдым покрытием:

"from";"to";"weight"
"Кр. Яр";"Ломеслуд";3,5
"Камышлы";"Ломеслуд";1,75
"_1";"Ломеслуд";4,55
"_1";"_2";1,75
"_3";"_2";1,05
"_3";"Полянское";2,1
"_3";"Бол. Уча";1,4
"Пазял-Зюмья";"Бол. Уча";2,1
…
"Водзя";"Чежесть-Какси";2,8

На карте результат процесса выглядит так:

Всего — 111 дуг, 96 вершин (в т. ч. 21 перекрёсток и 4 населённых пункта не относящихся к Можгинскому району), примерно 3 человеко-часа работы.

Код R, формирующий матрицу расстояний между всеми населёнными пунктами Можгинского района представлен ниже:

> require("graph")
> require("RBGL")
> mozhga <- read.csv("mozhga.csv", sep=";", dec=",")
> mozhga.BAM <- graphBAM(mozhga) # Формируем граф экспериментального класса graphBAM,
> mozhga.AM <- as(mozhga.BAM, "graphAM") # преобразуем его в стабильный класс graphAM,
> mozhga.dist <- floyd.warshall.all.pairs.sp(mozhga.AM) # вычисляем матрицу расстояний
> d <- mozhga.dist[grep("_", colnames(mozhga.dist), invert=T), # и избавляемся от лишних узлов
+                  grep("_", rownames(mozhga.dist), invert=T)] # (предварены символом "_").

Объект d — искомая матрица расстояний между населенными пунктами Можгинского района.

UPD:

Вместо пакетов graph и RBGL удобнее воспользоваться пакетом igraph. В этом случае код R, формирующий матрицу расстояний, будет выглядеть следующим образом:

> require("igraph")
> mozhga.graph <- graph.data.frame(mozhga, directed=FALSE) # Формируем граф,
> mozhga.dist <- shortest.paths(mozhga.graph) # вычисляем матрицу расстояний
> d <- mozhga.dist[grep("_", colnames(mozhga.dist), invert=T), # и избавляемся от лишних узлов
+                  grep("_", rownames(mozhga.dist), invert=T)] # (предварены символом "_").

Картограммы в R: используем OpenStreetMap в качестве фона

При помощи библиотеки OpenStreetMap, R-project может использовать растровые карты проекта OpenStreetMap.

Среди ее зависимостей находятся библиотеки rJava и rgdal. Про установку библиотеки rgdal я писал ранее. Установка библиотеки rJava также требует дополнительных действий — в ее зависимости входит Java JRE и Java JDK. Ubuntu 12.04 precise использует в качестве Java-окружения OpenJDK 6, так что для удовлетворения зависимостей rJava достаточно дополнительно установить пакет openjdk-6-jdk и познакомить с ним R:

$ sudo aptitude install openjdk-6-jdk
$ sudo R CMD javareconf

После такой подготовки системы, библиотека OpenStreetMap успешно устанавливается и мы можем переходить к её использованию. С тем, как получены объекты brks, part и udm.TIK можно познакомиться в предыдущих статьях.

> map <- openmap(c(59, 49), # Широта и долгота северо-западного угла карты.
+                           c(55, 56), # Широта и долгота юго-восточного угла карты.
+                           zoom=7,    # Масштаб карты, можно посмотреть в ее URL:
+ # http://www.openstreetmap.org/#map=7/57.035/53.553 (выделен красным).
+                           type="osm") # Тип скачиваемой карты (tile server). Есть и другие.
> colors <- rev(heat.colors(length(brks), alpha=1/3))
> udm.TIK.osm <- spTransform(udm.TIK, osm()) # Меняем проекцию объекта udm.TIK на используемую
# проектом OpenStreetMap.
> str(map) # Смотрим размеры карты в пикселах (выделены красным):
List of 2
 $ tiles:List of 1
  ..$ :List of 5
  .. ..$ colorData : chr [1:422878] "#F1EEE8" "#F1EEE8" "#F1EEE8" "#F1EEE8" ...
  .. ..$ bbox      :List of 2
  .. .. ..$ p1: num [1:2] 5454655 8180387
  .. .. ..$ p2: num [1:2] 6233891 7361866
  .. ..$ projection:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. .. ..@ projargs: chr "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0
+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
  .. ..$ xres      : int 667
  .. ..$ yres      : int 634
  .. ..- attr(*, "class")= chr "osmtile"
 $ bbox :List of 2
  ..$ p1: num [1:2] 5454655 8180387
  ..$ p2: num [1:2] 6233891 7361866
 - attr(*, "zoom")= int 7
 - attr(*, "class")= chr "OpenStreetMap"
> png("ER-osm.png", width=634, height=667, units="px") # Рисуем. Размер может быть любым,
#                            карта будет отмасштабирована с сохранением отношения сторон.
> plot(map)
> plot(udm.TIK.osm, col=colors[findInterval(part, brks, all.inside=T)],
+ border="transparent", add=TRUE)
> legend("left", inset=0.1, legend=leglabs(brks, under="<", over=">"), fill=colors, bty="n")
> title(main="Доля голосов за ЕР\nв муниципальных районах\nи городских округах Удмуртии")
> dev.off()

Буду благодарен за указание способа узнать разрешение растровой карты автоматически.

Картограммы в R: меняем проекцию

Для смены проекции карты в R можно использовать библиотеку rgdal:

> install.packages("rgdal")
> library("rgdal")

В Ubuntu для работы этой библиотеки требуются пакеты libgdal1-dev и libproj-dev с их зависимостями.

Мы можем скачать и прочитать шейп-файл, как указано ранее (это, кстати, существенно быстрее), однако продемонстрирую вариант подключения к открытой базе геоданных PostGIS при помощи того же rgdal:

> # Читаем список доступных слоёв базы GEN…
> ogrListLayers(dsn="PG:host=gis-lab.info user=guest dbname=gen password=guest")
 [1] "modis_tiles"              "oblasts"                  "ru_adm2_country"
 …
[26] "wrs2_landsat_rus"
> # … и загружаем слой границ России в объект rus.
> rus <- readOGR(dsn="PG:host=gis-lab.info user=guest dbname=gen password=guest",
+ layer="ru_adm2_country", verbose=T)

Проект GIS-Lab для отображения Сибирской части России и для отображения России целиком пользуется равновеликой конической проекцией Альберса с главными параллелями 52° и 64° с.ш. и главным меридианом 105° в.д.. Проекция объекта rus —

> proj4string(rus)
[1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

— что делает его изображение разорванным по 180-му меридиану, в районе Чукотки. Сменим проекцию и нарисуем карту:

> rus2 <- spTransform(rus, CRS("+proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=105
+ +x_0=18500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs"))
> plot(rus2, col="#31A35488", border=F)

Чтобы наложить на полученную карту сетку координат, необходимо эту сетку специально подготовить:

> grd <- gridlines(rus,
+ norths = seq(10, 90, by=5), # рисуем линии каждые 5° между 10 и 90° с.ш.
+ # Отобразятся только те из них, которые попадут в область рисуемой карты.
+ easts = seq(-179.9999, 180, by=5), # Если указать -180, этот меридиан не отобразится.
+ ndiscr=360) # Число опорных точек — если указать мало,
> # координатная сетка после перепроецирования может получиться ломаной.
> # Меняем проекцию координатной сетки так же, как у карты…
> grd2 <- spTransform(grd, CRS("+proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=105
+ +x_0=18500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs"))
> # … и добавляем её к рисунку.
> plot(grd2, add=T, col="#08519C88", lty=2)

Результат приведен в начале заметки.

Перед сменой проекции мы можем упростить контур объекта rus,

> rus.thiny <- thinnedSpatialPoly(rus,
+ tolerance=0.5, # Интенсивность упрощения, чем больше, тем проще контур.
+ minarea=0.1) # Минимальная площадь сохранямых островов.

объем данных уменьшается очень существенно, но машинное время, затраченное на этот процесс, превышает время на смену проекции и отрисовку исходного объекта вместе взятое. Результат с упрощенным контуром — ниже.

Рисуем точечные картограммы в R: расставляем точки

Окончание. Предыдущие части:
Рисуем картограммы в R,
Рисуем фоновые картограммы в R: красим районы,
Рисуем картодиаграммы в R: бросаем пироги.

В R мы можем иллюстрировать распределение исследуемых объектов по территории расставляя на карте пропорциональное их численности множество точек. В результате у нас получится так называемая dot-density map. Вот такая:

Получена она при помощи следующего кода:

> q <- 1000
> s <- .2
> png("dot-density.png", width=5, height=5, units="in", res=100, bg="transparent")
> par(mai=c(.1,.1,1,.1))
> plot(udm[udm$ADMIN_LVL == 2 & !is.na(udm$ADMIN_LVL),], col="white", border=gray(.25))
> points(dotsInPolys(udm.TIK, round(tik$ER/q, digits=0)), cex=s, col="#0000ff88")
> points(dotsInPolys(udm.TIK, round(tik$KPRF/q, , digits=0)), cex=s, col="#ff000088")
> title(main="Распределение голосов за ЕР и КПРФ\nв муниципальных районах
+ и городских округах Удмуртии")
> legend("right", legend=c("ЕР", "КПРФ"), cex=.7, pt.cex=s, pch=1,
+ col=c("#0000ff88", "#ff000088"), bty="n",
+ title="Одна точка —\nтысяча голосов за…")
> dev.off()

Координаты точек для карты их плотности формирует функция dotsInPolys из библиотеки maptools.

Вы можете ознакомиться и с другими примерами (англ.) карт в R.

Рисуем картодиаграммы в R: бросаем пироги

Продолжение. Предыдущие части:
Рисуем картограммы в R,
Рисуем фоновые картограммы в R: красим районы.

У карт с заливкой районов цветом в зависимости от доли партии на выборах есть существенный недостаток: на них никак не отражается плотность населения. Из-за этого искажается представление о значении того или иного региона для результатов голосования.

Избежать этого искажения можно изображая результаты выборов при помощи круговых диаграмм («пирогов») площадь которых пропорциональна числу проголосовавших:

Создано это изображение следующим образом:

> # Отсортируем районы в порядке убывания числа действительных бюлетеней
> # (большие «пироги» отрисуются раньше, это улучшает внешний вид):
> udm.TIK <- udm.TIK[order(tik$v10, decreasing=T),]
> tik <- tik[order(tik$v10, decreasing=T),]
> # Запрашиваем структуру первой строки объекта udm.TIK…
> str(udm.TIK[1,]) # … и изучаем её:
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 1 obs. of  3 variables:
  .. ..$ OSM_ID   : int -954515
  .. ..$ NAME     : Factor w/ 341 levels "Агрикольское",..: 107
  .. ..$ ADMIN_LVL: Factor w/ 7 levels "10","2","3","4",..: 5
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 53.2 56.9
  .. .. .. .. .. .. ..@ area   : num 0.0469
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:272, 1:2] 53 53 53 53.1 53.1 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 53.2 56.9
  .. .. .. ..@ ID       : chr "205"
  .. .. .. ..@ area     : num 0.0469
  ..@ plotOrder  : int 1
  ..@ bbox       : num [1:2, 1:2] 53 56.7 53.4 57
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

Обратите внимание на содержимое слота labpt, входящего в слот polygons (выделено красным) — это координаты, к которым прикрепляется на карте название региона. Для нас эта точка будет центром будущей круговой диаграммы.

> # Соберём долготы…
> x <- sapply(slot(udm.TIK, "polygons"), slot, "labpt")[1,]
> # … и широты для размещения круговых диаграмм
> y <- sapply(slot(udm.TIK, "polygons"), slot, "labpt")[2,]
> # Сформируем матрицу с долями партий — функция draw.pie требует именно матрицу:
> part <- matrix(
+ c(tik$ER/tik$v10, tik$KPRF/tik$v10, tik$LDPR/tik$v10, tik$SR/tik$v10,
+ 1-(tik$ER+tik$KPRF+tik$LDPR+tik$SR)/tik$v10),
+ ncol=5, dimnames=list(tik$x ,c("ЕР", "КПРФ", "ЛДПР", "СР", "др.")))
> # Закажем для партий цвета:
> color <- c("#0000ff88", "#ff000088", "#ffff0088", "#00ff0088", "transparent")
> #  "#ff000088" — красный полупрозрачный (88 — уровень прозрачности)
> # Закажем библиотеку, умеющую рисовать круговые диаграммы на карте:
> library("mapplots")
> # Подготовим фон:
> png("result.png", width=5, height=5, units="in", res=100, bg="transparent")
> par(mai=c(.1,.1,1,.1))
> plot(udm[udm$ADMIN_LVL == 2 & !is.na(udm$ADMIN_LVL),], col="white", border="gray")
> # Рисуем круговые диаграммы. Для того, чтобы площадь круга была пропорциональна
> # числу избирателей, радиус круга пропорционален квадратному корню из их числа.
> draw.pie(x, y, part, radius=sqrt(tik$v10)/1500, col=color)
> # Подпишем названия районов (мелким шрифтом и справа от центров «пирогов»)…
> text(x[grep("район", tik$x)]+0.05, y[grep("район", tik$x)],
+ labels=tik$x[grep("район", tik$x)], adj=0, cex=.5)
> # … и названия городов (нормальным шрифтом и слева от центров «пирогов»)
> text(x[grep("район", tik$x, invert=T)]-0.05, y[grep("район", tik$x, invert=T)],
+ labels=tik$x[grep("район", tik$x, invert=T)], adj=1)
> # Разместим легенду, заголовок и закроем png-устройство:
> legend("left", legend=c("ЕР", "КПРФ", "ЛДПР", "СР", "прочие"), fill=color, bty="n")
> title(main="Результаты выборов в Государственную Думу\nв муниципальных районах
+ и городских округах Удмуртии")
> dev.off()

UPD — Способ, позволяющий присоединять шкалу численности проголосовавших:

> # Сформируем матрицу с результатами партий, а не их долями:
> part <- matrix(c(tik$ER, tik$v10-tik$ER), ncol=2,
+ dimnames=list(tik$x ,c("ЕР", "прочие")))
> # Закажем для партий цвета:
> color <- c("#0000ff88", "transparent")
> plot(udm[udm$ADMIN_LVL == 2 & !is.na(udm$ADMIN_LVL),], col="white", border="gray")
> # Рисуем круговые диаграммы. Такой синтаксис обеспечивает пропорциональность
> # площади круга числу проголосовавших избирателей.
> draw.pie(x, y, part, radius=1/4, col=color)
> # Легенда для числа проголосовавших:
> legend.bubble("right", # Указываем размещение легенды, …
+ maxradius=1/4, # радиус самой большой окружности, …
+ z=round(max(tik$v10)/1000,0), # диапазон числовых подписей, …
+ txt.cex=0.5, # их размер, …
+ n=3, # и число…
+ bty="n") # а также отменяем отрисовку рамки.
> # Легенда для долей партий и их цветов:
> legend.pie("left", # Указываем размещение легенды, …
+ radius=1/4, # размер пирога, …
+ z=c(sum(tik$ER), sum(tik$v10-tik$ER)), # изображаемые на нём пропорции, …
+ col=color, # их цвета, …
+ labels=c("ЕР", "прочие"), # метки…
+ cex=0.5, # и их размер.
+ mab=2, # Корректируем отступы
+ bty="n") # а также отменяем отрисовку рамки.
> title(main="Результаты выборов в Государственную Думу\nв муниципальных районах
+ и городских округах Удмуртии")
> dev.off()

Рисуем фоновые картограммы в R: красим районы

Продолжение. Начало здесь: Рисуем картограммы в R.

Любоваться на нарисованную нами в R карту — занятие малополезное, на ней даже надписей никаких нет. Приведу более утилитарный пример:

> # Отберём муниципальные районы и городские округа…
> udm.TIK <- udm[udm$ADMIN_LVL == 6 & !is.na(udm$ADMIN_LVL), ]
> # … и сохраним их наименования в csv-файл.
> write.csv(udm.TIK$NAME, file="tik.csv", quote = FALSE)

Выглядит файл tik.csv так:

,x
1,Сюмсинский район
2,Кизнерский район
…
30,Камбарский район

Теперь его можно открыть любимым редактором таблиц и добавить, допустим, число голосов на выборах в Государственную Думу.

,x,v10,ER,SR,KPRF,LDPR
1,Сюмсинский район,6753,3628,702,1155,1104
2,Кизнерский район,10834,6607,937,1586,1458
…
30,Камбарский район,9605,4266,1654,1636,1801

Используем его для раскраски карты:

> tik <- read.csv("~/tik.csv")
> brks <- seq(from =.3, to = .7, by = .1) # Закажем границы цветовых интервалов…
> colors <- rev(heat.colors(length(brks))) # … и цвета для них.
> part <- tik$ER/tik$v10 # Определим долю голосов избирателей, поданых за ЕР.
> png("ER.png", width=5, height=5, units="in", res=100, bg="transparent")
> par(mai=c(.1,.1,1,.1))
> plot(udm.TIK, col=colors[findInterval(part, brks, all.inside=T)], border="white")
> legend("right", legend=leglabs(brks, under="<", over=">"), fill=colors, bty="n")
> title(main="Доля голосов за ЕР\nв муниципальных районах
+ и городских округах Удмуртии")
> dev.off()

Карту результатов КПРФ на тех же выборах мы нарисовали по аналогии, только цветовую палитру заказывали при помощи библиотеки RColorBrewer, вот так:

> library("RColorBrewer")
> colors <- brewer.pal(length(brks), "Reds")

Шейп-файлы с картой Удмуртии из проекта OpenStreetMap (лицензия CC-by-SA) предоставлены GIS-Lab.

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

Шейп-файлы с картой России (указаны границы административно-территориального деления вплоть до муниципальных районов и городских округов) можно получить на этом сайте сообщества GIS-Lab. Отдельные шейп-файлы для регионов РФ и стран — бывших республик СССР GIS-Lab выкладывает сюда (в файлах boundary-* указаны границы АТО вплоть до районов городских округов и территориальных органов управления).

> # Закажем необходимую для работы библиотеку maptools…
> library("maptools")
> #… и загрузим шейп-файл АТО Удмуртии в объект udm:
> udm <- readShapePoly("~/projects/r-maps/RU-UD/boundary-polygon.shp",
+ proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

Естественно, вместо "~/projects/r-maps/RU-UD/boundary-polygon.shp" следует указать имя Вашего файла, а вместо "+proj=longlat +ellps=WGS84 +datum=WGS84" — Ваши предпочтения по проекции и системе координат (Proj.4 представление). Значение proj4string можно изменять без перезагрузки шейп-файла, например, так:

> proj4string(udm) <- CRS("+proj=merc")

Собственно, уже можно начинать рисовать.

> # Создаём файл-устройство для рисования:
> png("boundary-all.png", width=300, height=400, units="px", bg="transparent")
> # Задаём поля отступов:
> par(mar=c(0.5,0.5,0.5,0.5))
> # Рисуем все АТО в объекте udm:
> plot(udm, col="white")
> # Закрываем файл-устройство:
> dev.off()

Результат:

Естественно, отображением можно управлять. Посмотрим на объект udm подробнее:

> summary(udm)
Object of class SpatialPolygonsDataFrame
Coordinates:
       min      max
x 51.12416 54.43483
y 55.85058 58.54584
Is projected: FALSE 
proj4string :
[+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs]
Data attributes:
     OSM_ID                    NAME       ADMIN_LVL  
 Min.   : -2076064   Октябрьское :  3   8      :311  
 1st Qu.: -1078857   Первомайское:  3   6      : 30  
 Median : -1066632   Каменское   :  2   9      :  5  
 Mean   :   521079   Камское     :  2   10     :  1  
 3rd Qu.: -1017007   Кильмезское :  2   2      :  1  
 Max.   :144692757   (Other)     :341   (Other):  2  
                     NA's        :  4   NA's   :  7

Нам интересны атрибуты OSM_ID, NAME и ADMIN_LVL. Так будет выглядеть та же карта с точностью до муниципальных районов и городских округов:

> png("boundary-6.png", width=300, height=400, units="px", bg="transparent")
> par(mar=c(0.5,0.5,0.5,0.5))
> plot(udm[udm$ADMIN_LVL == 6 & !is.na(udm$ADMIN_LVL), ], col="white")
> dev.off()

Конструкция "!is.na(udm$ADMIN_LVL)" отбрасывает те объекты, у которых не указан административный уровень. Если этого не делать, получим ошибку "NAs not permitted in row index". Запятая перед закрывающей квадратной скобкой (выделена красным) позволяет избежать ошибки "undefined columns selected".

Отбор отображаемых объектов можно делать и по их имени, вот так:

> plot(udm[udm$NAME == "Малопургинский район" & !is.na(udm$NAME), ], col="white")

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