#devtools::install_github("sfsheath/cawd")
library(cawd)
library(sp)21 Пространственные данные в R
Пространственные данные и геомоделирование в R – большая тема, которую в этом курсе мы затронем очень кратко. Для более глубого знакомства можно рекомендовать на русском языке пособие Тимофея Самсонова “Визуализация и анализ географических данных на языке R”, а на английском созданную коллективом автором работу “Geocomputation with R”. Дальше я во многом опираюсь на эти работы.
21.1 Данные: римские амфитеатры
Данные для этого урока происходят из пакета cawd (Collected Ancient World Data), который, в свою очередь, опирается на следующие ресурсы:
- Digital Atlas of the Roman Empire;
- Ancient World Mapping Center;
- Геопространственная сетевая модель Римской империи Orbis.
Мы заберем из пакета датафрейм с римскими амфитеатрами (подробнее о нем можно прочитать здесь) и карту Римской империи на 200 г. н.э. (в формате sp, который представляет собой немного устаревший, но легко конвертируемый формат хранения пространственных данных в R).
class(awmc.roman.empire.200.sp)[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"Объект sp имеет свой метод plot().
par(mai=c(0,0,0,0))
plot(awmc.roman.empire.200.sp)
Для начала выберем нужные столбцы из датафрейма с данными об амфитеатрах.
library(tidyverse)
ramphs <- cawd::ramphs |>
dplyr::select(label, longitude, latitude, capacity, type, prov.type)
ramphs21.2 Simple Features
Современный формат хранения векторных геоданных называется Simple Features. Основное отличие объектов sf от объектов sp в том, что данные хранятся в виде датафрейма со списком-колонкой для хранения геометрии (линии, точки или полигона). Эта колонка называется sfc (simple features geometry column), а сама геометрия внутри нее – sfg (simple feature geometry).
То, что объекты типа Simple Features реализованы в виде самых обычных фреймов данных, означает, что любая операция, применимая к фрейму данных, будет также применима к объекту типа
sf. Это очень важная особенность объектов типаsf, которой сильно не хватало в экосистеме исторического пакетаsp. – Источник.
library(sf)Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUEroman_map <- sf::st_as_sf(awmc.roman.empire.200.sp)
roman_mapПосмотрим внимательно на это описание.
bounding box: прямоугольная рамка, которая задает границы карты; здесь координаты по оси x соответствуют долготе, а по оси y - широте. Будьте внимательны, потому что мы пишем обычно сначала широту, а потом долготу. Но по аналогии с алгеброй x определяет сдвиг вправо-влево (долготу), в то время как y - вверх-вниз (широту). Получается, что перед нами кусочек северного полушария, в основном к востоку от нулевого меридиана.
library(tmaptools)
bb(roman_map) xmin ymin xmax ymax
-9.48732 22.89549 43.10774 55.10117 Geodetic CRS: +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 – это определение геодезической системы координат (Geodetic CRS):
+proj=longlat- указывает, что используется географическая система координат (широта и долгота).+datum=WGS84- определяет геодезическую основу, в данном случае это Всемирная геодезическая система 1984 года (World Geodetic System 1984).+ellps=WGS84- указывает, что используется эллипсоид, соответствующий системе WGS84.+towgs84=0,0,0- определяет параметры трансформации между используемым эллипсоидом и эллипсоидом системы WGS84. Значения “0,0,0” означают, что никаких трансформаций не требуется, так как данные уже находятся в системе WGS84.
Помимо географической системы координат, которые используют сферическую или эллипсоидальную поверхность Земли, бывают проекционные (плоские) системы, которые используют плоскую (двумерную) поверхность. Кроме того, они используют другие единицы измерения: не градусы широты и долготы, а линейные единицы (например, метры).
Функция st_is_valid() проверяет, является ли заданная пространственная геометрия (например, точка, линия, многоугольник) топологически корректной. В нашем случае есть одна ошибка.
st_is_valid(roman_map) [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[25] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[37] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE
[49] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[73] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[97] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[109] TRUE TRUE TRUE TRUEНадо починить, иначе дальше будет ошибка.
roman_map <- st_make_valid(roman_map)21.3 Пакет tmap
Есть множество пакетов для работы с пространственными данными в R; мы начнем с одного из наиболее простого и интуитивно понятного tmap.
# install.packages("tmap")
# install.packages("tmaptools")
library(tmap)Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')tmap_mode("plot")tmap mode set to plottingtm_shape(roman_map) +
tm_fill(col = "magenta") +
tm_borders(col = "white") 
Пакет tmap предлагает хороший выбор стилей для оформления карты.
tmap_style("classic")tmap style set to "classic"other available styles are: "white", "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "watercolor" tm_shape(roman_map) +
tm_fill() +
tm_borders() 
Вручную можно добавить, например, компасс, координатную сетку и шкалу масштаба.
tm_shape(roman_map) +
tm_polygons() +
tm_graticules() +
tm_compass(type = "8star", position = c("right", "top")) +
tm_scale_bar(breaks = c(0, 100, 200), text.size = 1, position = c("right", "top")) Scale bar set for latitude km and will be different at the top and bottom of the map.
21.4 tm_bubbles()
Тепрь нанесем на карту отдельные амфитеатры в виде точек. Для этого преобразуем датафрейм в объект sf. Обратите внимание, что геометрии теперь другие.
amph_points <- ramphs |>
st_as_sf(coords = c("longitude", "latitude"))
amph_pointstmap_style("classic")
tm_shape(roman_map) +
tm_polygons() +
tm_compass(type = "8star",
position = c("right", "top")) +
tm_scale_bar(breaks = c(0, 100, 200),
text.size = 1,
position = c("left", "bottom")) +
tm_shape(amph_points) +
# тут наши точки; размер кодирует вместимость
tm_bubbles(size = "capacity",
alpha = 0.8,
scale = 1,
col = "prov.type",
palette = c("red", "blue", "green")) +
tm_layout(legend.position = c("right", "bottom"),
legend.frame = TRUE
)
Как и в ggplot, разные геометрии могут использовать разные данные. В нашем случае – отфильтрованный список названий.
tm_shape(roman_map) +
tm_polygons() +
tm_compass(type = "8star",
position = c("right", "top")) +
tm_scale_bar(breaks = c(0, 100, 200),
text.size = 1,
position = c("left", "bottom")) +
tm_shape(amph_points) +
tm_bubbles(size = "capacity",
alpha = 0.8,
scale = 1,
col = "prov.type",
palette = c("red", "blue", "green")
) +
# фильтр для названий
tm_shape(amph_points |> filter(capacity > 30000 )) +
# текст
tm_text("label") +
tm_layout(legend.position = c("right", "bottom"),
legend.frame = TRUE
)Warning: Currect projection of shape amph_points unknown. Long-lat (WGS84) is
assumed.Warning: Currect projection of shape filter(amph_points, capacity > 30000)
unknown. Long-lat (WGS84) is assumed.Scale bar set for latitude km and will be different at the top and bottom of the map.
21.5 Пересечения между геометриями
Мы можем посчитать, сколько точек приходится на один многогранник. Но для этого координатные системы должны совпадать. Сейчас у точек нет никакой CRS, в чем легко убедиться.
st_crs(amph_points)Coordinate Reference System: NAst_intersects(roman_map, amph_points)Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, : st_crs(x) == st_crs(y) is not TRUEЧтобы избавиться от ошибки, необходимо назначить или трансформировать координатные системы. Четыре цифры ниже представляют собой код EPSG (European Petroleum Survey Group). Это один из способов задания (хранения) пространственной привязки. EPSG:4326 соответствует WGS84, а Web Mercator – EPSG:3857.
# если нужно трансформировать
roman_map <- st_transform(roman_map, 4326)
# если нужно назначить
amph_points <- st_set_crs(amph_points, 4326)Снова уточним пересечения.
# пересечения
inter <- st_intersects(roman_map, amph_points)
interSparse geometry binary predicate list of length 112, where the
predicate was `intersects'
first 10 elements:
1: 45, 46, 69, 78, 80, 88, 108, 109, 110, 151, ...
2: (empty)
3: (empty)
4: (empty)
5: (empty)
6: (empty)
7: 3, 21, 38, 39, 47, 74, 89, 114, 115, 124, ...
8: (empty)
9: (empty)
10: (empty)Добавим новый столбец в датафрейм с картой.
roman_map$count <- lengths(inter)Теперь его можно использоать для выбора цвета заливки.
tm_shape(roman_map) +
tm_polygons(col = "count") +
tm_layout(legend.position = c("right", "bottom"),
legend.frame = TRUE)
21.6 tm_basemap()
Парящая в вакууме империя не очень радует глаз; в таком случае стоит добавить растровое изображение ландшафта. Пока это доступно только для динамической карты, поэтому переключаемся в режим “view” (ниже представлен скриншот).
Будьте внимательны, совмещая исторические карты с современными! Убедитесь, что вы не показываете походы Цезаря в современную Швейцарию, как это произошло, например, здесь.
Фоновые карты для Рима доступны на сайте https://dh.gu.se/dare/.
tmap_mode("view")tmap mode set to interactive viewingtm_basemap("https://dh.gu.se/tiles/imperium/{z}/{x}/{y}.png") +
tm_shape(amph_points) +
tm_bubbles(size = "capacity", alpha = 0.8,
col = "white", border.col = "tomato")Legend for symbol sizes not available in view mode.21.7 Leaflet
Удобный способ создания интерактивных карт предлагает также пакет leaflet. Вызовем фон:
library(leaflet)
leaflet() |>
addTiles() Ccылка на галерею подложек (подсмотрена в курсе Георгия Мороза).
library(leaflet)
leaflet() |>
addProviderTiles("Esri.WorldImagery") Некоторые подложки потребуют аутентификации Для этого надо зарегистрироваться на сайте https://stadiamaps.com/ (это бесплатно), создать в личном кабинете Property и прописать доменное имя для карты. Например, акварельная подложка при публикации в Сети требует аутентификации. Многие другие работают без нее (и почти все – локально).
# цветовая палитра
ramphs$type <- factor(ramphs$type)
factpal <- colorFactor(palette = c("#DE7424FF", "#F5CA37FF", "#AD8D26FF", "#496849FF", "#654783FF"),
ramphs$type)
ramphs |>
leaflet() |>
addProviderTiles("Stadia.StamenWatercolor") |>
addCircles(lng = ~longitude,
lat = ~latitude,
color = ~factpal(type),
opacity = 0.7,
popup = ~paste0(
label,
"</br>",
capacity)
) |>
addLegend(pal = factpal,
values = ~type)Заменим кружки на маркеры и сгруппируем их. Наложим это все на снимок из космоса (ок, это просто демо, с освещением у них было не очень).
ramphs |>
leaflet() |>
addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") |>
addMarkers(lng = ~longitude,
lat = ~latitude,
popup = ~paste0(
label,
"</br>",
capacity),
clusterOptions = markerClusterOptions()
) Заменим маркеры на изображения амфитеатров.
my_icon <- makeIcon(
iconUrl = "./images/amphitheatre.png",
iconWidth = 31*215/230,
iconHeight = 31,
iconAnchorY = 16,
iconAnchorX = 31*215/230/2
)
ramphs |>
leaflet() |>
addProviderTiles("Esri.WorldTerrain") |>
addMarkers(icon = ~my_icon,
clusterOptions = markerClusterOptions())Assuming "longitude" and "latitude" are longitude and latitude, respectively21.8 Возможности ggplot2
Для статичных карт можно использовать привычный ggplot(), как показано, например, здесь.
ggplot() +
geom_sf(data = roman_map) +
geom_point(ramphs,
mapping = aes(longitude, latitude),
color = "steelblue",
alpha = 0.5) +
theme_bw()
Если точек много, то может быть уместней представить на карте плотность их распределения.
ggplot() +
geom_sf(data = roman_map, fill = "wheat") +
geom_point(ramphs, color = "steelblue", alpha = 0.5,
mapping = aes(longitude, latitude)) +
geom_density2d(data = ramphs,
mapping = aes(longitude, latitude,
color = after_stat(level)),
linewidth = 1, alpha = 0.5)
Еще один способ отразить области скопления точек – сотовая диаграмма. На такой диаграмме координатная плоскость разбивается на гексагоны, которые закрашиваются в соответствии с градиентом плотности попавших в них точек.
g <- ggplot() +
geom_sf(data = roman_map, fill = "wheat") +
geom_hex(data = ramphs,
mapping = aes(longitude, latitude),
bins = 25,
color = "royalblue") +
theme_bw() +
scale_fill_continuous(trans = "reverse")
g
21.9 Пакет plotly
Пакет plotly позволяет добавить интерактивности на карту.
library(plotly)
ggplotly(g)21.10 Данные: римские дороги
С точками и многоугольниками все понятно; но на карту можно нанести, в виде линий, и сетевые данные. Например, данные проекта Orbis о транспортном сообщении в Римской империи. Датасет можно забрать по ссылке https://purl.stanford.edu/mn425tz9757 или загрузить из репозитория курса (ребра и вершины).
orbis_e <- read_csv("../files/orbis_edges_0514.csv")
orbis_eorbis_n <- read_csv("../files/orbis_nodes_0514.csv")
orbis_nlibrary(igraph)
orbis_g <- graph_from_data_frame(orbis_e,
vertices = orbis_n,
directed = FALSE)Без учета расстояний в пространстве сеть дорог выглядит так. Красная точка – Рим.
library(ggraph)
set.seed(25092024)
ggraph(orbis_g) +
geom_edge_link(aes(color = type)) +
geom_node_point(color = "grey30") +
geom_node_point(color = "tomato", size = 3,
aes(filter = V(orbis_g)$label == "Roma"))Warning: Existing variables `x` and `y` overwritten by layout variables
21.11 Сеть на карте
Кооринаты узлов – это широта и долгота (главное не перепутать).
orbis_coord <- orbis_e |>
left_join(orbis_n, by = join_by(source == id)) |>
mutate(source = label, .before = target) |>
select(-label) |>
# координаты начала ребра
rename(x1 = y, y1 = x) |>
left_join(orbis_n, by = join_by(target == id)) |>
mutate(target = label, .after = source) |>
select(-label) |>
# координаты конца ребра
rename(x2 = y, y2 = x) |>
# отрезаем пуповины к центру мира
filter(x1 != 0, y1 !=0, x2 != 0, y2 != 0)# для простоты пока берем современную карту
world <- map_data("world")
ggplot(data = world, aes(long, lat)) +
geom_map(map = world, aes(map_id = region),
fill = "wheat", color = "grey") +
geom_point(data = orbis_coord, aes(x = x1, y = y1),
color = "steelblue", alpha = 0.5) +
coord_map(xlim = c(-10, 50),
ylim = c(23, 54)) +
geom_segment(data = orbis_coord,
aes(x = x1, y = y1, xend = x2, yend = y2,
color = type))
Параллельные линии создают шум в нижней правой четверти; есть несколько способов от этого избавиться, но мы пока просто отрежем часть городов (восточнее Берениса).
orbis_coord_pruned <- orbis_coord |>
filter(y1 > 28 & y2 > 28)
library(paletteer)
cols <- paletteer_d("basetheme::brutal")par(mar = rep(0,4))
set.seed(24092024)
ggplot(data = world, aes(long, lat)) +
geom_map(map = world, aes(map_id = region),
fill = "white", color = "wheat") +
geom_point(data = orbis_coord, aes(x = x1, y = y1),
color = cols[1], alpha = 0.5) +
geom_segment(data = orbis_coord_pruned,
aes(x = x1, y = y1, xend = x2, yend = y2,
color = type)) +
geom_label(data = orbis_coord |>
filter(source %in% c("Roma", "Alexandria", "Carthago", "Sirmium", "Corinthus", "Antiochia", "Londinium", "Tarraco", "Augusta Taurinorum", "Jerusalem")),
aes(x1, y1, label = source),
color = cols[5],
label.size = 0.15,
fontface = "bold") +
coord_map(xlim = c(-10, 45),
ylim = c(26, 54)) +
labs(x = NULL, y = NULL,
title = "Транспортное сообщение в Римской империи",
subtitle = "Данные проекта Orbis") +
theme_bw(base_family = "serif") +
theme(legend.position="bottom",
legend.box = "horizontal",
panel.background = element_rect(fill = "aliceblue"),
text = element_text(color = cols[5])) +
scale_color_manual("тип", values = sample(cols, 10))