R: хороплет-карта России с увеличенной европейской частью

Коротко о главном: прочитал недавно пост infotanka. Полез на сайт Татьяны Мисютиной и подсмотрел там хороплет-карту России с увеличенной европейской частью. И ведь, действительно, классная идея. Удобно, наглядно. Захотелось сделать себе шаблон под R для таких же графиков. Ведь хорошие идеи должны размножаться делением?

С сайта Global Administrative Areas берем полигоны границ регионов России. Там как раз есть нативный для R формат. Данные загружаются в виде SpatialPolygonsDataFrame.

Загружаем файл, пробуем отрисовать, что есть:

rusdf <- load('RUS_adm1.RData’) plot(gadm) 

Получилась вот такая странная штука:

В принципе, все правильно. Никакие проекции мы не использовали, по координатам все похоже на правду. Но «чукотский разброс» смущает. Конечно, самым правильным решением было бы использование цилиндрической проекции к тем полигонам, что уже есть. Чукотка «соберется», но граница по меридиану 180 долготы останется и будет видна. И по острову Врангеля, и по Чукотскому АО. Видимо, из-за погрешностей в один полигон они не соединяются. Плохо. Поэтому было принято решение расширить диапазон долгот over 180 град. Ага, «число Пи в военное время может достигать четырех».

for(i in 1:length(gadm@polygons)){   for(j in 1:length(gadm@polygons[[i]]@Polygons)){     gadm@polygons[[i]]@Polygons[[j]]@coords[,1]<- sapply(gadm@polygons[[i]]@Polygons[[j]]@coords[,1], function(x){         if(x < 0){           x<-359.999+x         }         else{x}       })   } } 

Что-то уже вырисовывается. Но что это? Читинская область, Агинский Бурятский, Корякский и Коми-Пермяцкий автономные округа. Как там у нас год на дворе? 2003 или 2013? В общем, надо начинать процесс административного объединения субъектов РФ. Старые регионы тоже на всякий случай сохраним. Вдруг пригодится для визуализации архивных данных.

# Zabaikalsky krai (Chitinskaya obl. + Aginskiy Buryatskiy AOk) united.reg[united.reg == 2 | united.reg == 13] <- 91  united.reg <- as.character(united.reg) rus.map <- unionSpatialPolygons(gadm, united.reg)  # Returning old regions (before unioning) old.regions <- list() old.regions <- c(old.regions, Polygons(gadm[gadm$ID_1==2,]@polygons[[1]]@Polygons, ID = '2')) old.regions <- c(old.regions, Polygons(gadm[gadm$ID_1==13,]@polygons[[1]]@Polygons, ID = '13'))  rus.map <- SpatialPolygons(c(slot(rus.map,'polygons'), old.regions))  

После объединения регионов кое-где на месте старых границ из-за непересечения образуются артефакты. Фильтруем артефакты по малой площади полигонов и параметру hole (функция clean.borders). С подготовкой карты, кажется, все.

Можно грузить таблицу со статистикой и рисовать.

В таблице проставлены ID регионов, по которым полигоны будут объединяться со статистикой. Там есть и устаревшие регионы, и существующие сейчас. Выбор полигонов для отрисовки происходит по колонке с данными: если в поле статистики стоит «NA» — этот субъект Федерации не отрисовывается совсем; если статистика «нулевая» — регион отрисовывается серым цветом. Остальные цвета карты настраиваются с помощью палитр RColorBrewer. Тут же можно выбрать тип палитры: (1) sequential — оттенки одного цвета, показывающие выраженность признака или (2) diverging — отклонение от среднего значения в плюс (один цвет, напр. зеленый) или минус (другой цвет, напр. красный).

Следующий этап — создание базовой карты. Последующая итоговая картинка будет просто комбинацией двух видов на этот базовый график. Убираем из базового графика все отступы, фон, подписи осей и метки.

p <- p + theme(axis.line=element_blank(),axis.text.x=element_blank(),                axis.text.y=element_blank(),axis.ticks=element_blank(),                axis.title.x=element_blank(),                axis.title.y=element_blank(),                legend.position = 'none',                panel.margin = unit(c(0,0,0,0), 'cm'),                axis.ticks.margin = unit(0, 'cm'),                axis.ticks.length = unit(0.001, 'cm'),                plot.margin = unit(c(0,0,0,0), 'cm'),                panel.grid = element_blank(),                panel.background = element_blank()               ) p <- p + labs(x=NULL, y = NULL) 

И небольшой «хак», связанный с порядком отрисовки полигонов. Ggplot2 при построении полигонов из SpatialPolygons никак не учитывает их свойство hole. В итоге, если Московская область будет отрисовываться после Москвы, даже несмотря на то, что в исходных данных был «исключающий» полигон-дырка — ggplot закрасит территорию Москвы цветом области. Та же проблема с Адыгеей. Поэтому эти два региона перерисовываем в последнюю очередь, после того как построены все остальные.

Теперь из одного графика нам надо получить два представления: увеличенная европейская часть (p1) и остальная страна (p2). На этом же этапе выбираем проекцию карты (с равными расстояниями), положение «камеры» и угол обзора. На один из графиков возвращаем легенду. Прелесть функции coord_map в том, что при отрисовке данные, не попавшие на картинку, не отфильтровываются и учитываются при построении видимой части графика. То есть, даже если диапазон значений статистики для европейской и азиатской частей России будет различаться — в легенде и цветовом кодировании сбоев не будет. Поля и углы обзора для этого графика я подбирал «методом научного тыка», перебирая различные варианты. Буду рад если кто-то подскажет более надежный и повторяемый способ.

Сборка финальной картинки делается при помощи библиотеки grid. Она позволяет разметить области для вставки отдельных элементов в общий рисунок (viewport’ы). Пока без «тяжелой» графики, с помощью placeholder-ов прикидываем подходящий вариант размещения элементов:

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

magnif.glass <- function(vport){   grid.circle(x=.6,y=.6,r=.3, gp=gpar(lwd=1.5, col='grey70'), vp = vport)   grid.lines(x=c(.6,.6), y=c(.5,.7), gp=gpar(lwd=1.5, col='grey70'), vp = vport)   grid.lines(x=c(.5,.7), y=c(.6,.6), gp=gpar(lwd=1.5, col='grey70'), vp = vport)   grid.lines(x=c(.1,.4), y=c(.1,.4), gp=gpar(lwd=1.5, col='grey70'), vp = vport)   grid.lines(x=c(.1,.3), y=c(.1,.3), gp=gpar(lwd=3, col='grey70'), vp = vport) }  

И добавляем текстовую выноску со средним значением или какой-нибудь другой огромной и важной цифрой. Вот вроде бы и все.

Весь код с тестовыми данными лежит на GitHub. В конце два вопроса сообществу:
1. Где взять shape-файл границ новой Москвы? Хочется проапдейтить этот момент.
2. Что точно надо переделать/улучшить? А то ведь я сварщик-то ненастоящий. Подскажете?