Часть 2
knitr::opts_chunk$set(echo = TRUE)
#в основном здесь нужен будет dplyr (для работы с дата фреймом) и ggplot для графики
library(tidyverse)
#пакет для чтения таблиц excel
library(readxl)
#пакет для того, чтобы загрузить данные с http://www.naturalearthdata.com/
#загрузка в формате shape файлов
#нужен будет, чтобы загрузить границы стран
library(rnaturalearth)
#пакет для работы с пространственными данными
#сейчас используется намного меньше, чем следующий пакет, но все же используется, потому что от него зависят некоторые другие пакеты
library(sp)
#основной и самый используемый сейчас пакет для работы с векторными пространственными данными
library(sf)
#пакет для создание тематических карт
library(cartography)
#сборник цветовых палитр
library(RColorBrewer)Первая часть: https://rpubs.com/baltti/thematic-maps-ggplot
Данные будут те же, что и в первой части.
Сначала подготовка данных. В файле, скачанном с http://happyplanetindex.org/ лучше поправить заголовки: часть заголовков идет в 2 абзаца в исходном документе, а также убрать пробелы в заголовке (из-за этого могут выдаваться ошибки)
countries<-ne_countries(scale = "medium", continent = "europe",
type = "countries", returnclass = "sf") %>%
st_crop(c(xmin = -10, xmax = 45,
ymin = 34, ymax = 71))%>%
st_transform(3857)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometrieshpi<- sp::merge(countries,
read_excel("hpi-data-2016.xlsx",col_names = TRUE,
sheet = 5,range = "B6:O146"),
by.x="sovereignt",by.y="Country")%>%
select(1,64:77)Тематические карты с помощью пакета cartography
Общая информация о пакете: http://riatelab.github.io/cartography/docs/index.html
Статья от авторов пакета, посвященная воспроизводимой картографии: https://riatelab.github.io/ReproducibleCartography/paper/paper.html
Также один из авторов ведет блог, где можно почитать про отдельные фукнции и обновления пакета https://rgeomatic.hypotheses.org/category/cartography
Логика создания тематических карт с помощью этого пакета очень простая: отрисовывается геометрия объектов, на которую накладывается нужный тематический слой, а также дополнительные элементы оформления (заголовок, легенда и прочее). Кроме того, можно добавлять подложку.
В пакете реализованы все основные виды тематических карт: картограмма, пропорциональные символы, грид, изоплеты, карта плотности точек, карты потоков. Совсем недавно появились функции добавления на карту waffleplot, облака слов, штриховки и еще целый ряд.
Пропорциональные символы или bubble map
Самый простой способ отображения числовых показателей.
Размер символа соответствует величине показателя, сами символы размещаются в геометрическом центре исходных полигонов.
Для создания такой карты используется функция propSymbolsLayer, аргументами которой являются объект, который нужно отобразить, - sf объект, объект SpatialPointsDataFrame или SpatialPolygonsDataFrame (это объекты, которые характерны для пакета sp); числовая характеристика, по которой задается величина символов, и параметры символов.
Здесь важно именно перед добавлением тематического слоя отрисовать геометрию объектов, потом уже добавлять все остальное. Добавим слой пропорциональных символов по умолчанию
par(mar = c(0, 0, 2, 0)) #отступы от края
#просто отрисовка полигонов
plot(hpi$geometry)
#добавление пропорциональных символов
propSymbolsLayer(x=hpi,
var = "Happy_Planet_Index") #переменная, на основе которой наносятся символы
Символы можно довольно хорошо настроить, а также добавить макет с заголовком, подзаголовком, автором, стрелкой на север, темой и прочее.
Палитра, использованная здесь для темы wine.pal - одна из палитр, встроенных в пакет. Со всеми палитрами можно ознакомиться командой display.carto.all()
par(mar = c(0, 0, 2, 0)) #отступы от края
#просто отрисовка полигонов
plot(hpi$geometry, col = "#45a16b", border = "black")
#добавление пропорциональных символов
propSymbolsLayer(hpi, inches = 0.3, #размер самого большого символа
lwd = 0.5, #ширина границы символа
col="#167b9c", #цвет символов
var = "Happy_Planet_Index", #переменная, на основе которой наносятся символы
symbols = "square", #тип символа
legend.pos = "topleft", #где расположена легенда
legend.values.rnd = 0, #количество знаков после запятой у чисел в легенде
legend.title.txt = "Индекс счастья", #заголовок легенды
legend.frame = F)
#слой макета
#задает расположение заголовка, подзаголовка, дополнительные символы
layoutLayer(title = "propSymbolsLayer()", #заголовок
#автор и источник данных либо дополнительны сведения
author = "https://t.me/geomess",
sources = "http://happyplanetindex.org/, http://www.naturalearthdata.com/",
frame = FALSE, #рамка
scale = NULL, #шкала масштаба
theme = "green.pal", #палитра цветов
north = TRUE) #направление на север
Отображение двух характеристик одновременно с помощью символов
Кроме простых пропорциональных символов (propSymbolsLayer), можно воспользоваться функцией propSymbolsChoroLayer, которая позволяет задать размер символов в зависимости от одной характеристики, а цвет в зависимости от другой. То есть получается по сути bivariate map, только закодированная не цветовой палитрой, а размером и цветом.
par(mar = c(0, 0, 2, 0)) #отступы от края
plot(hpi$geometry, col = "#8a96a8", border = "black")
propSymbolsChoroLayer(hpi,
var = "GDP", inches = 0.2, #параметр, по которому задается размер
var2 = "Happy_Planet_Index", #параметр, по которому задается цвет
col = carto.pal(pal1 = "orange.pal", n1 = 4),
symbols = "circle",
method = "equal", nclass = 4, #метод классификации и количество интервалов
border = "black", lwd = 0.5,
#параметры легенды для первой характеристики
legend.var.pos = "topleft",
legend.var.values.rnd = 0,
legend.var.title.txt = "ВВП/чел.",
legend.var.style = "e",
#параметры легенды для второй характеристики
legend.var2.pos = "bottomleft",
legend.var2.values.rnd = 0,
legend.var2.title.txt = "Индекс счастья")
# layout
layoutLayer(title = "propSymbolsChoroLayer()",
author = "https://t.me/geomess",
frame = FALSE, theme = "orange.pal",
scale = 500, north = FALSE)
Картограмма или choropleth
Картограмма она же choropleth map создается командой choroLayer
Здесь палитра задана через команду из пакета RColorbrewer, но можно пользоваться встроенными палитрами, а также создавать палитры вручную.
par(mar = c(0, 0, 2, 0)) #отступы от края
plot(hpi$geometry, col = "white", border = "black")
choroLayer(hpi, var = "Happy_Planet_Index",
col = brewer.pal(n = 4,name = "Pastel2"), #цветовая палитра
method = "equal", #метод классификации
nclass = 4, #количество групп
border = "black", lwd = 0.5, #цвет и ширина границ
legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
layoutLayer(title = "choroLayer()", author = "https://t.me/geomess",
sources = "http://happyplanetindex.org/, http://www.naturalearthdata.com/",
frame = FALSE, scale = NULL,
theme = "green.pal", north = FALSE)
На карту также можно добавить подписи. Они добавляются как отдельный слой Label Layer
par(mar = c(0, 0, 2, 0)) #отступы от края
plot(hpi$geometry, col = "white", border = "black")
choroLayer(hpi, var = "Happy_Planet_Index",
col = brewer.pal(n = 4,name = "Pastel2"), #цветовая палитра
method = "equal", #метод классификации
nclass = 4, #количество групп
border = "black", lwd = 0.5, #цвет и ширина границ
legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
labelLayer(hpi,txt="sovereignt", #колонка с названиями
col = "#4d4d4d",
overlap = FALSE, #допустимо ли наложение подписей
cex = 0.7) #размер подписей
layoutLayer(title = "Картограмма (choropleth map)", author = "https://t.me/geomess",
sources = "http://happyplanetindex.org/, http://www.naturalearthdata.com/",
frame = FALSE, scale = NULL,
theme = "green.pal", north = FALSE)
Одно из преимуществ этого пакета при отрисовке тематических карт над пакетом ggplot то, что не нужно отдельно рассчитывать классификацию параметра. При использовании пакета ggplot предварительно нужно рассчитать интервалы и создать для этого отдельный параметр, в пакете cartography нужно просто указать тип классификации в качестве одного из аргументов
par(mar = c(0, 0, 1.2, 0), mfrow = c(2, 2))
#разбивка на равные интервалы
plot(hpi$geometry, col = "white", border = "black")
choroLayer(hpi, var = "Happy_Planet_Index",
col = brewer.pal(n = 5,name = "Pastel2"), #цветовая палитра
method = "equal", #метод классификации
nclass = 5, #количество групп
border = "black", lwd = 0.5, #цвет и ширина границ
legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
layoutLayer(title = "Равные интервалы",
frame = FALSE, scale = NULL,
theme = "green.pal", north = FALSE)
#разбивка на квантили
plot(hpi$geometry, col = "white", border = "black")
choroLayer(hpi, var = "Happy_Planet_Index",
col = brewer.pal(n = 5,name = "Pastel2"), #цветовая палитра
method = "quantile", #метод классификации
nclass = 5, #количество групп
border = "black", lwd = 0.5, #цвет и ширина границ
legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
layoutLayer(title = "Квантили",
frame = FALSE, scale = NULL,
theme = "green.pal", north = FALSE)
#разбивка на основе стандартного отклонения
plot(hpi$geometry, col = "white", border = "black")
choroLayer(hpi, var = "Happy_Planet_Index",
col = brewer.pal(n = 5,name = "Pastel2"), #цветовая палитра
method = "msd",
nclass = 5, #количество групп
border = "black", lwd = 0.5, #цвет и ширина границ
legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
layoutLayer(title = "Стандартное отклонение",
frame = FALSE, scale = NULL,
theme = "green.pal", north = FALSE)
#разбивка на естественные интервалы
plot(hpi$geometry, col = "white", border = "black")
choroLayer(hpi, var = "Happy_Planet_Index",
col = brewer.pal(n = 5,name = "Pastel2"), #цветовая палитра
method = "fisher-jenks",
nclass = 5, #количество групп
border = "black", lwd = 0.5, #цвет и ширина границ
legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
layoutLayer(title = "Естественные интервалы",
frame = FALSE, scale = NULL,
theme = "green.pal", north = FALSE)
Кроме стандартных вариантов также можно сделать картограмму в виде “раскрашенных карандашами” полигонов. Для этого нужно предварительно рассчитать этот слой (getPencilLayer), а потом строить картограмму на его основе вместо исходного. Эта функция создает MULTILINESTRING объекты, заполняющие исходные полигоны.
par(mar = c(0, 0, 2, 0)) #отступы от края
plot(hpi$geometry, col = "white", border = "black")
hpi_pencil<-getPencilLayer(hpi,
size = 200, #частота штриховки
buffer = -2) #буфер вокруг границ
choroLayer(hpi_pencil, var = "Happy_Planet_Index",
col = brewer.pal(n = 4,name = "Spectral"), #цветовая палитра
method = "equal", #метод классификации
nclass = 4, #количество групп
border = "white", lwd = 0.5, #цвет и ширина границ
legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
layoutLayer(title = "getPencilLayer()", author = "https://t.me/geomess",
sources = "http://happyplanetindex.org/, http://www.naturalearthdata.com/",
frame = FALSE, scale = NULL,
theme = "green.pal", north = FALSE)
Штриховка
Совсем недавно появившаяся функция. Отличается от предыдущей тем, что можно выбрать из нескольких видов паттернов (в то числе задать паттерн на основе текста), а также тем, что штриховка является регулярной.
Подробнее про слой штриховки с примерами можно почитать здесь: https://dieghernan.github.io/cartographyvignette/
par(mar = c(0, 0, 2, 0)) #отступы от края
plot(hpi$geometry, col = "#90c9d5", border = "black")
hatchedLayer(hpi, pattern = "zigzag",
density = 5, #частота паттерна
add=TRUE) #нужно ли добавить слой к предыдущему
layoutLayer(title = "hatchedLayer()", author = "https://t.me/geomess",
frame = FALSE, scale = NULL,
theme = "grey.pal", north = FALSE)
Также в пакете есть возможность добавления подложки с открытых картографических серверов.
Добавление подложки нужно делать в 2 этапа: сначала загрузить ее (getTiles), а потом отрисовать (tilesLayer)
par(mar = c(0, 0, 2, 0)) #отступы от края
#загрузка подложки
tile<-getTiles(hpi,type="Stamen.Watercolor", crop=TRUE)
#отрисовка
tilesLayer(tile)
hatchedLayer(hpi, pattern = "diamond", density = 5,add=TRUE)
layoutLayer(title = "hatchedLayer()+tilesLayer()", author = "Stamen.Watercolor",
frame = FALSE, scale = NULL,
theme = "grey.pal", north = FALSE)
Так как отрисовка осуществляется слоями, то можно выбрать отдельные объекты и добавить штриховку только для них
post_communist<- hpi %>% filter(Region=="Post-communist")
#разбивка на естественные интервалы
plot(hpi$geometry, col = "white", border = "black")
choroLayer(hpi, var = "GDP",
col =rev (carto.pal(pal1 = "grey.pal",n1=5)), #цветовая палитра
method = "fisher-jenks",
nclass = 5, #количество групп
border = "black", lwd = 0.5, #цвет и ширина границ
legend.pos = "topleft",
legend.title.txt = "ВВП $/чел.", add = T)
hatchedLayer(post_communist, pattern = "left2right", density = 2,add=TRUE,
col="#f91f1f",lwd=2)
Карта плотности точек или dot density
Карта плотности точек задается функцией dotDensityLayer.
Для этого слоя обязателбно нужно указать, какой величине соответствует одна точка. Параметры типа точек (pch) соотвествуют параметрам базовой графической функции points.
par(mar = c(0, 0, 2, 0)) #отступы от края
plot(hpi$geometry)
dotDensityLayer(hpi, var="GDP",n=1000, #величина, которой равна 1 точка
cex = 0.6, type = "random", #размер и паттерн размещения точек
col = "#227745", pch = 19,
legend.txt = "ВВП 1000 $/чел.", legend.frame = FALSE)
layoutLayer(title = "dotDensityLayer()",
author = "https://t.me/geomess",
frame = FALSE, scale = NULL, theme = "green.pal", north = FALSE)
Комбинация тематических слоев
Так как тематические слои добавляются последовательно, то можно их комбинировать между собой.
Например, отрисовать пропорциональные символы поверх картограммы.
par(mar = c(0, 0, 2, 0)) #отступы от края
plot(hpi$geometry)
#нанесем на карту индекс счастья в виде картограммы
choroLayer(hpi, var = "Happy_Planet_Index",
col = carto.pal(pal1 = "sand.pal",n1 = 4), method = "fisher-jenks",
nclass = 4, border = "white", lwd = 0.5, legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
#и ВВП на душу населения пропорциональными символами
propSymbolsLayer(hpi, inches = 0.2, lwd = 1.25,
var = "GDP",
col = NA, border = "#6c2213", #символы будут отображаться только в виде контура
legend.pos = "bottomleft",
legend.values.rnd = 0,
legend.title.txt = "ВВП/чел.",
legend.frame = F)
layoutLayer(title = "propSymbolsLayer() + choroLayer()",
author = "https://t.me/geomess",
frame = FALSE, scale = NULL, theme = "sand.pal", north = FALSE)
Изоплеты или isopleth maps
На таких картах отображаются линии равного значения. Самый распространенный тип подобных карт - это карта рельефа с горизонталями.
Как правило, используются для отображения непрерывных величин.
Сама по себе карта строится функцией smoothLayer, которая осуществляет интерполяцию данных с учетом взаимодействия точек друг с другом.
Интерполяция здесь осуществляется не линейная, а на основе вычисдения потенциалов. Это сделано на базе потенциалов Стюарта.
Параметры интерполяции - span и beta. Первый показывает расстояние, на котором плотность вероятности функции достигает 0,5, то есть то расстояние, на котором точки перестают оказывать влияние друг на друга. Второй параметр коэффициент сопротивления для пространственной функции, то есть насколько среда проницаема для рассматриваемого параметра. Чем выше значение beta, тем выше сопротивление среды, то есть тем ближе горизонтали будут располагаться к исходной точке. И наоборот: чем меньше этот параметр, тем легче показатель “растекается” по территории.
par(mar = c(0, 0, 2, 0)) #отступы от края
smoothLayer(x = hpi, var = "GDP",
span = 100000, beta = 1,
mask = hpi, border = NA,
nclass = 9, #количество интервалов
col = rev(brewer.pal(n = 9,name = "YlGnBu")),
legend.title.txt = "ВВП $/чел.",
legend.pos = "topright", legend.values.rnd = 0)
layoutLayer(title = "smoothLayer()",
author = "https://t.me/geomess",
frame = FALSE, scale = NULL, theme = "blue.pal", north = FALSE)
Есть и другие варианты интерполяции и отображения изоплет, например, можно предварительно сделать грид, по которому рассчитать интерполированные значения, и отрисовать.
Грид
Для грида значения отображаются согласно ячейкам регулярной сетки (как правило, прямоугольник/квадрат или шестиугольник). Также как и изоплеты их рекомендуют использовать для непрерывных величин.
Фактически сама сетка со значениями здесь строится как слой картограммы (choropleth), но предварительно нужно построить эту сетку и рассчитать значения в каждой ячейке функцией getGridLayer.
grid <- getGridLayer(x = hpi,
cellsize = median(as.numeric(st_area(hpi)))/4, #размер ячейки
var = "Population",type = "regular")
par(mar = c(0, 0, 2, 0)) #отступы от края
#нанесем на карту страны
plot(hpi$geometry, col = "#bff2f0", border = "black")
# Plot the population density
choroLayer(x = grid, var = "Population", method = "fisher-jenks", nclass=9,
col = brewer.pal(n = 9,name = "PuBuGn"), border = "grey80",
lwd = 0.5, legend.pos = "topright", add = TRUE,
legend.title.txt = "Население",
legend.values.rnd =0)
layoutLayer(title = "Grid",
frame = FALSE, north = FALSE,
theme = "turquoise.pal")
Плиточные карты (hexbin grid map, 2d histogram map)
В плиточных картах объекты преобразуются в ячейки регулярной сетки, то есть один объект (в нашем случае - одна страна) - это одна ячейка. Так же как в в простом гриде основные типы ячеек - это прямоугольник/квадрат и шестиугольник. Такие карты рекомендуют использовать тогда, когда нужно отобразить абсолютный параметр так, чтобы люди не отвлекались на размер объектов.
Опять же как и в случае с гридом сама карта создается как обычная картограмма или choropleth, но предварительно объекты нужно преобразовать. Для этого нужно воспользоваться пакетом geogrid. Почитать про этот пакет можно у его автора вот здесь https://github.com/jbaileyh/geogrid
Для начала рекомендуется просмотреть возможные варианты комбинаций ячеек
#пакет для создания грида
library(geogrid)
#сначала лучше предварительно рассчитать, чтобы понять лучшую конфигурацию
#грид из шестиугольников
#сразу просчитываются возможные конфигурации
par(mfrow = c(2, 3), mar = c(0, 0, 2, 0))
for (i in 1:6) {
#функция преобразования полигонов в шестиугольники
new_cells <- calculate_grid(shape = hpi, grid_type = "hexagonal", seed = i)
plot(new_cells, main = paste("Seed", i, sep = " "))
}
#грид из квадратов
#сразу просчитываются возможные конфигурации
par(mfrow = c(2, 3), mar = c(0, 0, 2, 0))
for (i in 1:6) {
#преобразование полигонов в квадраты
new_cells <- calculate_grid(shape = hpi, grid_type = "regular", seed = i)
plot(new_cells, main = paste("Seed", i, sep = " "))
}
Собственно создание ячеек в соотвествии с выбранным вариантом, а также присоединение атрибутов из исходных объектов
#расчет грида
europe_hex <- calculate_grid(shape = hpi, grid_type = "hexagonal", seed = 1)
#присоединение атрибутов
resulthex <- assign_polygons(hpi, europe_hex)
#расчет грида
europe_reg <- calculate_grid(shape = hpi, grid_type = "regular", seed = 3)
#присоединение атрибутов
resultreg <- assign_polygons(hpi, europe_reg)Отрисуем три варианта для сравнения: обычная картограмма и две плиточные карты (палитра и классификация на всех картах одинаковые)
par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 3))
plot(hpi$geometry, col = "#62abe3", border = "black")
#нанесем на карту ВВП
choroLayer(hpi, var = "Happy_Planet_Index",
col = brewer.pal(n = 5,name = "Pastel2"), method = "fisher-jenks",
nclass = 5, border = "white", lwd = 0.5, legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
layoutLayer(title = "choroLayer()",
frame = FALSE, scale = NULL, theme = "sand.pal", north = FALSE)
plot(resulthex$geometry, col = "#62abe3", border = "black")
#нанесем на карту ВВП
choroLayer(resulthex, var = "Happy_Planet_Index",
col = brewer.pal(n = 5,name = "Pastel2"), method = "fisher-jenks",
nclass = 5, border = "white", lwd = 0.5, legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
labelLayer(resulthex,txt="admin", #колонка с названиями
col = "blue",
overlap = TRUE, #допустимо ли наложение подписей
cex = 0.4)
layoutLayer(title = "hexgrid map",
frame = FALSE, scale = NULL, theme = "sand.pal", north = FALSE)
plot(resultreg$geometry, col = "#62abe3", border = "black")
#нанесем на карту ВВП
choroLayer(resultreg, var = "Happy_Planet_Index",
col = brewer.pal(n = 5,name = "Pastel2"), method = "fisher-jenks",
nclass = 5, border = "white", lwd = 0.5, legend.pos = "topleft",
legend.title.txt = "Индекс счастья", add = T)
labelLayer(resultreg,txt="admin", #колонка с названиями
col = "blue",
overlap = TRUE, #допустимо ли наложение подписей
cex = 0.4)
layoutLayer(title = "2d histogram map",
frame = FALSE, scale = NULL, theme = "sand.pal", north = FALSE)