Press "Enter" to skip to content

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

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

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

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

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

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

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

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

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

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

1
2
3
4
5
6
7
8
9
10
11
# Скачиваем файл
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')

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

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

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

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

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

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

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

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

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

28
29
30
31
32
33
# Раздвигаем метки по вертикали.
# Координаты для сдвига получены опытным путём.
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-файл.

34
35
36
37
38
39
40
41
# Читаем массив с данными по площади и населению районов
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

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

42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
# Функция по "прилизываю" графика.
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)
}

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

60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
# Создаём переменную класса 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')

4 Comments

  1. Сергей
    Сергей 2014-11-12

    Добрый день, Александр!
    R выдает ошибки при загрузке файла
    download.file(‘http://gis-lab.info/data/yav/adm/20120101/adm6_district.7z’, ‘adm6_district.7z’)
    и установки пакета rgdal.
    При загрузке файла пишет:
    download.file(‘http://gis-lab.info/data/yav/adm/20120101/adm6_district.7z’, ‘adm6_district.7z’)
    пробую URL ‘http://gis-lab.info/data/yav/adm/20120101/adm6_district.7z’
    Ошибка в download.file(“http://gis-lab.info/data/yav/adm/20120101/adm6_district.7z”, :
    не могу открыть URL ‘http://gis-lab.info/data/yav/adm/20120101/adm6_district.7z’
    Вдобавок: Предупреждение
    In download.file(“http://gis-lab.info/data/yav/adm/20120101/adm6_district.7z”, :
    не могу открыть: статус HTTP был ‘404 Not Found’

    При загрузке пакета пишет:
    > install.packages(“rgdal”)
    Installing package into ‘/home/serg/R/i686-pc-linux-gnu-library/3.1’
    (as ‘lib’ is unspecified)
    пробую URL ‘http://cran.rstudio.com/src/contrib/rgdal_0.9-1.tar.gz’
    Content type ‘application/x-gzip’ length 1624496 bytes (1.5 Mb)
    открытие URL
    ==================================================
    downloaded 1.5 Mb

    * installing *source* package ‘rgdal’ …
    ** пакет ‘rgdal’ удачно распакован, MD5 sums проверены
    configure: CC: gcc -std=gnu99
    configure: CXX: g++
    configure: rgdal: 0.9-1
    checking for /usr/bin/svnversion… no
    configure: svn revision: 518
    configure: gdal-config: gdal-config
    checking gdal-config usability…
    ./configure: line 2119: gdal-config: command not found
    no
    Error: gdal-config not found
    The gdal-config script distributed with GDAL could not be found.
    If you have not installed the GDAL libraries, you can
    download the source from http://www.gdal.org/
    If you have installed the GDAL libraries, then make sure that
    gdal-config is in your path. Try typing gdal-config at a
    shell prompt and see if it runs. If not, use:
    –configure-args=’–with-gdal-config=/usr/local/bin/gdal-config’
    with appropriate values for your installation.

    ERROR: configuration failed for package ‘rgdal’
    * removing ‘/home/serg/R/i686-pc-linux-gnu-library/3.1/rgdal’
    Warning in install.packages :
    installation of package ‘rgdal’ had non-zero exit status

    Вы не могли бы помочь разобраться в этих ошибках.
    R 1.1.1 ос ubuntu 14.04

    • Александр Матрунич
      Александр Матрунич 2014-11-19

      Добрый день, Сергей!

      > Вдобавок: Предупреждение
      > In download.file(“http://gis-lab.info/data/yav/adm/20120101/adm6_district.7z”, :
      > не могу открыть: статус HTTP был ‘404 Not Found’

      Что-то поменялось на сайте gis-lab.info, и теперь этого файла нет. Надо смотреть, что там теперь у них. Где-то у них есть проект по административным границам, там должно быть написано, где теперь хранятся файлы.

      > Error: gdal-config not found
      Помимо самого gdal в системе должна быть установлена версия для разработки. Она нужна, чтобы скомпилировать расширение, которое зависит от этой программы. Поищите, что есть в Убунте: gdal-dev или gdal-devel, что-то в этом духе. Устновите и попробуйте заново. Для Федоры этого хватает.

      Не понял, а R какой версии стоит? 1.1.1?

  2. Сергей
    Сергей 2014-11-27

    Добрый день, Александр!

    На счет версии R, я ошибся, стоит 3.1.2.

    С установкой пакета разобрался, нужно было ввести сначала
    sudo apt-get install libgdal-dev
    потом только ставить пакет.

    На счет загрузки карты, я загрузил другую, по Пензенской обл.
    Спасибо!

Leave a Reply

Your email address will not be published. Required fields are marked *

Защита от спама *