Shape-файлы

На этом занятии мы обсудим, как работать со статическими картами – картами, которые загружаются из файлов (без подключения к сервисам типа Яндекс-карты или Google Maps) и не поддерживают интерактив (нельзя уменьшать/увеличивать определенные области, двигать карту и прочее).

Один из широко распространенных форматов хранения географической информации – shape-файл. В этом файле хранятся географические (а точнее, геометрические) объекты разных типов: точки, линии, многоугольники и другие. Несмотря на то, что shape-файл является основным файлом, нужным для построения карты, одного его недостаточно: рядом с ним всегда есть необходимые вспомогательные файлы. В shp-файле хранится информация о геометрических объектах, в dbf – об их атрибутах, а shx-файл случит для связи между первыми двумя. Наверное, эта информация выглядит сейчас немного туманно, но главное следующее:

Чтобы понять, как выглядят нужные для построения карт файлы, скачаем набор файлов для России с сайта Global Administrative Areas. (Еще ресурсы, где можно найти файлы для карт России бесплатно: GIS-Lab, GisGeo). Так как уровень детализации может быть разный (вся страна, регионы или районы), файлы в папке тоже разные (_adm0, _adm1, _adm2, _adm3). Пока нам нужны только регионы (границы субъектов) – файл RUS_adm1.shp.

Построение карт (и их раскраска)

Для работы с картами нам понадобятся три библиотеки: maptools, RColorBrewer, classInt. Строго говоря, для самой карты нам нужен пакет maptools, а два остальных нужны для разбивки значений и выбора палитры цветов, в которые будет окрашиваться карта.

Установим их:

install.packages(c("maptools", "RColorBrewer", "classInt"))

Обратимся к ним:

library(maptools)
library(RColorBrewer)
library(classInt)

А теперь загрузим shape-файл (файл с расширением .shp):

# Для удобства сделаем папку со всеми географическими файлами рабочей.
setwd("/home/oem/Рабочий стол/RUS_adm_shp/")

rus <- readShapePoly("RUS_adm1.shp")

На загруженный файл можно посмотреть – как на обычную таблицу:

View(rus)

А теперь самое главное – построим карту.

plot(rus)

Похоже на Россию, правда? Только с Чукоткой беда: она не поместилась полностью, и поэтому одна ее часть осталась на востоке, а другая переместилась на запад. Это не проблема конкретного shape-файла, он не “кривой”, такая проблема будет возникать и с другими картами России, потому что Россия вытянута по широте. Но с этой трудностью можно справиться – вернемся к этому на следующей лекции, а пока сделаем что-то менее замысловатое.

Загрузим файл, который содержит детализацию по районам (муниципальным образованиям).

ru_map <- readShapePoly("RUS_adm2.shp")
View(ru_map)

Выберем Калужскую область и будем работать только с ней. С shp-файлом можно обращаться так же, как и с базой данных (выбирать нужные строки, добавлять столбцы и прочее), единственное, функции из библиотеки dplyr с ним корректно работать не будут. Поэтому для выбора области воспользуемся обычным subset():

kaluga_map <- subset(ru_map, NAME_1 == "Kaluga")
View(kaluga_map)

А теперь построим карту:

plot(kaluga_map)

Карта симпатичная, но строили мы ее для того, чтобы визуализировать какие-то данные с географической привязкой! Возьмем базу, содержащую демографические данные по Калужской области (данные Росстата с прошлого семинара):

info <- read.csv("kaluga.csv")
View(info)

На основе shp-файла создадим базу данных (просто сконвертируем):

kaluga_df <- as.data.frame(kaluga_map)
View(kaluga_df)

Теперь нужно как-то объединить карту и данные. Для того, чтобы это стало возможным, нужно, чтобы районы в обеих базах имели одинаковые названия. Кроме того, так как склеивать “географический” и “демографический” файл мы будем по столбцу с названиями районов, нужно, чтобы этот столбец тоже был назван одинаково в двух базах данных. Как можно заметить, с названиями районов на кириллице в shp-файле все не очень – они почему-то не в едином формате. Можно произвести замену, но в данном случае проще будет создать столбец ID_2 для районов в базе info:

info$ID_2 <- c(575, 576, 577, 579, 578, 601, 602, 585, 581, 586, 587, 589, 590, 
               591, 592, 593, 595, 596, 597, 598, 599, 580, 583, 600, 582, 594)

В наших демографических данных (info) нет деления на крупные города и районы (кроме Калуги и Обнинска), в частности, нет отдельных строк для городов Кирова и Людинова. Так как файл с географическими объектами для нас главный, его изменять не будем, а просто добавим две лишние строки в базу info. Раз Киров неотличим от Кировского района, а Людиново от Людиновского района, то строки можно просто скопировать:

kirov <- subset(info, district == info$district[8]) # Кировский
ludinovo <- subset(info, district == info$district[12]) # Людиновский
info <- rbind(info, kirov, ludinovo)

Но при этом изменить id:

info[27, "ID_2"] <- 584
info[28, "ID_2"] <- 588

Склеим две базы данных: “географическую” и “демографическую”. Понятно, что склеивать их имеет смысл по общему столбцу, в данном случае по столбцу с id районов. Склеим с помощью функции merge():

full <- merge(info, kaluga_df, by = "ID_2")
View(full)

Предположим, мы хотим нанести на карту данные по доле трудоспособного населения в районах: чем темнее цвет района на карте, тем больше доля трудоспособного населения в нем. Вопрос: каким образом мы будем разбивать значения доли, чтобы получать районы разных оттенков? Предложение такое: разбивать по процентилям. Упорядочить все значения доли трудоспособного населения по возрастанию и объединить их в группы: те, что входят в первые 20% – это первая группа (районы закрашиваются самым светлым цветом), во вторые 20% – это вторая группа (районы закрашиваются цветом потемнее), и так далее.

Разобьем значения доли трудоспособного населения на 5 групп.

brks <-classIntervals(full$wa_perc, n = 5, style = "quantile")

Посмотрим:

brks
## style: quantile
##   one of 12,650 possible partitions of this variable into 5 classes
## [52.43714,53.27346) [53.27346,54.29038) [54.29038,55.76317) 
##                   6                   5                   6 
## [55.76317,56.49207) [56.49207,59.39125] 
##                   5                   6
str(brks)
## List of 2
##  $ var : num [1:28] 56.3 54.5 59.2 53.5 56.5 ...
##  $ brks: num [1:6] 52.4 53.3 54.3 55.8 56.5 ...
##  - attr(*, "style")= chr "quantile"
##  - attr(*, "nobs")= int 26
##  - attr(*, "call")= language classIntervals(var = full$wa_perc, n = 5, style = "quantile")
##  - attr(*, "intervalClosure")= chr "left"
##  - attr(*, "class")= chr "classIntervals"

Выберем сами точки (границы) разбиения – brks:

brks <- brks$brks
length(brks) # их 6, то есть на 1 больше, чем самих групп
## [1] 6

А теперь выберем палитру цветов (из библиотеки RColorBrewer). В нашем случае выбран набор из 5 оттенков красного (та же история про разбивку на 5 групп):

colors <- brewer.pal(5, "Reds") 
length(colors)
## [1] 5

Какие палитры еще бывают, можно узнать, вызвав help:

help(brewer.pal)

Наконец, построим саму карту!

kaluga_map$wa_perc <- full$wa_perc
plot(kaluga_map, col = colors[findInterval(full$wa_perc, brks, all.inside = TRUE)], 
     axes = FALSE)

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

plot(kaluga_map, col = colors[findInterval(full$wa_perc, brks, all.inside = TRUE)], 
     axes = FALSE)
text(coordinates(kaluga_map), labels = kaluga_map$ID_2, col = "blue") # подписи районов
legend(x ='bottomleft', # положение легенды
       legend = leglabs(round(brks, 2)),  # подписи - округленные границы интервалов разбивки 
       fill = colors, # чем заполнены клетки в легенде
       title = 'Working people (%)',
       title.col = 'blue')

Обсудим, как можно сохранить карту в файл (вообще это касается любых графиков, не только карт). Понятно, что можно просто сделать это вручную: во вкладке Plots выбрать Export и экспортировать картинку в файл нужного формата. Но иногда этот способ не подходит: бывает, что картинка при ручном экспорте деформируется (особенно, если она большая). Поэтому посмотрим, как сохранять картинку с помощью кода. Сохраним в pdf:

# сначала обязательно строка с названием файла
pdf('map_Kaluga.pdf', width = 15, height = 10)

plot(kaluga_map, col = colors[findInterval(full$wa_perc, brks, all.inside = TRUE)], 
     axes=FALSE)

text(coordinates(kaluga_map), labels = kaluga_map$ID_2, col = "blue") # подписи районов

legend(x ='bottomleft', # положение легенды
       legend = leglabs(round(brks, 2)),  # подписи - округленные границы интервалов разбивки 
       fill = colors, # чем заполнены клетки в легенде
       title = 'Working people (%)')

# в конце обязательно "выйдем" из графика, иначе он не сохранится
# без следующей строки созданные pdf-файл будет пустым (даже без страниц)
dev.off()
## png 
##   2

Теперь вроде бы все хорошо. Но не очень: номера районов плохо видны на темно-красном цвете. Можно, конечно, просто поиграть с цветами и подобрать какой-то более удачный цвет, но можно сделать все более концептуально: выделять текст на темном фоне светлым цветом, а на более бледном – темным. Предлагаю подумать об этом сейчас немного и написать код, который бы создавал вектор цветов, в котором

Подсказка: загляните в вектор findInterval(full$wa_perc, brks, all.inside=TRUE)

# Ваш код, маэстро
# решение 

# номер цветовой группы, чем больше число, тем больше значения -> тем темнее цвет
coldrops <- findInterval(full$wa_perc, brks, all.inside = TRUE) 
# самые темные районы - с большими значениями доли трудоспособного населения
dark <- c(4, 5) 
# если цвет района темный, то цвет текста будет белый, иначе - черный
coltext <- ifelse(coldrops %in% dark, "white", "black")

Попробуем построить карту с новыми красивыми подписями:

plot(kaluga_map, col = colors[findInterval(full$wa_perc, brks, all.inside = TRUE)], 
     axes=FALSE)

text(coordinates(kaluga_map), labels = kaluga_map$ID_2, col = coltext) # подписи районов

legend(x ='bottomleft', 
       legend = leglabs(round(brks, 2)),  
       fill = colors, 
       title = 'Working people (%)')

Получилось!

Бонус: для тех, кому вдруг придется добавлять карты (и вообще картинки) в Rmd-файл: обратите внимание на опции fig.width и fig.height в ячейках с кодом {r}. С их помощью можно настроить размеры картинки в готовом html-файле – график будет иметь приличный вид, а не деформированный (“ужатый”), как бывает при просмотре в окне Plots (Zoom).