Создание картограмм при помощи R

Одним из интересных и весьма актуальных аспектов работы с данными является их визуализация в привязке к географическим координатам. В программе R имеется богатый арсенал инструментов для работы с пространственными данными (подробный обзор можно найти на странице CRAN Task Views: Analysis of Spatial Data). Например, в одном из предыдущих сообщений я уже писал о создании пользовательских карт Google средствами пакета googleVis. В этом сообщении я покажу, как можно создать картограммы (англ. choropleth maps) при помощи R.

На картограммах отдельные территориальные единицы залиты цветом, интенсивность которого соответствует величине отображаемой количественной переменной (плотность населения, уровень дохода, процент проголосовавших за  определенного кандидата на выборах, и т.п.). В R подобные карты можно построить несколькими способами. Ниже будет рассмотрен один из них - при помощи пакета sp. Для воспроизведения примеров просто повторяйте описанные ниже шаги.





Шейп-файлы

Все карты мы будем строить на примере Республики Беларусь. Для начала нам потребуются т.н. "шейп-файлы" (англ. shapefile) для разных административно-территориальных единиц страны. Шейп-файл представляет собой распространенный векторный формат хранения информации о геометрическом положении и определенных атрибутах географических объектов. Несмотря на свое название, в действительности шейп-файл - это не один, а набор из нескольких файлов с одинаковым именем, но разными расширениями. Обязательными при этом являются файлы с расширениями .shp, .dbf и .shx. Я не буду описывать здесь, что содержится в каждом из этих файлов - об этом можно найти много информации в Сети (см., например, статью в Википедии и официальное техническое описание формата). Для наших целей достаточным будет представление о шейп-файле как о файле, по данным из которого можно воспроизвести контурную карту страны и/или отдельных ее территориальных единиц.

Естественно, встает вопрос о том, где можно раздобыть шейп-файлы для той или иной интересующей нас территории? Одним из бесплатных (для некоммерческого использования!)  источников является Global Administrative Areas (GADM) - Глобальная база данных административных областей. В этой базе данных доступны шейп-файлы трех разных уровней пространственного разрешения, например, на уровне страны, областей и районов. Соответствующие файлы для Беларуси на сайте GADM можно получить двумя разными путями. Во-первых, мы можем сделать это "вручную", зайдя на сайт GADM и выполнив следующие шаги.

1. В главном меню сайта выберите опцию Download (Скачать; отмечена стрелкой на рисунке ниже):

2. На открывшейся странице из выпадающих меню выбрать желаемую страну (Country) и формат скачиваемых файлов (File Format).

В меню File Format, как видно из рисунка ниже, доступны несколько опций, включая нужный нам шейп-файл (Shapefile). Обратите внимание: народ, работающий над созданием и поддержкой GADM, отлично осведомлен о распространенности и возможностях R, и, как следствие, в качестве одной из опций в меню File Format мы видим R (SpatialPolygonsDataFrame) - файл "рабочего пространства" (англ. workspace), который уже содержит шейп-файл в необходимом формате и может быть непосредственно загружен в R. Описание того, как это сделать, приведено ниже. Пока же выберем опцию Shapefile:

3. На открывшейся странице увидим карту страны, под которой находится ссылка download для скачивания zip-архива с нужными нам шейп-файлами для всех трех уровней пространственного разрешения. Сохраните этот архив (BLR_adm.zip) в удобном для вас месте.

4. Разархивируйте BLR_adm.zip в текущую рабочую папку R. В результате в этой папке появятся следующие файлы:

Все перечисленные файлы важны и должны находиться в одной директории. Обратите внимание на цифры в именах файлов - они помогают понять, какому уровню пространственного разрешения соответствует каждый файл. Так, файлы, в названии которых имеется 0, описывают контуры страны, 1 - контуры страны и областей, а 2 - контуры страны, областей и районов.

Для загрузки шейп-файлов следует воспользоваться функцией readShapePoly() из пакета maptools (установите его предварительно при помощи команды install.packages("maptools")). Так, шейп-файл для первого уровня пространственного разрешения (страна, разделенная на области) можно загрузить следующим образом:

library(maptools)
Regions <- readShapePoly("BLR_adm1.shp")

Второй способ получения шейп-файлов с сайта GADM состоит в их загрузке по URL непосредственно из R. Я предпочитаю именно этот способ и ниже буду предполагать, что им вы и воспользовались. Так, для загрузки шейп-файла с контурами страны и областей следует выполнить следующую команду:

load(url("http://gadm.org/data/rda/BLR_adm1.RData"))

В результате в рабочей среде R появится объект с именем gadm. Во избежание путаницы в дальнейшем, стоит сохранить этот объект, присвоив ему какое-либо более описательное имя, например Regions, а исходный объект затем просто удалить:

Regions <- gadm
rm(gadm)

С точки зрения своего внутреннего устройства и представления в R, созданный нами объект Regions принадлежат к т.н. объектам класса S4 (подробнее об этом классе можно узнать здесь (англ. яз.)). Вкратце можно сказать, что такие объекты имеют сложную структуру, напоминающую структуру обычных списков R. Отдельные компоненты S4-объектов называются слотами (англ. slots) и могут сами являться сложными объектами любого другого класса. В отличие от обычных списков, обращение к слотам S4-объектов происходит при помощи оператора @, а не $. Кроме того, к слотам можно обращаться только по их именам, тогда как отдельные элементы списков мы можем индексировать также по их порядковым номерам. Имена слотов, входящих в состав Regions легко выяснить при помощи базовой R-функции slotNames():

slotNames(Regions)
[1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

Наибольший интерес для нас представляют слоты "data" и "polygons". Первый включает обычную таблицу данных (data frame), содержание которой станет понятным из приведенного ниже вывода команды Regions@data (примечание: на моем компьютере русские названия территориальных единиц не распознались - отсюда знаки ??? в примере ниже):

Regions@data
  ID_0 ISO  NAME_0 ID_1    NAME_1 NL_NAME_1
0   21 BLR Belarus    1     Brest      <NA>
1   21 BLR Belarus    2   Homyel'      <NA>
2   21 BLR Belarus    3    Hrodna      <NA>
3   21 BLR Belarus    4  Mahilyow      <NA>
4   21 BLR Belarus    5     Minsk      <NA>
5   21 BLR Belarus    6 Vitsyebsk      <NA>
                                               VARNAME_1
0                       ?????|Brestskaya Voblasts'|Brиst
1       ??????|Gomel|Gomel'|Homyel'skaya Voblasts'|Homje
2                   ??????|Grodno|Hrodzenskaya Voblasts'
3       ???????|Mahiljow|Mogilev|Mahilyowskaya Voblasts'
4  ?????|??????? ???????|Minsk Oblast|Minskaya Voblasts'
5 ???????|Vicebsk|Vitebsk|Vitsyebskaya Voblasts'|Witebsk
     TYPE_1 ENGTYPE_1
0 Voblasts'    Region
1 Voblasts'    Region
2 Voblasts'    Region
3 Voblasts'    Region
4 Voblasts'    Region
5 Voblasts'    Region

В свою очередь слот "polygons" содержит список с координатами точек, задающих контуры соответствующих территориальных единиц (в формате системы координат WGS84), а также некоторые другие атрибуты:

str(Regions@polygons)
List of 6
 $ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. ..@ Polygons :List of 1
  .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. ..@ labpt  : num [1:2] 25.5 52.4
  .. .. .. .. ..@ area   : num 4.3
  .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. ..@ coords : num [1:2100, 1:2] 25.9 25.9 26 26 26 ...
  .. ..@ plotOrder: int 1
  .. ..@ labpt    : num [1:2] 25.5 52.4
  .. ..@ ID       : chr "0"
  .. ..@ area     : num 4.3
... # остальная часть выведенной информации опущена для экономии места

Разобравшись немного со структурой шейп-файлов и их представлением в R, приступим собственно к построению картограмм.


Функция spplot() из пакета sp

Пакет sp содержит функции, реализующие ряд методов и классов для пространственных данных (импорт, преобразование и экспорт данных, графические методы, и т.п.). В частности, в нем имеется функция spplot(), которая позволяет отображать на графике данные, заключенные в шейп-файле. В основе этой функции лежат высокоуровневые команды lattice - одного из популярных графических R-пакетов.  Не забудьте установить sp перед его первым использованием (install.packages("sp")).

Для начала построим картограмму Беларуси, на которой каждая область будет залита своим цветом (цвет на этой карте будет служить просто для идентификации области - никакой другой информационной нагрузки он нести не будет). Посмотрим еще раз на названия областей (приведены оригинальные белорусские названия в написании латинскими буквами):

Regions@data$NAME_1
[1] Brest     Homyel'   Hrodna    Mahilyow  Minsk     Vitsyebsk
Levels: Brest Homyel' Hrodna Mahilyow Minsk Vitsyebsk

Таким образом, первой идет Брестская область, а за ней - Гомельская, Гродненская, Могилевская, Минская и Витебская области. Переменная NAME_1, как видим, является фактором с 6 уровнями. Каждому их этих уровней при построении карты необходимо будет присвоить уникальный цвет. Хорошо отличимые друг от друга цвета можно подобрать, воспользовавшись одной из стандартных палитр R, например rainbow():

spplot(Regions,
       "NAME_1", # отображаемая переменная
       scales = list(draw = T), # опция, включающая отображение координатных осей
       col.regions = rainbow(n = 6) ) # опция, задающая заливку цветом

Можно попробовать и другие стандартные палитры R, например:

Однако я больше предпочитаю хорошо подобранные и потому приятные на вид палитры, входящие в состав пакета RColorBrewer (см. также здесь). Помимо использования одной из палитр RColorBrewer, в приведенном примере также добавлена опция par.settings, при помощи которой отключаются координатные оси вокруг карты:

spplot(Regions, "NAME_1",
       col.regions = brewer.pal(6, "Set3"), 
       par.settings = list(axis.line = list(col = NA)))

Теперь придадим цвету смысл - построим картограмму, на которой цветовая шкала будет соответствовать числу мужчин, проживающих в каждой области. Данные взяты из этого отчета (стр. 10), составленного по результатам переписи населения Беларуси за 2009 г. Значения численности мужчин достаточно добавить в таблицу данных, хранящуюся в слоте data объекта Regions (см. выше):

# Порядок чисел важен! Он должен совпадать с порядком, в котором
# таблице data перечислены соответствующие области:
Regions@data$Population = c(186716, 152169, 131817,
                            105417, 253427, 135348)

Перед построением картограммы нам необходимо определиться с цветовой шкалой, т.е. какими цветами шкала должна быть представлена и на сколько градаций разбита. Сделать это можно разными способами, в зависимости от стоящей задачи, а также от разброса значений отображаемой на картограмме переменной. Я воспользуюсь одним из таких способов - определением цветовой палитры при помощи базовой R-функции colorRampPalette(). В качестве аргументов этой функции следует указать "опорные" цвета, между которыми требуется произвести интерполяцию (т.е. плавный переход от, например, красного к желтому). Функция colorRampPalette() возвращает другую функцию, которая и выполняет такую интерполяцию автоматически. Для нашего примера мы создадим цветовую шкалу, в которой будет происходить плавный переход от цвета seagreen к whitesmoke (для ознакомления со списком стандартных цветов R выполните команду colors() ):

mypalette <- colorRampPalette(c("seagreen", "whitesmoke"))
 
# Объект mypalette - это функция:
mypalette
function (n) 
{
    x <- ramp(seq.int(0, 1, length.out = n))
    rgb(x[, 1L], x[, 2L], x[, 3L], maxColorValue = 255)
}

Теперь функцию mypalette мы можем подать на функцию spplot(). При этом в качестве (единственного) аргумента mypalette следует указать общее количество цветов (n, опорные плюс промежуточные), которое должно присутствовать в шкале. Чем больше разброс значений отображаемой переменной, тем выше должно быть это число - так будет достигаться плавный переход между цветами. В нашем случае разброс значений численности мужчин в отдельных областях невелик, и поэтому я выбрал n = 20. В приведенной ниже команде для демонстрации еще одной опции функции spplot() добавлен аргумент col со значением "transparent" (прозрачный) - с его помощью отключается отображение линий, ограничивающих каждую область:

spplot(Regions, "Population",
       col.regions = mypalette(20), # определение цветовой шкалы
       col = "transparent", # отключение контурных линий на карте
       par.settings = list(axis.line = list(col = NA)))

Как видим, наибольшая численность мужчин в Беларуси имеет место в Минской области и в двух южных областях - Брестской и Гомельской.

Рассмотренный принцип добавления пользовательских данных в слот data шейп-файла c последующим их отображением на картограмме можно применить и к более дробным территориальным единицам - например, районам. Загрузим соответствующий шейп-файл с сайта GADM и сохраним его в виде объекта с именем Counties:

load(url("http://gadm.org/data/rda/BLR_adm2.RData"))
Counties <- gadm
rm(gadm)

Вместо каких-либо реальных данных, я воспользуюсь генератором случайных чисел и создам вектор Fake из нормально распределенных значений. Этот вектор будет добавлен в слот data объекта Counties. Длина вектора Fake будет соответствовать числу районов в Беларуси:

set.seed(1234) # для воспроизводимости результата
Counties@data$Fake <- rnorm(dim(Counties@data)[1], 100, 35)

Для создаваемой картограммы я применю еще один способ определения цветовой шкалы. В основу этой шкалы лягут цвета из пакета RColorBrewer (см. выше), в частности из палитры RdBu (см. сайт Colorbrewer). При работе с количественными переменными возникает проблема с тем, чтобы разбить их значения на отдельные промежутки, которым можно сопоставить соответствующие цвета. В R такое разбиение не составляет труда, причем есть несколько способов это сделать. Мы воспользуется возможностями пакета classInt (не забудьте установить его перед первым использованием!):

require(classInt)
 
# Функция classIntervals разбивает имеющуюся совокупность чисел на
# классы с определенным шагом, задаваемым при помощи аргумента
# style. Мы постараемся разбить наш вектор Fake на 6 классов одинакового размера:
brks.eq = classIntervals(Counties$Fake, n = 6, style = "equal")
 
# Вот результат разбиения этого вектора на 6 одинковых классов:
brks.eq
style: equal
  one of 167,549,733 possible partitions of this variable into 6 classes
[17.90058,46.45293) [46.45293,75.00528) [75.00528,103.5576) 
                  4                  31                  43 
  [103.5576,132.11)   [132.11,160.6623) [160.6623,189.2147] 
                 24                  11                   5
 
# Задаем палитру из 6 цветов:
plotcol <- brewer.pal(6,"RdBu")
 
# Рисуем карту:
spplot(Counties, "Fake",
       col.regions = plotcol,
       at = brks.eq$brks, # задает границы классов
       par.settings = list(axis.line = list(col = NA)))

В следующем сообщении я опишу, как можно создать картограммы при помощи другого пакета R - популярного ggplot2.