Коротко о главном: прочитал недавно пост 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. Что точно надо переделать/улучшить? А то ведь я сварщик-то ненастоящий. Подскажете?
Автор: sobach
