Тема 5: Географические данные в R

Лекция 8. Работа с географическими данными в R

Большие данные и аналитика
Бизнес-информатика 38.03.05 (бакалавриат)

Источники

Книги

Обзор библиотек

Примеры

OpenStreetMap

Анализ достижимости

Облако точек (LiDAR)

Jean-Romain Roussel, Tristan R.H. Goodbody, Piotr Tompalski The lidR package

3D-модели

3D-модели

library(rayshader)
library(rayvista)

HopkinsNZ <- plot_3d_vista(lat = -44.042238, 
                           long = 169.860985, 
                           radius = 5000, 
                           overlay_detail = 14,
                           elevation_detail = 13, 
                           zscale = 5, 
                           theta = 25, phi = 25, zoom = 0.6,
                           windowsize = 1200, 
                           solid = T, 
                           background = "white")

render_depth(focus = 0.6, focallength = 15, clear=TRUE)

Leaflet

Leaflet

library(leaflet)

lng_ITMO <- 30.3108933
lat_ITMO <- 59.9565684

leaflet() |>
  addTiles() |> 
  # центр и масштаб
  setView(lng = lng_ITMO, 
          lat = lat_ITMO, zoom = 14) |>
  # маркер
  addMarkers(lng = lng_ITMO, 
             lat = lat_ITMO,
             popup = "Университет ИТМО") |>
  # линейка масштабирования
  addScaleBar()

MapGL

library(mapgl)

mapboxgl() |> 
  fly_to(
    center = c(37.530728, 
               55.703134),
    zoom = 15.4,
    bearing = 0, 
    pitch = 60
  )

Тематические карты

Типы географических данных

Модели пространственных данных

  • Векторная модель пространственных данных включает описание координатных данных пространственных объектов и, опционально, топологических отношений между ними (точки, линии, полигоны, объемные тела…).
  • Растровая модель описывает не объекты, а пространственное распределение некоторой (выбранной исследователем) характеристики. Пространство разбивается регулярной сеткой ячеек, в каждой ячейке фиксируется значение исследуемого параметра.

Самсонов Т.Е. Визуализация и анализ географических данных на языке R
(глава Пространственные данные). М.: Географический факультет МГУ, 2024.
DOI: 10.5281/zenodo.901911

Типы файлов

  • Векторные данные
    • Shapefile (.shp): избегаем! (см. Switch from Shapefile)
    • GeoPackage (.gpkg): общий тип данных
    • GeoJSON (.geojson): для веб-приложений
    • GeoParquet (.parquet): большие данные
  • Растровые данные
    • GeoTIFF (.tif): наиболее используемый формат
    • Zarr (.zarr): данные, оптимизированные для облачных вычислений
  • Облако точек
    • LAS (.las): наиболее используемый формат
    • LAZ (.laz): сжатая версия LAS

Векторные данные

Линейно связные модели пространственных объектов

Тип Описание
POINT нульмерная геометрия, содержащая одну точку
LINESTRING последовательность точек, соединенных прямыми, несамопересекающимися отрезками
POLYGON двумерная геометрия с положительной площадью; последовательность точек, отрезки между которыми формируют замкнутое кольцо без самопересечений
MULTIPOINT множество точек
MULTILINESTRING множество линий
MULTIPOLYGON множество полигонов
GEOMETRYCOLLECTION множество геометрий произвольного типа

Это 7 основных моделей из 17 возможных.

Формат WKT

Well-Known Text (WKT) — стандарт представления геометрии в виде множества списков координат, в которых координаты вершин разделены пробелами, вершины разделены запятыми, а компоненты полигонов и мультигеометрий заключены в круглые скобки и также разделены запятыми.

Примеры представлений:

# POINT (0.5 0.5)
# LINESTRING (0 1, 0.5 1.5, 1.2 1.2, 2 1.3, 3 2)
# POLYGON ((0.5 0.5, 2 0, 3 2, 1.5 4, 0 3, 0.5 0.5), (1 1, 0.8 2, 2 2.2, 1.4 1.1, 1 1))
# MULTIPOINT ((0.5 0.5), (1 3), (2 1), (0.2 2), (2 3), (1.5 1.5))
# MULTILINESTRING ((0.5 1.5, 1.2 1.2, 2 1.3), (0 1.5, 0.5 2, 1.2 1.7), (2 1.8, 3 2.5))

Пример данных

В качестве исходных данных выберем North Carolina SIDS data.

nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), 
                  quiet = TRUE)
Simple feature collection with 5 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS:  NAD27
   AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
1 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
2 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
3 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
4 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
5 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
  NWBIR74 BIR79 SID79 NWBIR79                       geometry
1      10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
2      10   542     3      12 MULTIPOLYGON (((-81.23989 3...
3     208  3616     6     260 MULTIPOLYGON (((-80.45634 3...
4     123   830     2     145 MULTIPOLYGON (((-76.00897 3...
5    1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...

Базовая карта

library(sf)

ggplot(nc) +
  geom_sf(aes(fill = AREA),
          color = "black",
          linewidth = 0.2) +
  labs(fill = "площадь") +
  hrbrthemes::theme_ipsum() +
  viridis::scale_fill_viridis() + 
  guides(fill = guide_colorbar(title.position = "top", 
                               title.hjust = 0.5,
                               barwidth = unit(20, "lines"), 
                               barheight = unit(0.7, "lines"))) +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14))

Базовая карта

Пусть цвет соответствует площади полигонов графств в градусных единицах.

Базовая карта

library(mapview)
mapview(nc, zcol = "AREA")

Natural Earth

Natural Earth — это открытые мелкомасштабные картографические данные высокого качества. Данные доступны для трех масштабов: 1:10М, 1:50М и 1:110М, доступные с помощью библиотеки {rnaturalearth}.

Сделаем оценку населения по странам.

library(rnaturalearth)
world <- ne_countries(returnclass = "sf")

ggplot(world) + 
  # отображение объектов типа sf
  # цвет = логарифму оценке населения
  geom_sf(aes(fill = log10(pop_est))) + 
  # картографическая проекция
  coord_sf(crs = "+proj=eqearth") +
  viridis::scale_fill_viridis(option = "turbo", 
                              direction = -1) +
  hrbrthemes::theme_ipsum() +                            
  theme(legend.position = "none")

Natural Earth

Spatial Data Sources in R (Cédric Scherer)

Растровые данные

Растр представляет из себя матрицу значений. Каждой ячейке матрицы соответствует прямоугольная пространственная область фиксированного размера, которая называется пикселом (см. курс Основы геоинформатики).

Растровые данные

Растровые карты могут представлять непрерывные данные, такие как высота над уровнем моря, температура, количество осадков.

Карты также могут представлять категорные данные, такие как виды почвы или классы растительного покрова.

Растровые данные

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

Пример растровой карты

Рассмотрим среднемесячную температуру в градусах Цельсия глобальных данных о температуре из базы данных WorldClim, указав:

страну (country = "Liechtenstein")

переменную — среднюю температуру (var = "tavg")

разрешение (res = 10)

путь загрузки данных в качестве временного файла (path = tempdir()).

library(terra)
library(geodata)

r <- geodata::worldclim_country(country = "Liechtenstein", 
                                var = "tavg",
                                res = 10, 
                                path = tempdir())
plot(r)

Пример растровой карты

Географические проекции

Проекции

Проекция — это математический способ перехода от географических координат на эллипсоиде к плоским прямоугольным координатам карты.

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

Страница List of all projection images содержит все проекции.

Сайт The True Size of… позволяет сравнивать реальные размеры стран.

Системы отсчета координат (CRS)

  • Система координатной привязки (CRS) пространственных данных определяет начало и единицу измерения пространственных координат.
  • В географической CRS значения широты и долготы используются для определения местоположений на трехмерной эллипсоидной поверхности Земли. Проецируемые CRS используют декартовы координаты для ссылки на местоположение на двумерном представлении Земли.
  • Все проекции искажают поверхность Земли в той или иной степени и не могут одновременно сохранять все свойства площади, направления, формы и расстояния.

Проекция Миллера

my_map + coord_sf(crs = "+proj=mill")

Проекция Атласа Times

my_map + coord_sf(crs = "+proj=times")

Проекция Equal Earth

my_map + coord_sf(crs = "+proj=eqearth")

Проекция Эккерта III

my_map + coord_sf(crs = "+proj=eck3")

Проекция Мольвейде

my_map + coord_sf(crs = "+proj=moll")

EPSG-коды

  • Наиболее распространенные CRS можно указать, указав их коды EPSG (European Petroleum Survey Group) — европейская рабочая группа нефтегазовой области, которая ведет реестр систем координат с уникальными цифровыми кодами вида EPSG:xxxxxx.
  • Выбрать проекцию можно с помощью сайта Projection Wizard. Все распространенные пространственные проекции можно найти на странице Spatial Reference List.
  • Например, WGS 84 (англ. World Geodetic System 1984) — всемирная система геодезических параметров Земли 1984 года, в число которых входит система геоцентрических координат, соответствует коду EPSG 4326.

Данные OpenStreetMap

Начало работы с OSM

OpenStreetMap (OSM) — самая крупная открытая база ГИС-данных.

Данные OSM можно скачать и затем использовать локально, воспользовавшись сервером загрузки Geofabrik.

Для работы с OSM данными можно использовать библиотеки {osmdata} и {osrm}.

# работа с географическими данными
library(osmdata)
library(osrm)
library(sf)

# шкала масштаба на карте
library(ggspatial)

Административные границы

# загрузим все границы, связанные с объектом
my_place <- "Выборг"

all_boundaries <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "boundary",
                  value = "administrative") |>
  osmdata_sf() |>
  unname_osmdata_sf() %>%
  .$osm_multipolygons
1
делает Overpass-запрос
2
добавляет функции в запрос Overpass, которые задаются парами key-value; подробное описание значений представлено на wiki OSM
3
возвращает sf-объект
4
удаляет имена из osmdata геометрических объектов

Все границы

all_boundaries |>
  select(name)
Simple feature collection with 16 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 27.50471 ymin: 60.00596 xmax: 30.10531 ymax: 61.32925
Geodetic CRS:  WGS 84
First 10 features:
                                      name                       geometry
299498  Каменногорское городское поселение MULTIPOLYGON (((28.83912 60...
312143     Гончаровское сельское поселение MULTIPOLYGON (((28.99518 60...
312144                              Выборг MULTIPOLYGON (((28.82928 60...
391526       Советское городское поселение MULTIPOLYGON (((28.99051 60...
1162550                   Выборгский район MULTIPOLYGON (((29.56322 61...
1290892    Селезнёвское сельское поселение MULTIPOLYGON (((27.50471 60...
2136991             Калининский микрорайон MULTIPOLYGON (((28.6809 60....
2136992               Кировский микрорайон MULTIPOLYGON (((28.8467 60....
2136993           Петербургский микрорайон MULTIPOLYGON (((28.87301 60...
2136994              Петровский микрорайон MULTIPOLYGON (((28.66611 60...

Визуализация границ

Выделим из полученного множества границ только нужную нам, соответствующую границе городского округа.

boundary_Vyborg <- all_boundaries |>
  dplyr::filter(osm_id == 312144)

ggplot() +
  geom_sf(data = all_boundaries, 
          alpha = 0.001) +
  geom_sf(data = boundary_Vyborg, 
          color = "red",
          alpha = 0.001)

Визуализация границ

Типы загружаемых дорог

# типы загружаемых дорог и толщина линий на карте согласно OpenStreetMap
highway_sizes <- tibble::tribble(
  ~highway, ~highway_group,
  "motorway",      "large",
  "trunk",         "large",
  "motorway_link", "large",
  "primary",       "large",
  "primary_link",  "large",
  "secondary",     "medium",
  "secondary_link","medium",
  "tertiary",      "medium",
  "tertiary_link", "medium",
  "trunk_link",    "medium",
  "residential",   "small",
  "living_street", "small",
  "unclassified",  "small",
  "service",       "small",
  "road",          "small",
)

Загрузка дорог

# streets (улицы)
# учитывается тип дорожного покрытия
streets_osm <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "highway", 
                  value = highway_sizes$highway) |>
  osmdata_sf() |> 
  unname_osmdata_sf()

streets <- streets_osm$osm_lines %>%
  mutate(length = as.numeric(st_length(.))) |>
  left_join(highway_sizes, by = "highway")

Водные объекты

# water (водные объекты)
water_osm <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "natural", 
                  value = "water") |>
  osmdata_sf() |> 
  unname_osmdata_sf()

river_osm <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "waterway", 
                  value = c("river", 
                            "riverbank")) |>
  osmdata_sf() |> 
  unname_osmdata_sf()

water <- 
  c(water_osm, river_osm) %>%
  .$osm_multipolygons |>
  dplyr::select(osm_id, name)

Железные дороги

# railways (железные дороги)
railways_osm <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "railway", 
                  value = "rail") |>
  osmdata_sf() |> 
  unname_osmdata_sf()

railways <- railways_osm$osm_lines |> 
  dplyr::select()

Базовая карта

base_map <- 
  ggplot() + 
  geom_sf(data = boundary_Vyborg,
          alpha = 0.001, linewidth = 0.8, color = "red") +
  geom_sf(data = streets |>
            dplyr::filter(highway_group == "large"),
          linewidth = 0.4, color = "grey30") +
  geom_sf(data = streets |>
            dplyr::filter(highway_group == "medium"),
          linewidth = 0.2, color = "grey35") +
  geom_sf(data = streets |>
            dplyr::filter(highway_group == "small"),
          linewidth = 0.1, color = "grey40") +
  geom_sf(data = water,
          fill = "#A5CFE9",
          lwd = 0, alpha = 0.5) +
  geom_sf(data = railways,
          color = "grey30", linewidth = 0.3, linetype = "dotdash", alpha = 0.6)
1
границы объекта
2
дорожная сеть
3
водные объекты
4
железные дороги

Базовая карта

Ограничим карту на прямоугольник-подмножество.

base_map +
  # тема
  hrbrthemes::theme_ipsum() +
  # географические границы части города
  coord_sf(xlim = c(28.6, 28.9), 
           ylim = c(60.65, 60.82),
           expand = FALSE)

Базовая карта

Точки интерса (POI)

attractions_tibble <- tibble(
  attractions = c("Выборгский замок", "Парк Монрепо", "Часовая башня", 
                  "Круглая башня", "Библиотека Алвара Аалто"),
  lon = c(28.7295, 28.7364, 28.7321, 28.7333, 28.7389),
  lat = c(60.7159, 60.7234, 60.7136, 60.7142, 60.7112)
)

# создание sf-объекта
attractions_tibble_sf <- attractions_tibble |>
  st_as_sf(coords = c("lon", "lat"), 
           crs = 4326)

Точки интерса (POI)

attractions_tibble
# A tibble: 5 × 3
  attractions               lon   lat
  <chr>                   <dbl> <dbl>
1 Выборгский замок         28.7  60.7
2 Парк Монрепо             28.7  60.7
3 Часовая башня            28.7  60.7
4 Круглая башня            28.7  60.7
5 Библиотека Алвара Аалто  28.7  60.7

Точки интерса (POI)

attractions_tibble_sf
Simple feature collection with 5 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 28.7295 ymin: 60.7112 xmax: 28.7389 ymax: 60.7234
Geodetic CRS:  WGS 84
# A tibble: 5 × 2
  attractions                      geometry
* <chr>                         <POINT [°]>
1 Выборгский замок        (28.7295 60.7159)
2 Парк Монрепо            (28.7364 60.7234)
3 Часовая башня           (28.7321 60.7136)
4 Круглая башня           (28.7333 60.7142)
5 Библиотека Алвара Аалто (28.7389 60.7112)

Карта объектов

base_map +
  geom_sf(data = attractions_tibble_sf, 
          shape = 8, size = 2.5, stroke = 1.1, color = "red") +
  # географические границы части города
  coord_sf(xlim = c(28.72, 28.77), 
           ylim = c(60.705, 60.73),
           expand = FALSE) +
  # масштабная линейка
  annotation_scale(location = "tr", 
                   width_hint = 0.5, 
                   style = "ticks") +
  # подписи
  ggrepel::geom_label_repel(data = attractions_tibble, 
                            aes(lon, lat, label = attractions), 
                            fontface = "bold",
                            size = 4, alpha = 0.7) +
  labs(x = "", y = "")

Карта объектов

Что мы изучили

Рассмотренные вопросы

  • Основные модели географических данных (векторные, растровые)
  • Географические проекции
  • Работа с сетевыми данными OpenStreetMap

Спасибо за внимание!