12 февр. 2013 г.

Построение страновых хороплет-карт (choropleth maps) с помощью пакета rworldmap в статистическом пакете R


Довольно часто социально-экономические данные удобно визуализировать с помощью карт. Стандартным является представление в форме так называемых choropleth maps (похоже, на русском они так и называются - хороплет-карты). Хороплет - это такая карта, на которой с помощью визуального оформления (обычно цвета, но может быть и к примеру, точками или линиями разной густоты/толщины) отображается интенсивность какого-либо показателя для географических регионов. Такие карты очень легко позволяют наглядно отобразить данные по большому количеству географических объектов.
Вот для примера такая карта - правда, наглядно? На ней изображена доля курильщиков среди мужчин по почти 150 странам. Создание такой карты в R заняло 8 строк кода, включая загрузку данных из Сети.
Smoking prevalence, males (% of adults)

Как создавать подобные карты?
Существует множество различных вариантов. В принципе, любая ГИС-программа типа MapInfo или ArcGIS умеет строить хороплет, но покупать эти дорогостоящие системы только для того, чтобы рисовать цветные карты, смысла не имеет. Можно рисовать карты хорорлет и в любимом Excel. Существует несколько "рецептов", как можно это делать.  Как правило, они требуют запуска VBA-макроса - для того, чтобы сопоставить исходный показатель с географическими регионами на карте. Еще одно условние - необходимо иметь исходную карту (линии) в графическом формате SVG. Это обычно не проблема, и другие гео-форматы можно тоже преобразовать в графический SVG, но это требует определенных манипуляций Вот один из предлагаемых рецептов, как построить карту хорорлет в Excel. Обратите внимание на то, сколько телодвижений надо сделать для того, чтобы нарисовать достаточно простую карту. Причем если даже процесс автоматизирован при смене данных (то есть, меняешь данные - карта перерисовывается), то при смене самой карты все ломается и требуется "доводка напильником" под новую карту.
 Наконец - обещанное. Конечно, хороплет-карты можно построить в статистическом пакете R. Возможности по работе с гео-данными в R, в том числе и графические, сопоставимы, а в некоторых случаях даже превосходят возможности полноценных ГИС-систем. И все это - в одном флаконе.

Сейчас я расскажу о простом и часто встречающемся случае - визуализация данных по странам как в случае нарисованного выше графика. Как обычно в R, для каждой специальной задачи уже есть созданный кем-то пакет. В нашем случае мы воспользуемся пакетом - rworldmap, созданный и поддерживамый Andy South. Основная его задача как раз в этом - построение хороплет-карт для страновых данных. Фактически этот пакет  - надстройка над другими, которые обеспечивают работу с гео-данными, но нас это не должно слишком волновать. Для нас пакет rworldmap сильно экономит время и строки для построения типовых страновых карт.

Краткая процедура построения выглядит следующим образом:

  1. Получение данных для отображения. В нашем случае мы используем данные из базы данных World Development Indicators (WDI), которую поддерживает Всемирный банк. В этой базе, как обещано, 331 показатель по 221 стране. Возьмем показатель "Smoking prevalence, males (% of adults)", который мне показался достаточно интересным. Загрузим его в R с помощью специального пакета, который носит название "WDI", как и следовало ожидать. Этот пакет позволяет обратиться напрямую к базе данных из R и получить данные сразу в необходимой разбивке "страна"-"год"-"показатель". Можно сразу загружать несколько показатель для диапазона лет для всех или некоторых стран. Очень удобно! В функции WDI() из одноименного пакета я рекомендую использовать параметр extra = TRUE для того, чтобы получаемые данные содержали код страны по ИСО (двух и трехзначный). 
  2. Сопоставление карты и данных. Для этого уже используется пакет rworldmap и его функция с говорящим названием joinCountryData2Map. Здесь все просто выбираем ключ, по которосу мы сопоставляем данные и карту. Это может быть название страны на английском языке, но лучше всего использовать коды стран по ISO (ISO2 или ISO3). Соответственно это параметр "joinCode". 
  3. Выбор показателя, который мы хотим отобразить на карте с помощью функции mapCountryData. Здесь все понятно  nameColumnToPlot соответственно название показателя, который хотим отобразить. Важными являются также параметры catMethod (определяет методы для категоризации данных, то есть деления их на группы одного цвета), а также colourPalette, который определяет цветовую палитру. Для палитры удобно использовать имеющиеся варианты из пакета RColorBrewer. Как выглядят сами палитры удобно смотреть на специальном сайте
  4. Дополнительно. Функция mapCountryData сразу же отображает карту, но если вы хотите внести некоторые изменения - в нашем случае немного поправить легенду с помощью addMapLegend. 
Вот код , который уже нарисует долю курильщиков-женщин по странам мира:

<!-- Styles for R syntax highlighter
require(WDI)
require(rworldmap)
require(RColorBrewer)

smoker <- WDI(country = "all", indicator = "SH.PRV.SMOK.FE", start = 2009, end = 2009, 
    extra = TRUE, cache = NULL)
# получить данные из WDI
fMap <- joinCountryData2Map(smoker, joinCode = "ISO3", nameJoinColumn = "iso3c", 
    nameCountryColumn = "country")
# cопоставить данные с картой по коду ISO3 (трехзначный код страны)
mapParams <- mapCountryData(fMap, nameColumnToPlot = "SH.PRV.SMOK.FE", catMethod = "quantiles", 
    missingCountryCol = gray(0.8), colourPalette = rev(brewer.pal(7, "RdYlBu")), 
    addLegend = FALSE, mapTitle = "")
# добавить легенду
do.call(addMapLegend, c(mapParams, legendLabels = "all", labelFontSize = 0.8, 
    legendWidth = 0.5, sigFigs = 0))

--> А вот и результат его выполнения:
Smoking prevalence, females (% of adults)
По картам можно сразу же сделать несколько наблюдений. Среди мужчин больше всего курят Россия, Китай, Восточная Европа и Малайзия с Индонезией. Среди женщин все по другому - здесь в лидерах жительницы стран Западной Европы (Австрия на первом месте), Чили и  Папуа-Новая Гвинея.