21  Пространственные данные в R

Пространственные данные и геомоделирование в R – большая тема, которую в этом курсе мы затронем очень кратко. Для более глубого знакомства можно рекомендовать на русском языке пособие Тимофея Самсонова “Визуализация и анализ географических данных на языке R”, а на английском созданную коллективом авторов работу “Geocomputation with R”. Далее я во многом опираюсь на эти руководства.

21.1 Данные: римские амфитеатры

Данные для этого урока происходят из пакета cawd (Collected Ancient World Data), который, в свою очередь, опирается на следующие ресурсы:

Мы заберем из пакета датафрейм с римскими амфитеатрами (подробнее о нем можно прочитать здесь) и карту Римской империи на 200 г. н.э. (в формате sp, который представляет собой немного устаревший, но легко конвертируемый формат хранения пространственных данных в R).

#devtools::install_github("sfsheath/cawd")
library(cawd)
library(sp)
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)

ramphs

21.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.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
roman_map <- sf::st_as_sf(awmc.roman.empire.200.sp)
roman_map
# Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
# Simple feature collection with 112 features and 8 fields
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: -9.48732 ymin: 22.89549 xmax: 43.10774 ymax: 55.10117
# Geodetic CRS:  +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# First 10 features:
#   OBJECTID         AREA  PERIMETER NEWDIO_ NEWDIO_ID ID Shape_Leng
# 0        1 19.612702708 35.3870861       2         0  0 35.1149460
# 1        2  0.080670307  1.2122280       3         0  0  1.2122280

Посмотрим внимательно на это описание.

  • 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):

  1. +proj=longlat - указывает, что используется географическая система координат (широта и долгота).
  2. +datum=WGS84 - определяет геодезическую основу, в данном случае это Всемирная геодезическая система 1984 года (World Geodetic System 1984).
  3. +ellps=WGS84 - указывает, что используется эллипсоид, соответствующий системе WGS84.
  4. +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)
tmap_mode("plot")
ℹ tmap mode set to "plot".
tmap_style("white") # default
style set to "white" (tmap default)
other available styles are: "gray", "natural", "cobalt", "albatross", "beaver", "bw", "classic", "watercolor"
tmap v3 styles: "v3" (tmap v3 default), "gray_v3", "natural_v3", "cobalt_v3", "albatross_v3", "beaver_v3", "bw_v3", "classic_v3", "watercolor_v3"
tm_shape(roman_map) +
  tm_fill(fill = "magenta") +
  tm_borders(col = "white") 

Пакет tmap предлагает хороший выбор стилей для оформления карты.

tmap_style("classic")
style set to "classic"
other available styles are: "white" (tmap default), "gray", "natural", "cobalt", "albatross", "beaver", "bw", "watercolor"
tmap v3 styles: "v3" (tmap v3 default), "gray_v3", "natural_v3", "cobalt_v3", "albatross_v3", "beaver_v3", "bw_v3", "classic_v3", "watercolor_v3"
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_scalebar(
    breaks = c(0, 500, 1000, 1500),
    text.size = 1, 
    position = c("left", "bottom")) 
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_points
tmap_style("classic")
style set to "classic"
other available styles are: "white" (tmap default), "gray", "natural", "cobalt", "albatross", "beaver", "bw", "watercolor"
tmap v3 styles: "v3" (tmap v3 default), "gray_v3", "natural_v3", "cobalt_v3", "albatross_v3", "beaver_v3", "bw_v3", "classic_v3", "watercolor_v3"
tm_shape(roman_map) +
  tm_polygons() +
  tm_compass(
    type = "8star",
    position = c("right", "top")
  ) +
  tm_scalebar(position = c("left", "bottom")) +
  tm_shape(amph_points) +
  tm_bubbles(
    size = "capacity", 
    size.scale = tm_scale_continuous(values.scale = 1),  
    fill = "prov.type", 
    fill.scale = tm_scale(values = c("red", "blue", "green")),
    fill_alpha = 0.8  
  ) +
  tm_layout(
    legend.position = c("right", "bottom"),
    legend.frame = TRUE
  )
Scale bar set for latitude km and will be different at the top and bottom of the map.

Как и в ggplot, разные геометрии могут использовать разные данные. В нашем случае – отфильтрованный список названий.

tmap_mode("plot")  # Убедимся, что мы в режиме рисования
ℹ tmap mode set to "plot".
tm_shape(roman_map) +
  tm_polygons() +
  tm_compass(type = "8star", position = c("right", "top")) +
  tm_scalebar(breaks = c(0, 100, 200), text.size = 1, position = c("left", "bottom")) +

  # Слой с пузырями
  tm_shape(amph_points) +
  tm_bubbles(
    size = "capacity",
    size.scale = tm_scale_continuous(values = c(0.5, 2)), 
    fill = "prov.type", 
    fill.scale = tm_scale(values = c("red", "blue", "green")),
    fill_alpha = 0.8
  ) +

  # Слой с подписями для амфитеатров с capacity > 30000
  tm_shape(amph_points |> filter(capacity > 30000)) +
  tm_text(
    text = "label",
    options = opt_tm_text(point.label = TRUE)
    ) +

  # Настройка легенды
  tm_layout(
    legend.position = c("right", "bottom"),
    legend.frame = TRUE
  )
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: NA
st_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)

inter
Sparse 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(fill = "count") +
  tm_layout(legend.position = c("right", "bottom"), 
            legend.frame = TRUE)

21.6 tm_basemap()

Парящая в вакууме империя не очень радует глаз; в таком случае стоит добавить растровое изображение ландшафта. Пока это доступно только для динамической карты, поэтому переключаемся в режим “view” (ниже представлен скриншот).

Будьте внимательны, совмещая исторические карты с современными! Убедитесь, что вы не показываете походы Цезаря в современную Швейцарию, как это произошло, например, здесь.

# Устанавливаем режим просмотра (интерактивная карта)
tmap_mode("view")
ℹ tmap mode set to "view".
# Строим интерактивную карту
tm_shape(amph_points) +
  tm_basemap("Stadia.StamenTerrainBackground") +
  tm_symbols(
    size = "capacity",     
    fill = "white",
    col = "steelblue", 
    fill_alpha = 0.8
  )

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, respectively

21.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 Данные: римские дороги

roman_roads <- cawd::darmc.roman.roads.major.sp |> 
  st_as_sf()

roman_roads
ggplot() +
  geom_sf(data = roman_map, fill = "wheat") +
  geom_sf(data = roman_roads,
          color = "steelblue",
          alpha = 0.5) +
  geom_tile()

Добавим подложку.

library(ggspatial)
# Reproject to EPSG:3857 (Web Mercator)
roman_map_3857 <- st_transform(roman_map, 3857)
roman_roads_3857 <- st_transform(roman_roads, 3857)

ggplot() +
  # Add spatial tile background
  annotation_map_tile(
    type = "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/${z}/${y}/${x}.jpg",
    zoomin = -1) +

  # Add your spatial features
  geom_sf(data = roman_map_3857, 
          fill = "wheat", alpha = 0.4) +
  geom_sf(data = roman_roads_3857, 
          color = "darkblue", alpha = 0.8, size = 1) +

  # Set coordinate system to Web Mercator to match tiles
  coord_sf(crs = st_crs(3857)) 
Zoom: 3