7  Работа с геоданными

7.1 Векторная и растровая графика

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

  • Растровые изображения представляют собой набор упорядоченных пикселей, про каждый из которых хранится информация о цвете. Векторное изображение нельзя бесконечно увеличивать — в какой-то момент станут видны пиксели, которые в каком-то смысле являются пределом увеличения. Наиболее популярные форматы растровых изображений: JPEG, GIF, PNG, BMP, TIFF и другие.
  • В векторных изображениях информация хранится как собрание точек, линий и полигонов в некоторой системе координат, что позволяет бесконечно увеличивать такие изображения не теряя в качестве. Наиболее популярные форматы векторных изображений: PDF, SVG, EPS и другие.

Современные технологии позволяют соединять растровые и векторные изображения, а также трансформировать их друг в друга. Картографические данные могут попадать в разные типы: точки (столицы всех стран), линии (улицы в каком-нибудь городе), полигоны (границы стран и меньших регионов) обычно имеют некоторую геопривязку (для простоты давайте считать таким все, что имеет широту и долготу), так что могут быть представлены векторно, однако существует достаточно много информации, которую невозможно представить никак по-другому, кроме как векторно: спутниковые снимки, существующие физические/политические/климатические/исторические и т. п. карты, выдача картографических сервисов, таких как Google Maps. Кроме того, занимаясь любыми типами визуализации, следует помнить о разнице статической визуализации, которую после создания нельзя изменить, и динамической визуализации, которая позволяет настроить отоброжение по желанию пользователя (увеличивать и уменьшать, кликать на собрание точек и видеть их значения и т. п.). В данной главе, в отличие от предыдущих, мы сосредоточимся на пакете для динамического картографирования leaflet. Достаточно много тем останется за пределами этой главы: изменение проекции, манипуляции с географическими данными, работа с растровыми изображениями и другие (см., например, (Lovelace, Nowosad, and Muenchow 2019), доступная он-лайн).

7.2 Картографические примитивы

В картографии существуют свои элементарные единицы:

Географические примитивы из (Lovelace, Nowosad, and Muenchow 2019) (CC BY-NC-ND 4.0)

Эти единицы поддерживают популярные пакеты для манипуляции с георграфическими объектами: sp, sf и другие. В данном разделе мы не будем учиться операциям с этими объектами (объединение, вычитание и т. п., подробности смотрите в документации к пакету sp или в (Lovelace, Nowosad, and Muenchow 2019)).

7.3 Пакет leaflet

Мы пойдем необычным путем и начнем с инструмента, который создает динамические карты — пакета leaflet, который является оберткой для одноименного популярного пакета для визуализации карт в интернете на JS.

Для начала включим библиотеки:

library("leaflet")
library("tidyverse")

Здесь доступен cheatsheet, посвященный пакету leaflet.

7.3.1 .csv файлы

Источником географических данных могут быть обычные привычные нам csv файлы. Например, вот здесь хранится датасет из проекта The Unwelcomed Мохамада А. Вэйкда (Mohamad A. Waked), содержащий информацию о месте и причинах смерти мигрантов и беженцев по всему миру с января 2014 года по июнь 2019 года.

unwelcomed <- read_csv("https://raw.githubusercontent.com/agricolamz/daR4hs/main/data/w6_death_of_migrants_and_refugees_from_the_Unwelcomed_project.csv")
  • id — идентификационный номер;
  • date — дата происшедшего;
  • total_death_missing — количество погибших/пропавших;
  • location — место происшедшего;
  • lat — широта;1
  • lon — долгота;
  • collapsed_region — обобщенная информация о регионе;
  • region — информация о регионе;
  • collapsed_cause — обобщенная информация о причине смерти;
  • cause_of_death — информация о причине смерти.

Самый простой способ нанести на карту координаты, это использовать комбинацию функций leaflet() |> addCircles():

unwelcomed |>  
  leaflet() |>  
  addCircles(lng = ~lon, # обратите внимание на особый синтаксис с тильдой
             lat = ~lat)

7.3.2 Формат .geojson

Существует несколько форматов, в которых принято распространять картографические данные, и если точки удобно хранить в .csv формате, то с полигонами и линиями tidy подход одно наблюдение – одна строчка не подходит. Наиболее распространенными являются .geojson и .shp. Формат .geojson можно прочитать при помощи функции read_json() из пакета jsonlite (я вызываю эту функцию, не загружая пакета, так как пакет jsonlite конфликтует с tidyverse):

moscow_districts <- jsonlite::read_json("https://raw.githubusercontent.com/agricolamz/daR4hs/main/data/w7_moscow.geojson")

leaflet() |> 
  addTiles() |> 
  addGeoJSON(geojson = moscow_districts)

Необходимо зазумиться, так как при отображении полигонов зум не происходит автоматически, этого можно добиться при помощи функции setView():

leaflet() |> 
  addTiles() |> 
  addGeoJSON(geojson = moscow_districts) |> 
  setView(zoom = 8, lng = 37.35, lat = 55.65)

Кроме того, .geojson и .shp можно прочитать функцией st_read() из пакета sf.

7.3.3 Функции пакета leaflet()

Чтобы точки не “висели в воздухе” можно добавить подложку:

unwelcomed |> 
  leaflet() |>  
  addTiles() |>  
  addCircles(lng = ~lon,
             lat = ~lat)

Функция addCircles() имеет массу аргументов, которые отвечают за отображение:

  • radius
  • color
  • opacity
  • fill
  • fillColor
  • label
  • popup

К сожалению, в пакете leaflet нет такого удобного автоматического раскрашивания по некоторой переменной, поэтому для решения такой задачи нужно сначала создать свою функцию раскрашивания. Это делается при помощи функций colorNumeric(), colorFactor(), colorBin() или colorQuantile().

pal_cat <- colorFactor("Set3", domain = unwelcomed$collapsed_cause)
pal_cat(unwelcomed$collapsed_cause[1])
[1] "#D9D9D9"

Теперь в переменную pal_cat записана функция, которая возварщает цвета в зависимости от значения. В качестве первого аргумента в фукнций colorNumeric(), colorFactor(), colorBin() или colorQuantile() отправляется палитра, которую пользователь может задать сам или использовать уже имеющуюся (их можно посмотреть при помощи команды RColorBrewer::display.brewer.all()):

RColorBrewer::display.brewer.all()

Теперь мы готовы сделать нашу первую осмысленную карту

unwelcomed |> 
  filter(str_detect(date, "2014")) |> 
  leaflet() |> 
  addTiles() |> 
  addCircles(lng = ~lon,
             lat = ~lat,
             label = ~total_death_missing, # пусть возникает подпись с количеством
             color  = ~pal_cat(collapsed_cause), # это обобщенная причина
             opacity = 0.9,
             popup = ~cause_of_death) |>  # а это конкретная причина, появляется при клике мышкой
  addLegend(pal = pal_cat,
            values = ~collapsed_cause,
            title = "")

Вообще цветовая схема не очень сочетается с подложкой, так что можно поменять подложку при помощи функции addProviderTiles() (галлерею подложек можно посмотреть вот здесь):

unwelcomed |> 
  filter(str_detect(date, "2014")) |> 
  leaflet() |> 
  addProviderTiles("Esri.WorldPhysical") |> 
  addCircles(lng = ~lon,
             lat = ~lat,
             label = ~total_death_missing, # пусть возникает подпись с количеством
             color  = ~pal_cat(collapsed_cause), # это обобщенная причина
             opacity = 0.9,
             popup = ~cause_of_death) |>  # а это конкретная причина, появляется при клике мышкой
  addLegend(pal = pal_cat,
            values = ~collapsed_cause,
            title = "")

7.3.4 Комбинация карт: leafsync

Карты, как и все объекты в R, тоже можно записать в переменную:

unwelcomed |> 
  filter(str_detect(date, "2014")) |> 
  leaflet() |> 
  addTiles() |> 
  addCircles(lng = ~lon,
             lat = ~lat,
             label = ~total_death_missing, # пусть возникает подпись с количеством
             color  = ~pal_cat(collapsed_cause), # это обобщенная причина
             opacity = 0.9,
             popup = ~cause_of_death) |>  # а это конкретная причина, появляется при клике мышкой
  addLegend(pal = pal_cat,
            values = ~collapsed_cause,
            title = "2014") ->
  m_2014

Теперь, если вызвать переменную m_2014, появится карта, которую мы сделали. Но что, если мы хотим отобразить рядом карты 2014 года и 2015 года? Как сделать фасетизацию? К сожалению, функции для фасетизации в пакете не предусмотрено, но мы можем сделать ее самостоятельно. Для начала создадим вторую карту:

unwelcomed |> 
  filter(str_detect(date, "2015")) |> 
  leaflet() |> 
  addTiles() |> 
  addCircles(lng = ~lon,
             lat = ~lat,
             label = ~total_death_missing, # пусть возникает подпись с количеством
             color  = ~pal_cat(collapsed_cause), # это обобщенная причина
             opacity = 0.9,
             popup = ~cause_of_death) |>  # а это конкретная причина, появляется при клике мышкой
  addLegend(pal = pal_cat,
            values = ~collapsed_cause,
            title = "2015") ->
  m_2015

Включим библиотеку:

library(leafsync)

И теперь соединим две карты:

sync(m_2014, m_2015)

7.3.5 Работа с полигонами и полилиниями

Для добавления на карты линии следует использовать функцию addPolylines():

moscow_metro <- read_csv("https://raw.githubusercontent.com/agricolamz/daR4hs/main/data/w7_moscow_metro.csv")

moscow_metro |> 
  filter(line_name == "сокольническая") |> 
  leaflet() |> 
  addTiles() |> 
  addPolylines(lng = ~longitude,
               lat = ~latitude, 
               color = "tomato") |> 
  addCircles(lng = ~longitude,
             lat = ~latitude)

Для добавления на карты полигона следует использовать функцию addPolygons():

moscow_metro |> 
  filter(line_name == "кольцевая") |> 
  leaflet() |> 
  addTiles() |> 
  addPolygons(lng = ~longitude,
              lat = ~latitude, 
              color = "brown") |> 
  addCircles(lng = ~longitude,
             lat = ~latitude,
             label = ~name)

К сожалению, leaflet не такой удобный, как, скажем, ggplot2, поэтому для того, чтобы нарисовать много отдельных линий или полигонов, нужно использовать цикл:

moscow_metro |> 
  leaflet() |> 
  addTiles() ->
  moscow_map

walk(unique(moscow_metro$line_name), function(i){
  moscow_map <<- moscow_map |> 
    addPolylines(lng = ~longitude,
                 lat = ~latitude, 
                 label = ~line_name, 
                 weight = 2,
                 data = moscow_metro |> filter(line_name == i)) |> 
    addCircles(lng = ~longitude,
               lat = ~latitude,
               label = ~name, 
               data = moscow_metro |> filter(line_name == i))
})

moscow_map

7.4 ggplot2, maps и другие пакеты

Карты можно строить и в статическом ggplot2. Для этого обычно используют пакет maps, в котором хранится датасет с полигонами стран.

library(maps)

Функция map_data() позволяет достать полигоны стран в табличном формате. Я обычно убираю Антарктику, так как на ней не так часто что-то происходит.

map_data("world") |> 
  filter(region != "Antarctica") ->
  world

world |> 
  ggplot(aes(long, lat))+
  geom_map(map = world, aes(map_id = region), color = "grey80", fill = "grey95")+
  geom_point(data = unwelcomed, aes(lon, lat),  alpha = 0.5, size = 0.2)+
  coord_quickmap()+
  theme_void()

Для того, чтобы рисовать карту, я использую geom_map(), который позволяет нанести полигоны стран (обратите внимание на аргумент map в функции geom_map и аргумент map_id в функции aes()). Функция coord_quickmap() позволяет сделать стандартную проекцию (см. ниже). Теперь мы можем легко использовать стандартные средства ggplot2, например, фасетизацию.

world |> 
  ggplot(aes(long, lat))+
  geom_map(map = world, aes(map_id = region), color = "grey80", fill = "grey95")+
  geom_point(data = unwelcomed, aes(lon, lat),  alpha = 0.5, size = 0.2)+
  coord_quickmap()+
  theme_void()+
  facet_wrap(~collapsed_cause)

7.4.1 Проекции

Земля имеет очень сложную форму, поэтому редукция ее до 2D пространства — это всегда достаточно сложная математическая операция по замене геоида (в случае Земли) на какую-то другую геометрическую фигуру и последующую развертку получившегося на плоскость. При любой проекции есть искажения длин, углов, площадей или форм. Есть большая классификация проекций, в которую мы не будем вдаваться. Мы ограничимся лишь информацией, что изменение проекции происходит в функции coord_map() и в справке к этой функции можно посмотреть доступные проекции (coord_map).

states <- map_data("state")

states |> 
  ggplot(aes(long, lat)) +
  geom_polygon(aes(group = group)) +
  coord_map("albers",  lat0 = 45.5, lat1 = 29.5)


  1. Информация о широте и долготе иногда записывают в градусах, минутах и секундах, а иногда в десятичной записи, в R обычно используется десятичная запись. В интернете легко найти конвертеры из одного формата в другой и обратно.↩︎