Часть 1

Введение

Про работу с пространственными данными в R есть подробные и хорошие книги

на английском Geocomputation with R https://geocompr.robinlovelace.net/

и на русском Визуализация и анализ географических данных на языке R https://tsamsonov.github.io/r-geo-course/

Для создания карт и визуализации пространственных данных может использоваться распространенный пакет ggplot с некоторыми пакетами дополнениями - виньетками. Но для тематического картографирования существует два широко используемых пакета cartography и tmap.

knitr::opts_chunk$set(echo = TRUE)

#сборник распространенных пакетов работы с данными
#в основном здесь нужен будет dplyr (для работы с дата фреймом) и ggplot для графики
library(tidyverse)
#графический пакет для добавления элементов и комбинирования графиков в ggplot
library(cowplot)
#пакет-виньетка для работы с пространственными данными в ggplot
library(ggspatial)

#пакет для чтения таблиц excel
library(readxl)

#пакет для того, чтобы загрузить данные с http://www.naturalearthdata.com/
#загрузка в формате shape файлов
#нужен будет, чтобы загрузить границы стран
library(rnaturalearth)

#пакет для работы с пространственными данными
#сейчас используется намного меньше, чем следующий пакет, но все же используется, потому что от него зависят некоторые другие пакеты
library(sp)

#основной и самый используемый сейчас пакет для работы с векторными пространственными данными
library(sf)

#сборник цветовых палитр
library(RColorBrewer)

Подготовка данных

Основным пакетом для работы с пространственными векторными данными сейчас является пакет sf. Он преобразует простраственные данные в дата-фрейм с колонкой geometry - объекты simple features (прочесть про них подробнее можно здесь https://github.com/r-spatial/sf). Такой подход и хранение данных позволяют использовать для анализа как методы специфичные для пространственных данных (например, поиск пересечений между объектами, расчет центроидов и прочее), так и методы работы с дата-фреймами, принятые в tidyverse.

До него использовался пакет sp, который и сейчас используется, так как некоторые пакеты все еще опираются на его работу.

Для создания объектов simple features (sf) можно использовать векторные данные, с которыми принято работать в геоинформационных системах, например, shape-файлы, geojson. Чтение таких файлов и создание из них объектов simple features осуществляется командой st_read

Здесь я воспользовалась данными с ресурса http://www.naturalearthdata.com/. Данные оттуда можно загружать в виде shape-файлов, но есть специальный пакет для работы с ними rnaturalearth. Этот пакет дает возможность загружать контура полигонов напрямую.

countries<-ne_countries(scale = "medium", continent = "europe",
                        type = "countries", returnclass = "sf")

plot(countries$geometry)

Перепроецирование объектов

Очевидно, что нужно дополнительно поработать с контурами, чтобы отобразить их более понятно Попробуем перепроецировать контура. Для этого используется команда st_transform. В качестве аргумента этой функции нужно указать целевую проекцию. Проекцию можно указать в виде номера по номенклатуре ESPG (ознакомиться с ней и подобрать подходящую проекцию можно здесь https://epsg.io/)

countries_reproj <- countries %>% 
  st_transform(3035) #проекция, предназначенная для континентальной Европы

plot(countries_reproj$geometry)

countries_mercator <- countries %>% 
  st_transform(3857) #проекция трансформируется в стандартную проекцию google maps (псевдо-Меркатор)

plot(countries_mercator$geometry)

Для того, чтобы карта была более читаемая немного обрежем контура стран. Так как здесь нет задачи целиком отобразить Россию, оставим только европейскую часть.

Для обрезки нужно указать координаты, ограничивающие нужный фрагмент. При обрезке здесь выходит предупреждение об атрибутивных данных: по умолчанию считается, что атрибуты одинаковы в любой точке объекта.

Также перепроецируем объекты, чтобы карта выглядела более привычно.

Если бы нужно было отобразать целиком всю территорию Европы и России, то в качестве проекции стоило бы выбрать либо проекцию для континентальной Европы, либо проекцию для территории России. Здесь я воспользуюсь проекцией псевдо-Меркатора или сферической проекцией Меркатора, которую используют для картографических сервисов Google maps, OSM и прочие.

countries <- st_crop(countries, c(xmin = -10, xmax = 45, 
                                  ymin = 34, ymax = 71)) %>% 
  st_transform(3857)

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

plot(countries$geometry)

Тематические карты будут составляться на основе индекса счастья и связанных с ним показателей.

Эти данные, а также методика расчета показаны http://happyplanetindex.org/ Данные в формате xlsx нужно предварительно скачать в свою рабочую папку отсюда http://happyplanetindex.org/resources

Полученные характеристики стран будут присоединены к объектам simple features с помощью функции пространственного присоединения

#чтение данных из файла
hpi_data<-read_excel("hpi-data-2016.xlsx",col_names = TRUE,
                     sheet = 5,range = "B6:O146")

hpi<- sp::merge(countries, #объект simple features
                hpi_data, 
        by.x="admin",by.y="Country")%>% #названия колонок, по которым происходит объединение
  select(1,64:77) #так как колонок получается очень много, здесь происходит выборка характеристик, которые будут использоваться дальше

Создание тематических карт с помощью пакета ggplot

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

Чтобы отрисовать sf-объект в самом пакете ggplot есть встроенная функция geom_sf(). Эта функция отсутствовала в более ранних версиях пакета.

Для векторных пространственных данных дополнительные пакеты - это прежде всего пакеты ggspatial и ggmap.

Первый позволяет добавлять пространственные векторные и растровые данные, а также направление на север, линейку с масштабом, аннотации, подписи.

Пакет ggmap предназначен для загрузки подложки с Google Maps и Stamen Maps, но этот пакет требует обязательного получения ключа для работы с API Google maps

Отрисуем контура стран

hpi %>%
  ggplot()+
  geom_sf()

Зададим сплошной цвет заливки и сменим тему оформления карты

hpi %>%
  ggplot()+ 
  geom_sf(fill="white")+ 
  theme_minimal()

Можно добавить подложку

hpi %>%
  ggplot()+ 
  annotation_map_tile(type = "stamenwatercolor") +
  geom_sf(alpha=0.5)+
  theme_minimal()

## Zoom: 2

Поменяем заливку на градиентную в зависимости от значения параметра Индекса счастья В результате получена карта choropleth или картограмма, в которой цветовая заливка полигонов зависит от значения числовой характеристики.

hpi %>%
  ggplot()+ 
  #добавление объектов с градиентной заливкой в зависимости от параметра
  geom_sf(aes(fill=`Happy Planet Index`))+
  theme_minimal()

Картограмма (choropleth map)

Очень хороший материал по созданию тематических карт только с использованием ggplot https://timogrossenbacher.ch/2016/12/beautiful-thematic-maps-with-ggplot2-only/ На этом материале основан код для следующих двух карт

Создадим картограмму с классификацией стран по величине индекса счастья на квантили. При классификации на основе квантилей в каждой группе будет содержаться одно и то же число объектов.

#разобьем все страны на несколько групп
no_classes <- 4 #количество групп 
labels <- c() #это будет вектор со значениями квантилей для легенды

#разбивка по квантилям
quantiles <- quantile(hpi$'Happy Planet Index', 
                      probs = seq(0, 1, length.out = no_classes + 1))

#округление границ 
labels <- c() 
for(idx in 1:length(quantiles)){
  labels <- c(labels, paste0(round(quantiles[idx], 0), 
                             " – ", 
                             round(quantiles[idx + 1], 0)))}
#удаление последнего элемента 
#потому что он будет подобным "66.62 - NA"
labels <- labels[1:length(labels)-1]
#дополнительная переменная с номером квантиля
hpi$hpi_quant <- cut(hpi$'Happy Planet Index', 
                                     breaks = quantiles, 
                                     labels = labels, 
                                     include.lowest = T)

#отрисуем карту с новой палитрой и разбивкой стран на группы 
hpi %>%
  ggplot() + #создание объекта ggplot
  #добавление объектов с цветовой заливкой в зависимости от параметра
  geom_sf(aes(fill=hpi_quant))+ 
  #палитра
  scale_color_brewer(type = "div",palette="RdYlGn",  
                     aesthetics = "fill",
                     name = "Индекс счастья", #заголовок легенды
                     #параметры легенды
                     guide = guide_legend(keyheight = unit(5, units = "mm"),
                       title.position = 'top', #расположение заголовка
                       reverse = T))+ #обратить последовательность 
  theme_minimal()+
  labs(caption = "source: http://happyplanetindex.org/")

Создадим картограмму с классификацией стран по величине индекса счастья на равные интервалы. При классификации на основе квантилей в каждой группе будет содержаться одно и то же число объектов.

#разбивка на равные интервалы
equal.interval <- seq(min(hpi$'Happy Planet Index'), 
                     max(hpi$'Happy Planet Index'), 
                     by = (max(hpi$'Happy Planet Index')-min(hpi$'Happy Planet Index'))/4)
labels <- c() 
for(idx in 1:length(equal.interval)){
  labels <- c(labels, paste0(round(equal.interval[idx], 1), 
                             " – ", 
                             round(equal.interval[idx + 1], 1)))}
#удаление последнего элемента 
#потому что он будет подобным "66.62 - NA"
labels <- labels[1:length(labels)-1]
hpi$hpi.equal <- cut(hpi$'Happy Planet Index', breaks=equal.interval, 
                     labels=labels,include.lowest = TRUE)
hpi %>%
  ggplot() +
  geom_sf(aes(fill=hpi.equal))+
  scale_color_brewer(type = "div",palette="RdYlGn",
                     aesthetics = "fill",
                     name = "Индекс счастья",
                     guide = guide_legend(
                       keyheight = unit(5, units = "mm"),
                       title.position = 'top',
                       reverse = T))+
  theme_bw()+
  #удаление сетки, подписей координат по осям
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(),axis.text =element_blank(),
        axis.ticks=element_blank())+
  labs(caption = "source: http://happyplanetindex.org/")

Bivariate map

Bivariate map - карта с одновременным отображением двух признаков. На этой карте признаки будут закодированы цветом.

За основу взят это материал https://timogrossenbacher.ch/2019/04/bivariate-maps-with-ggplot2-and-sf/ В качестве признаков будут использованы Индекс счастья и ожидаемая продолжительность жизни

Цветовую схему можно подобрать здесь https://observablehq.com/@benjaminadk/bivariate-choropleth-color-generator

#создаем 3 класса для продолжительности жизни
quantiles_ale <- hpi %>%
  pull(`Inequality-adjusted Life Expectancy`) %>%
  quantile(probs = seq(0, 1, length.out = 4))
# и 3 класса для ожидаемой продолжительности жизни
quantiles_hpi <- hpi %>%
  pull('Happy Planet Index') %>%
  quantile(probs = seq(0, 1, length.out = 4))

# нужно создать цветовую схему для кодировки переменных
# желтый цвет - ожидаемая продолжительность
# синий - индекс счастья
#каждой комбинации будет соответствовать свой цвет 
#цвета заданы в кодировке HEX - шестнадцатеричное представление RGB
bivariate_color_scale <- tibble(
  "3 - 3" = "#0d8000", # высокая продолжительность, высокий индекс счастья
  "2 - 3" = "#81b700",
  "1 - 3" = "#e8ca00", # низкая продолжительность, высокий индекс счастья 
  "3 - 2" = "#0d808b",
  "2 - 2" = "#81b78b", # средняя продолжительность, средний индекс счастья
  "1 - 2" = "#e8dc8b",
  "3 - 1" = "#0d80d9", # высокая продолжительность, низкий индекс счастья
  "2 - 1" = "#81b7e1",
  "1 - 1" = "#e8e8e8" # низкая продолжительность, низкий индекс счастья
) %>%
  gather("group", "fill")
# делим объекты на группы и присоединяем отдельным атрибутом 
hpi %<>%
  mutate(ale_quantiles = cut(`Inequality-adjusted Life Expectancy`,
      breaks = quantiles_ale,
      include.lowest = TRUE),
    hpi_quantiles = cut(`Happy Planet Index`,
      breaks = quantiles_hpi,
      include.lowest = TRUE),
    # вносим номера групп по обоим переменным в один столбец
    # это нужно, чтобы присоединить цветовую кодировку
    group = paste(as.numeric(ale_quantiles), "-",
      as.numeric(hpi_quantiles))
    ) %>%
  # теперь присоединяем цветовую кодировку
  # теперь у каждой страны есть свой цветовой код на основе значений
  #ВВП и ожидаемой продолжительности жизни
  left_join(bivariate_color_scale, by = "group")

#создание карты
map<-hpi%>%
  ggplot() +
  geom_sf(aes(fill=fill))+ 
  scale_fill_identity() +
  #подписи осей и заголовок
  labs(x = NULL,
       y = NULL,
       title = "Индекс счастья и ожидаемая продолжительность жизни в Европе")+
  theme_minimal()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(),axis.text =element_blank())+
  labs(caption = "source: http://happyplanetindex.org/")


#для создания легенды нужно выделить отдельно продолжительность и индекс счастья
#так как при подготовке цветовой схемы каждый цвет был соотнесен с группой
#группы закодированы в виде "1-2" 
bivariate_color_scale %<>%
  separate(group, into = c("Life_Expectancy", "Happy_Planet_Index"), sep = " - ") %>%
  mutate(Life_Expectancy = as.integer(Life_Expectancy),
         Happy_Planet_Index = as.integer(Happy_Planet_Index))

#подготовка легенды
legend <- ggplot() +
  geom_tile(data = bivariate_color_scale,
    mapping = aes(
      x = Life_Expectancy,
      y = Happy_Planet_Index,
      fill = fill)) +
  scale_fill_identity() +
  labs(x = "Увеличение продолжительности жизни-->",
       y = "Повышение уровня счастья-->") +
  theme_minimal() +
  # размер шрифта
  theme(axis.title = element_text(size = 6)) +
  # квадратные ячейки
  coord_fixed()
  
ggdraw() +
  draw_plot(map, 0, 0, 1, 1) +
  draw_plot(legend, 0.05, 0.075, 0.25, 0.25)

Анаморфная карта (cartogram)

Создадим анаморфную карту или cartogram. На такой карте отображается не реальный размер объектов, а искаженный в зависимости от величины показателя, то есть чем больше значение параметра, тем больше будет объект на карте.

Для этого понадобится одноименный пакет cartogram

library(cartogram)

#функция рассчитывает границы на основе способо "резинового листа"
#в качестве исходного берется sf объект или полигоны и параметр, на основе которого вводится искажение
#также можно задать условие для остановки: количество итераций (как здесь) или среднее значение ошибки
eu_cartogram<-cartogram_cont(hpi,'Happy Planet Index',
                             itermax = 8)

## Mean size error for iteration 1: 6.01643345309566

## Mean size error for iteration 2: 3.75315155946759

## Mean size error for iteration 3: 2.55633352194125

## Mean size error for iteration 4: 1.91697731091375

## Mean size error for iteration 5: 1.65372303707895

## Mean size error for iteration 6: 1.44012882332117

## Mean size error for iteration 7: 1.30789859135038

## Mean size error for iteration 8: 1.22729221298555

eu_cartogram %>%
  ggplot() +
  annotation_map_tile(type = "cartolight") + #добавим подложку, чтобы было понятно искажение
  geom_sf(aes(fill=hpi.equal,alpha=0.2))+
  scale_color_brewer(type = "div",palette="RdYlGn",
                     aesthetics = "fill",
                     name = "Индекс счастья",
                     guide = guide_legend(
                       keyheight = unit(5, units = "mm"),
                       title.position = 'top',
                       reverse = T))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(),axis.text =element_blank(),
        legend.position = "none",axis.ticks=element_blank())

## Zoom: 2