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

Плотность населения районов Псковской области. Картограмма сделана в 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')