Основы картографии

class: center, middle, inverse, title-slide

.title[
# Основы картографии
]
.subtitle[
## Визуализация и анализ географических данных на языке R
]
.author[
### Тимофей Самсонов
]
.institute[
### МГУ имени Ломоносова, Географический факультет
]
.date[
### 2023-11-07
]

---

## Используемые пакеты

```r
library(sf)
library(stars)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(ggnewscale)
library(rnaturalearth)
library(rmapshaper)
library(RColorBrewer)
```

---

## Источники данных для картографирования

- Сверхмелкие масштабы (10M и мельче): __Natural Earth__ (мир)

- Мелкие масштабы (1 - 10M): __ВСЕГЕИ__ (Россия)

- Средние масштабы (100k - 1M): .red[__?__]

- Крупные масштабы (крупнее 100k): __OpenStreetMap (OSM)__ (мир)

---

## Загрузка данных Natural Earth

Данные NE могут быть выгружены онлайн:

```r
countries = ne_countries(scale = 110, returnclass = 'sf')
coast = ne_coastline(scale = 110, returnclass = 'sf')

ocean = ne_download(scale = 110,
                    type = 'ocean',
                    category = 'physical',
                    returnclass = 'sf')
## OGR data source with driver: ESRI Shapefile 
## Source: "/private/var/folders/5s/rkxr4m8j24569d_p6nj9ld200000gn/T/Rtmpl6QJAh", layer: "ne_110m_ocean"
## with 2 features
## It has 3 fields

cities = ne_download(scale = 110,
                     type = 'populated_places',
                     category = 'cultural',
                     returnclass = 'sf')
## OGR data source with driver: ESRI Shapefile 
## Source: "/private/var/folders/5s/rkxr4m8j24569d_p6nj9ld200000gn/T/Rtmpl6QJAh", layer: "ne_110m_populated_places"
## with 243 features
## It has 137 fields
## Integer64 fields read as strings:  POP_MAX POP_MIN POP_OTHER MAX_POP10 MAX_POP20 MAX_POP50 MAX_POP300 MAX_POP310 MAX_NATSCA MIN_AREAKM MAX_AREAKM MIN_AREAMI MAX_AREAMI MIN_PERKM MAX_PERKM MIN_PERMI MAX_PERMI POP1950 POP1955 POP1960 POP1965 POP1970 POP1975 POP1980 POP1985 POP1990 POP1995 POP2000 POP2005 POP2010 POP2015 POP2020 POP2025 POP2050 WOF_ID NE_ID GEONAMESID
```

---

## Загрузка данных Natural Earth

В то же время, более эффективно использовать локальную копию данных:

```r
ne = '/Volumes/Data/Spatial/Natural Earth/natural_earth_vector.gpkg'
rivers = st_read(ne, 'ne_110m_rivers_lake_centerlines', quiet = TRUE)
lakes = st_read(ne, 'ne_110m_lakes', quiet = TRUE)
land = st_read(ne, 'ne_110m_land', quiet = TRUE)
borders = st_read(ne, 'ne_110m_admin_0_boundary_lines_land', quiet = TRUE)

lyr110 = lst(ocean, land, coast, countries, 
             rivers, lakes, cities, borders)
```

---

## Визуализация средствами ggplot2

- `geom_sf()` вызывает `stat_sf()` и `coord_sf()` чтобы отобразить объекты типа `sf` в нужной системе координат;
- `geom_stars()` отображает объекты типа `stars`;
- `coord_sf()` обеспечивает поддержку картографических проекций и позволяет отображать данные в нужной системе координат на лету;
- `stat_sf()` отвечает за отображение переменных данных на графические переменные для пространственных данных;
- `geom_sf_label()` позволяет отображать подписи объектов на плашке;
- `geom_sf_text()` позволяет размещать подписи объектов без плашки.

---

## Визуализация средствами ggplot2

Пример — политическая карта мира

```r
ggplot() +
  geom_sf(data = lyr110$countries, color = NA, 
          mapping = aes(fill = as.factor(mapcolor7)), 
          show.legend = FALSE) +
  scale_fill_manual(values = brewer.pal(7, 'Set2')) +
  geom_sf(data = lyr110$borders, size = 0.2) +
  geom_sf(data = lyr110$ocean, fill = 'azure', color = NA) +
  geom_sf(data = lyr110$coast, size = 0.4, color = 'steelblue4') +
  theme_void()
```

---

## Визуализация средствами ggplot2

Пример — политическая карта мира

![](10_BaseMaps_files/figure-html/unnamed-chunk-4-1.png)

---

## Визуализация средствами ggplot2

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

```r
lyr110$megacities = lyr110$cities |> 
  filter(SCALERANK == 0, ! NAME %in% c('Washington, D.C.', 
        'Paris', 'Riyadh', 'Rome', 'São Paulo', 'Kolkata'))

basemap = list(
  geom_sf(data = lyr110$countries, color = NA, 
          mapping = aes(fill = as.factor(mapcolor7)), 
          show.legend = FALSE),
  scale_fill_manual(values = brewer.pal(7, 'Set2')),
  geom_sf(data = lyr110$borders, size = 0.2),
  geom_sf(data = lyr110$ocean, fill = 'azure', color = NA),
  geom_sf(data = lyr110$coast, size = 0.4, color = 'steelblue4'),
  geom_sf(data = lyr110$megacities, shape = 21, 
          fill = 'white', stroke = 0.75, size = 2)
)
```

---

## Подписи объектов

Полученный список теперь можно использовать в качестве карты-подложки. Нанесем точки населенных пунктов:

```r
ggplot() +
  basemap +
  geom_sf_text(data = lyr110$megacities, 
               mapping = aes(label = NAME_EN),
               size = 3, nudge_y = 5, 
               family = 'Open Sans', fontface = 'bold') +
  theme_void()
```

---

## Подписи объектов

Полученный список теперь можно использовать в качестве карты-подложки. Нанесем точки населенных пунктов:

![](10_BaseMaps_files/figure-html/unnamed-chunk-6-1.png)

---

## Подписи объектов

Полученный список теперь можно использовать в качестве карты-подложки. Нанесем точки населенных пунктов:

```r
ggplot() +
  basemap +
  geom_text_repel(data = lyr110$megacities, stat = "sf_coordinates",
                  size = 3, aes(label = NAME, geometry = geometry), 
                  family = 'Open Sans', fontface = 'bold') +
  theme_void()
```

---

## Подписи объектов

Полученный список теперь можно использовать в качестве карты-подложки. Нанесем точки населенных пунктов:

![](10_BaseMaps_files/figure-html/unnamed-chunk-7-1.png)

---

## Улучшение читаемости подписей

Полученный список теперь можно использовать в качестве карты-подложки. Нанесем точки населенных пунктов:

```r
basemap0 = list(
  geom_sf(data = lyr110$countries, color = NA, alpha = 0.5,
          mapping = aes(fill = as.factor(mapcolor7)), show.legend = FALSE),
  scale_fill_manual(values = brewer.pal(7, 'Set2')),
  geom_sf(data = lyr110$borders, alpha = 0.5, size = 0.2),
  geom_sf(data = lyr110$ocean, fill = 'azure', color = NA),
  geom_sf(data = lyr110$coast, alpha = 0.5, size = 0.4, color = 'steelblue4'),
  geom_sf(data = lyr110$megacities, shape = 21, fill = 'white', stroke = 0.75, size = 2)
)
ggplot() +
  basemap0 +
  geom_text_repel(data = lyr110$megacities, stat = "sf_coordinates",
                  size = 3, aes(label = NAME, geometry = geometry), 
                  family = 'Open Sans', fontface = 'bold') +
  theme_void()
```

---

## Улучшение читаемости подписей

Полученный список теперь можно использовать в качестве карты-подложки. Нанесем точки населенных пунктов:

![](10_BaseMaps_files/figure-html/unnamed-chunk-8-1.png)

---

## Улучшение читаемости подписей

Альтернативное решение — можно сделать полупрозрачную плашку под подписью:

```r
ggplot() +
  basemap +
  geom_label_repel(data = lyr110$megacities, stat = "sf_coordinates",
                  aes(label = NAME, geometry = geometry), 
                  size = 3, 
                  label.size = NA, 
                  label.padding=.1, 
                  fill = alpha("white", 0.7), 
                  family = 'Open Sans', fontface = 'bold') +
  theme_void()
```

---

## Улучшение читаемости подписей

Полученный список теперь можно использовать в качестве карты-подложки. Нанесем точки населенных пунктов:

![](10_BaseMaps_files/figure-html/unnamed-chunk-9-1.png)

---

## Проекции

Замечательной особенностью ggplot2 является умени проецировать данные на лету непосредственно при отображении карты. Запишем непроецированную карту в отдельную переменную:

```r
map = ggplot() +
  basemap +
  geom_label_repel(data = lyr110$megacities, stat = "sf_coordinates",
                  aes(label = NAME, geometry = geometry), 
                  size = 3, 
                  label.size = NA, 
                  label.padding=.1, 
                  fill = alpha("white", 0.7), 
                  family = 'Open Sans', fontface = 'bold') +
  theme_void()
```

---

## Проекции

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

```r
map + coord_sf(crs = "+proj=moll")
```

---

## Проекции

Проекция __Эккерта III__:

```r
map + coord_sf(crs = "+proj=eck3")
```

---

## Проекции

Проекция __Equal Earth__:

```r
map + coord_sf(crs = "+proj=eqearth")
```

---

## Проекции

Проекция __Атласа Times__:

```r
map + coord_sf(crs = "+proj=times")
```

---

## Проекции

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

```r
map + coord_sf(crs = "+proj=mill")
```

---

## Градусные сетки

При создании карт мира целесообразно заготовить сетку параллелей и меридианов, а также контур Земного шара, который будет опоясывать карту. Для того чтобы контур проецировался в виде гладкой кривой, следует уплотнить его вершины. Это можно сделать посредством `smoothr::densify()`:

```r
lons = seq(-180, 180, by = 30)
lats = seq(-90, 90, by = 30)

# Сетка координат
grat = st_graticule(lon = lons, lat = lats)

# Контур Земли
box = st_bbox(c(xmin = -180, xmax = 180, 
                ymax = 90, ymin = -90), 
              crs = st_crs(4326)) |> 
  st_as_sfc() |> 
  smoothr::densify(max_distance = 1) 
```

---

## Градусные сетки

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

```r
degree_labels = function(grat, vjust, hjust, 
                         size, lon = T, lat = T) {
  pts = grat |>  
    st_cast('POINT') |> 
    group_by(degree, type, degree_label) |> 
    filter(row_number() == 1)
    
  list(
    if (lon) geom_sf_text(data = filter(pts, type == 'E'), 
                vjust = vjust, size = size,
                mapping = aes(label = degree_label), parse = TRUE),
    if (lat) geom_sf_text(data = filter(pts, type == 'N'), 
                hjust = hjust, size = size,
                mapping = aes(label = degree_label), parse = TRUE)
  )  
}
```

---

## Градусные сетки

Присоединим новые параметры карты к уже существующим:

```r
basemap_grat = c(basemap, list(
  geom_sf(data = grat, size = 0.1),
  geom_sf(data = box, size = 0.5, fill = NA)
))

map_grat = ggplot() +
  basemap_grat +
  geom_label_repel(data = lyr110$megacities, stat = "sf_coordinates",
                  aes(label = NAME, geometry = geometry), 
                  size = 3, 
                  label.size = NA, 
                  label.padding=.1, 
                  fill = alpha("white", 0.7), 
                  family = 'Open Sans', fontface = 'bold') +
  theme_void()
```

---

## Градусные сетки

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

```r
map_grat + coord_sf(crs = "+proj=moll") +
  degree_labels(grat, vjust = +1.5, hjust = +1.5, size = 3)
```

---

## Градусные сетки

Проекция __Эккерта III__:

```r
map_grat + coord_sf(crs = "+proj=eck3") +
  degree_labels(grat, vjust = +1.5, hjust = +1.5, size = 3)
```

---

## Градусные сетки

Проекция __Equal Earth__:

```r
map_grat + coord_sf(crs = "+proj=eqearth") +
  degree_labels(grat, vjust = +1.5, hjust = +1.5, size = 3)
```

---

## Градусные сетки

Проекция __Атласа Times__:

```r
map_grat + coord_sf(crs = "+proj=times") +
  degree_labels(grat, vjust = +1.5, hjust = +1.5, size = 3)
```

---

## Градусные сетки

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

```r
map_grat + coord_sf(crs = "+proj=mill") +
  degree_labels(grat, vjust = +1.5, hjust = +1.5, size = 3)
```

---

### Отображение растровых данных

На общегеографических картах довольно часто присутствует изображение рельефа. Чтобы добавить его на карту, можно использовать специальный тип геометрии `geom_stars`:

```r
dem = read_stars('../r-geo-course/data/world/gebco.tif') # Цифровая модель рельефа

ggplot() +
  geom_stars(data = dem) +
  geom_sf(data = lyr110$coast, size = 0.4, color = 'white') +
  coord_sf() +
  theme_void()
```

---

### Отображение растровых данных

![](10_BaseMaps_files/figure-html/unnamed-chunk-24-1.png)

---

## Отображение растровых данных

Если рельеф показывается и на сушу и на море, его придется раздлелить на 2 ЦМР, иначе цветовая окраска будет некорректной:

```r
pal = c('navyblue', 'steelblue', 'azure', 'darkslategray', 
        'olivedrab', 'lightyellow', 'firebrick', 'pink', 'white')
val = c(min(dem[[1]]), -4000, -200, 0, 100, 300, 
        1000, 2500, max(dem[[1]])) |> scales::rescale()
hydro_lyrs = list(
  geom_sf(data = lyr110$coast, size = 0.4, color = 'steelblue4'),
  geom_sf(data = lyr110$rivers, size = 0.3, color = 'steelblue4'),
  geom_sf(data = lyr110$lakes, size = 0.3, color = 'steelblue4', fill = 'azure')
)
ggplot() +
  geom_stars(data = dem) +
  scale_fill_gradientn(colours = pal, 
                       values = val) +
  hydro_lyrs +
  coord_sf() + theme_void()
```

---

## Отображение растровых данных

Если рельеф показывается и на сушу и на море, есго придется раздлелить на 2 ЦМР, иначе цветовая окраска будет некорректной:

![](10_BaseMaps_files/figure-html/unnamed-chunk-25-1.png)

---

## Отображение растровых данных

Если рельеф показывается и на сушу и на море, его придется разделить на 2 ЦМР, иначе цветовая окраска будет некорректной:

![](10_BaseMaps_files/figure-html/unnamed-chunk-26-1.png)

---

## Отображение растровых данных

Решение: кропнуть растр полигонами суши и океана, отобразить независимо:

```r
sf_use_s2(FALSE)
dem_land = dem[lyr110$land]
dem_ocean = dem[lyr110$ocean]

map = ggplot() +
  geom_stars(data = dem_ocean) +
  scale_fill_gradientn(
    colours = c('navyblue', 'steelblue4', 'skyblue2', 'azure', 'azure'),
    values = scales::rescale(
      c(min(dem_ocean[[1]], na.rm = T), 
        -4000, -200, 0, 
        max(dem_ocean[[1]], na.rm = T))
    ),
    na.value = NA
  ) +
  new_scale_fill() + # обнуляем цветовую щкалу!
  geom_stars(data = dem_land) +
  scale_fill_gradientn(
    colours = c('darkslategray', 'darkslategray', 'olivedrab', 
                'lightyellow', 'firebrick', 'pink', 'white'), 
    values = scales::rescale(
      c(min(dem_land[[1]], na.rm = T), 
        -50, 100, 300, 1500, 3500, 
        max(dem_land[[1]], na.rm = T)
      )
    ), 
    na.value = NA
  ) +
  hydro_lyrs +
  coord_sf() +
  theme_void()
```

---

## Отображение растровых данных

```r
map + coord_sf(xlim = c(10, 75), ylim = c(20, 50)) + anno
```

![](10_BaseMaps_files/figure-html/unnamed-chunk-28-1.png)

---

## Детализация

При создании карт очень важно обеспечить адекватный уровень детализации данных. Иначе изображение будет загромождено деталиями или наоборот будет неинформативным. Иногда проблему можно решить путем выбора картографической основы подходящей детализации:

```r
cnt010 = st_read(ne, 'ne_10m_admin_0_countries')
## Reading layer `ne_10m_admin_0_countries' from data source 
##   `/Volumes/Data/Spatial/Natural Earth/natural_earth_vector.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 255 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS:  WGS 84
cnt050 = st_read(ne, 'ne_50m_admin_0_countries')
## Reading layer `ne_50m_admin_0_countries' from data source 
##   `/Volumes/Data/Spatial/Natural Earth/natural_earth_vector.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 241 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS:  WGS 84
cnt110 = st_read(ne, 'ne_110m_admin_0_countries')
## Reading layer `ne_110m_admin_0_countries' from data source 
##   `/Volumes/Data/Spatial/Natural Earth/natural_earth_vector.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 177 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84

prj = '+proj=laea +lat_0=50 +lon_0=10'

box = st_bbox(c(xmin = -10, xmax = 33, 
                ymin = 33, ymax = 60),
              crs = st_crs(4326)) |> 
  st_as_sfc() |> 
  st_transform(prj) |> 
  st_bbox()

cnts = list(cnt010, cnt050, cnt110)
scales = c(10, 50, 110)
```

---

## Детализация

![](10_BaseMaps_files/figure-html/unnamed-chunk-30-1.png)

---

## Детализация

![](10_BaseMaps_files/figure-html/unnamed-chunk-31-1.png)

---

## Детализация

![](10_BaseMaps_files/figure-html/unnamed-chunk-32-1.png)

---

## Отбор

Довольно часто плотность объектов превышает необходимую:

```r
cities_eu = st_read(ne, 'ne_10m_populated_places', quiet=T) |> 
  st_transform(prj) |> 
  st_crop(box)

cnt_eu = cnt050 |> 
  st_transform(prj) |> 
  st_crop(box)

ggplot() +
  geom_sf(data = cnt_eu, size = 0.25) +
  geom_sf(data = cities_eu, size = 0.5, color = 'darkviolet') +
  geom_sf_text(data = cities_eu, 
               mapping = aes(label = NAME),
               size = 1.5, nudge_y = 30000) +
  theme_bw()
```

---

## Отбор

![](10_BaseMaps_files/figure-html/unnamed-chunk-33-1.png)

---

## Отбор

Проблему можно решить путём отбора объектов по атрибутам:

```r
capitals = filter(cities_eu, FEATURECLA == 'Admin-0 capital')

ggplot() +
  geom_sf(data = cnt_eu, size = 0.25) +
  geom_sf(data = capitals, size = 1.2, color = 'darkviolet') +
  geom_text_repel(data = capitals, stat = "sf_coordinates",
                size = 2.5, aes(label = NAME, geometry = geom), 
                fontface = 'bold') +
  theme_bw()
```

---

## Отбор

![](10_BaseMaps_files/figure-html/unnamed-chunk-34-1.png)

---

## Классификация

Для классификации объектов необходимо создать соответствующий атрибут, который ранжирует их по величине. Например, по численности населения:

```r
capitals = capitals |>
  mutate(RANK = cut(POP_MAX, c(0, 100000, 1000000, 15000000)),
         SIZE = as.numeric(RANK))

ggplot() +
  geom_sf(data = cnt_eu, size = 0.25) +
  geom_sf(data = capitals, mapping = aes(size = SIZE), colour = "black",
          fill = "white", shape = 21,  stroke = 0.5) +
  scale_size(guide = 'none', range = c(1, 1.5)) +
  ggnewscale::new_scale('size') +
  geom_text_repel(data = capitals, stat = "sf_coordinates",
                aes(label = NAME, geometry = geom, size = SIZE),
                fontface = 'bold') +
  scale_size(guide = 'none', range = c(1.5, 2))+
  labs(size = NULL) +
  
  theme_bw()
```

---

## Классификация

![](10_BaseMaps_files/figure-html/unnamed-chunk-35-1.png)