9 мар. 2014 г.

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)] # (предварены символом "_").

5 нояб. 2013 г.

Картограммы в 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()

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

3 нояб. 2013 г.

Решение задачи коммивояжера в R

Для решения задачи коммивояжера в R-project можно использовать пакет TSP:

> library("TSP")
# Заполним матрицу расстояний между городами…
> m <- matrix(c(0,   222, 652, 52,  447, 160, 445, 134, 380, 112, 536, 
+               222, 0,   594, 170, 565, 278, 223, 252, 322, 230, 654,
+               652, 594, 0,   600, 995, 708, 817, 682, 272, 660, 1084,
+               52,  170, 600, 0,   395, 108, 393, 82,  328, 60,  484,
+               447, 565, 995, 395, 0,   503, 340, 313, 723, 455, 519,
+               160, 278, 708, 108, 503, 0,   501, 190, 436, 48,  592,
+               445, 223, 817, 393, 340, 501, 0,   475, 545, 453, 877,
+               134, 252, 682, 82,  313, 190, 475, 0,   410, 142, 402,
+               380, 322, 272, 328, 723, 436, 545, 410, 0,   388, 812,
+               112, 230, 660, 60,  455, 48,  453, 142, 388, 0,   544,
+               536, 654, 1084, 484, 519, 592, 877, 402, 812, 544, 0),
+             nrow=11, byrow=TRUE)
> colnames(m) <- c("Воткинск", "Глазов", "Екатеринбург", "Ижевск", "Казань",
+                  "Камбарка", "Киров", "Можга", "Пермь", "Сарапул", "Уфа")
> rownames(m) <- colnames(m)
> m
             Воткинск Глазов Екатеринбург Ижевск Казань Камбарка Киров Можга Пермь Сарапул  Уфа
Воткинск            0    222          652     52    447      160   445   134   380     112  536
Глазов            222      0          594    170    565      278   223   252   322     230  654
Екатеринбург      652    594            0    600    995      708   817   682   272     660 1084
Ижевск             52    170          600      0    395      108   393    82   328      60  484
Казань            447    565          995    395      0      503   340   313   723     455  519
Камбарка          160    278          708    108    503        0   501   190   436      48  592
Киров             445    223          817    393    340      501     0   475   545     453  877
Можга             134    252          682     82    313      190   475     0   410     142  402
Пермь             380    322          272    328    723      436   545   410     0     388  812
Сарапул           112    230          660     60    455       48   453   142   388       0  544
Уфа               536    654         1084    484    519      592   877   402   812     544    0
# … создадим TSP-объект…
> tsp <- TSP(m)
# … и найдем решение методом "2-opt":
> tour <- solve_TSP(tsp, method="2-opt")
> tour
object of class ‘TOUR’ 
result of method ‘2-opt’ for 11 cities
tour length: 3080
# Получаем замкнутый маршрут, заданный списком городов по порядку обхода.
> labels(tour)
 [1] "Камбарка"     "Сарапул"      "Ижевск"       "Екатеринбург" "Пермь"        "Глазов"      
 [7] "Киров"        "Казань"       "Уфа"          "Можга"        "Воткинск"
> 48+60+600+272+322+223+340+519+402+134+160
[1] 3080

Существует альтернативный способ создания таблицы расстояний, с заполнением только нижней ее половины:

> m <- structure(c(222, 652, 52,  447, 160, 445, 134, 380, 112, 536,
+                  594, 170, 565, 278, 223, 252, 322, 230, 654,
+                  600, 995, 708, 817, 682, 272, 660, 1084,
+                  395, 108, 393, 82,  328, 60,  484,
+                  503, 340, 313, 723, 455, 519,
+                  501, 190, 436, 48,  592,
+                  475, 545, 453, 877,
+                  410, 142, 402,
+                  388, 812,
+                  544), 
+                Labels = c("Воткинск", "Глазов", "Екатеринбург", "Ижевск", "Казань",
+                           "Камбарка", "Киров", "Можга", "Пермь", "Сарапул", "Уфа"), 
+                 Size = 11L, class = "dist", Diag = FALSE, Upper = FALSE)
# Результат:
> m
             Воткинск Глазов Екатеринбург Ижевск Казань Камбарка Киров Можга Пермь Сарапул
Глазов            222                                                                     
Екатеринбург      652    594                                                              
Ижевск             52    170          600                                                 
Казань            447    565          995    395                                          
Камбарка          160    278          708    108    503                                   
Киров             445    223          817    393    340      501                          
Можга             134    252          682     82    313      190   475                    
Пермь             380    322          272    328    723      436   545   410              
Сарапул           112    230          660     60    455       48   453   142   388        
Уфа               536    654         1084    484    519      592   877   402   812     544

Отличий в использовании такой таблицы нет.

31 мар. 2013 г.

Живой журнал

Завёл себе живой журнал. Уже давно — 16 сентября 2012 года.

30 мар. 2013 г.

Картограммы в 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) # Минимальная площадь сохранямых островов.

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

24 нояб. 2012 г.

Рисуем точечные картограммы в 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.

22 нояб. 2012 г.

Рисуем картодиаграммы в 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()