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

# Скачиваем файл
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=='Псковская область', ]

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

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

# получаем координаты центроида каждого района,
# чтобы разместить названия районов
# пользуемся функцией 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

# Читаем массив с данными по площади и населению районов
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)

# Функция по "прилизываю" графика.
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')