Растровый анализ

class: center, middle, inverse, title-slide

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

---

## Пакеты

Для выполнения кода данной лекции вам понадобятся следующие пакеты:

``` r
library(sf)
library(stars)
library(geosphere)
library(tidyverse)
library(terra)
library(tidyterra)
library(whitebox)
library(ggnewscale)
options(scipen = 999)
```

---

## Исходные данные

``` r
# Чтение данных
(bed = rast('../r-geo-course/data/etopo1_bed.tif'))
## class       : SpatRaster 
## dimensions  : 360, 720, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : etopo1_bed.tif 
## name        : etopo1_bed
(ice = rast('../r-geo-course/data/etopo1_ice.tif'))
## class       : SpatRaster 
## dimensions  : 360, 720, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : etopo1_ice.tif 
## name        : etopo1_ice
countries = read_sf('../r-geo-course/data/countries.gpkg')
```

---

## Исходные данные

``` r
ggplot() +
  geom_spatraster(data = bed) +
  scale_fill_hypso_tint_c(palette = "gmt_globe") +
  ggtitle('ETOPO Bedrock')
```

---

## Исходные данные

``` r
ggplot() +
  geom_spatraster(data = ice) +
  scale_fill_hypso_tint_c(palette = "gmt_globe") +
  ggtitle('ETOPO Ice surface')
```

---

## Растровая алгебра

__Типы операций:__

- _Локальные_ — анализируется одна ячейка растра или совпадающие в пространстве ячейки нескольких растров

- _Фокальные_ — анализируются все ячейки в окрестности. Фокальные методы алгебры карт также называются _методами анализа соседства_.

- _Зональные_ — анализируются все ячейки в пределах зон, определяемых извне (например, вторым растровым слоем).

- _Глобальные_ — анализируются все ячейки растра.

---

## Локальные операции

Вычисление мощности покровного оледенения

.pull-left[
.code-small[

``` r
ice.depth = ice - bed
ice.depth[ice.depth == 0] = NA

ggplot() +
  geom_spatraster(data = bed) +
  scale_fill_gradient(low = 'black', 
                      high = 'white',
                      guide="none") +
  new_scale_fill() +
  geom_spatraster(data = ice.depth) +
  scale_fill_gradient(low = 'white', 
                      high = 'navyblue', 
                      na.value = "transparent") +
  geom_sf(data = countries, fill = NA) +
  labs(title = 'Мощность покровного оледенения',
         fill = '[м]')
```
]
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-5-1.png)
]

---

## Фокальные операции

В общем случае фиксированная окрестность может иметь различную форму, однако наиболее часто используется квадратная окрестность размером `$3\times3$`:

<div class="figure">
<img src="../r-geo-course/images/raster_neighborhoods.png" alt="Виды растровых окрестностей. Темной точкой выделена анализируемая ячейка" width="50%" />
<p class="caption">Виды растровых окрестностей. Темной точкой выделена анализируемая ячейка</p>
</div>

---

## Сглаживание

``` r
# Вырежем кусок из ЦМР
dem = crop(ice, ext(-120, -75, 10, 40))

# Среднее
wgt = matrix(c(1, 1, 1,
               1, 1, 1,
               1, 1, 1) / 9, 
              nrow = 3)
# на самом деле проще написать так:
# wgt = matrix(1/9, 3, 3), но полная форма записана для наглядности

# выполним обработку ЦМР с помощью фокального фильтра
filtered = focal(dem, w = wgt)
```

---

## Сглаживание

``` r
ggplot() +
  geom_spatraster(data = c(dem, filtered)) +
  scale_fill_hypso_tint_c(palette = "etopo1") + 
  facet_wrap(~lyr)
```

---

## Гауссово сглаживание

``` r
filtered = focal(dem, focalMat(dem, 0.5, "Gauss"))

ggplot() +
  geom_spatraster(data = c(dem, filtered)) +
  scale_fill_hypso_tint_c(palette = "etopo1") + 
  facet_wrap(~lyr)
```

---

## Распознавание границ

.pull-left[

``` r
# Матрица Собеля:
wgt = matrix(c(1, 2, 1,
                0, 0, 0,
               -1,-2,-1) / 4, 
              nrow=3)
filtered = focal(dem, wgt)

# Это поверхность производных:

ggplot() +
  geom_spatraster(data = filtered) +
  scale_fill_gradient2(
    low = 'navyblue', 
     high = 'darkred'
  ) +
  ggtitle('Производная поверхности')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-10-1.png)
]

---

## Распознавание границ

.pull-left[

``` r
# Отберем все ячейки, обладающие 
# высокими значениями производных
faults = abs(filtered) > 1500
faults[faults == 0] = NA

ggplot() +
  geom_spatraster(data = dem) +
  scale_fill_hypso_tint_c(
    palette = "etopo1"
  ) + 
  new_scale_fill() +
  geom_spatraster(data = faults) +
  scale_fill_discrete(
    type = 'black', 
    na.value = "transparent"
  ) +
  ggtitle('Уступы континента')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-11-1.png)
]

---

## Морфометрия рельефа

.pull-left[
Высоты:

``` r
setwd('/Users/tsamsonov/GitHub/r-geo-course/data')
dem = rast('dem_fergana.tif')

ggplot() +
  geom_spatraster(data = dem) +
  scale_fill_hypso_c() +
  labs(title = 'Ферганская долина', 
       fill = 'Высота, [м]')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-12-1.png)
]

---

## Морфометрия рельефа

.pull-left[
Углы наклона:

``` r
slope = terrain(dem, 'slope', 
                unit = 'degrees')

ggplot() +
  geom_spatraster(data = slope) +
  scale_fill_gradient(
    low = 'lightcyan', 
    high = 'darkred'
  ) +
  labs(title = 'Углы наклона', 
       fill = 'Градусы [°]')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-13-1.png)
]

---

## Морфометрия рельефа

.pull-left[
Экспозиция:

``` r
aspect = terrain(dem, 'aspect', 
                 unit = 'degrees')

ggplot() +
  geom_spatraster(data = aspect) +
  scale_fill_gradientn(
    colors = rainbow(9), 
    values = 0:8 / 8
  ) +
  labs(title = 'Экспозиция', 
       fill = 'Градусы [°]')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-14-1.png)
]

---

## Отмывка рельефа

Отмывка рельефа средствами terra получается при комбинировании растров углов наклона и экспозиции склона. Для повышения контрастности можно предварительно умножить высоты на определенный масштабный коэффициент:

``` r
# углы наклона для классической отмывки
slope_rad10 = terrain(dem * 10, 'slope', unit = 'radians')

# углы наклона для вертикальной отмывки
slope_rad5 = terrain(dem * 5, 'slope', unit = 'radians')

# экспозиция для любой отмывки
aspect_rad = terrain(dem, 'aspect', unit = 'radians')
```

---

## Отмывка рельефа

.pull-left[
Классическая отмывка при северо-западном освещении:

``` r
hill = shade(
  slope = slope_rad10, 
  aspect = aspect_rad, 
  angle = 45, 
  direction = 315
)

ggplot() +
  geom_spatraster(data = hill) +
  scale_fill_gradient(
    low = 'black', 
    high = 'white'
  ) +
  ggtitle('Отмывка классическая')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-16-1.png)
]

---

## Отмывка рельефа

.pull-left[
Отвесная отмывка при вертикальном освещении (подчеркивает линии хребтов):

``` r
hill_vert = shade(
  slope = slope_rad5, 
  aspect = aspect_rad, 
  angle = 90
)

ggplot() +
  geom_spatraster(data = hill_vert) +
  scale_fill_gradient(
    low = 'black', 
    high = 'white'
  ) +
  ggtitle('Отмывка отвесная')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-17-1.png)
]

---

## Отмывка рельефа

.pull-left[
Кобинированная отмывка:

``` r
hill_comb = hill * hill_vert

ggplot() +
  geom_spatraster(data = hill_comb) +
  scale_fill_gradient(
    low = 'black', 
    high = 'white'
  ) +
  ggtitle('Отмывка комбинированная')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-18-1.png)
]

---

## Гидрологический анализ ЦМР

**Гидрологический анализ ЦМР** --- совокупность методов обработки ЦМР, связанных с анализом и моделированием геометрических условий распределения поверхностного стока.

Распространенные задачи гидрологического анализа включают:

1.  Определение направлений стока.
2.  Вычисление площади водосбора.
3.  Построение сети потенциальных водотоков.
4.  Построение границ бассейнов.

> Для гидрологического анализа играют очень важную роль процедуры предварительной подготовки ЦМР, такие как *устранение ошибочных замкнутых локальных понижений*.

---

## Направление и аккумуляция стока

``` r
setwd('/Users/tsamsonov/GitHub/r-geo-course/data')
wbt_flow_accumulation_full_workflow(
  dem       = 'dem_fergana.tif',
  out_dem   = 'fill.tif', # ЦМР без замкнутых понижений
  out_pntr  = 'dir.tif',  # Направления стока по D8
  out_accum = 'acc.tif',  # Аккумуляция стока
  out_type = 'ca'         # Тип аккумуляции: площадь водосбора
)

fill = rast('fill.tif')
dir = rast('dir.tif')
acc = rast('acc.tif')
```

---

## Направление и аккумуляция стока

.pull-left[

``` r
levels(dir) = data.frame(
  value = c(0, 2^(0:7)), 
  desc = c('NA', 'NE', 'E', 'SE', 
           'S', 'SW', 'W', 'NW', 'N')
)

ggplot() +
  geom_spatraster(data = dir) +
  scale_fill_manual(values = rev(rainbow(9))) +
  labs(title = 'Направление стока', 
       fill = '')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-20-1.png)
]

---

## Направление и аккумуляция стока

.pull-left[

``` r
ggplot() +
  geom_spatraster(data = acc / 10^9) +
  scale_fill_gradient(low = 'white', high = 'black', trans = 'log') +
  labs(title = 'Аккумуляция стока (площадь)', 
       fill = expression(Тыс.~км^2))
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-21-1.png)
]

---

## Потенциальные водотоки

Выделяются по растрам направления и аккумуляции стока

``` r
setwd('/Users/tsamsonov/GitHub/r-geo-course/data')
wbt_extract_streams(
  flow_accum = 'acc.tif',
  output = 'str.tif',
  threshold = 1e8 ## 100 кв.км.
)
wbt_raster_streams_to_vector(
  streams = 'str.tif',
  d8_pntr = 'dir.tif',
  output = 'streams.shp',
)
streams = read_sf('streams.shp') |> 
  st_set_crs(crs(dem, proj=TRUE))
```

---

## Потенциальные водотоки

.pull-left[

``` r
ggplot() +
  geom_spatraster(data = dem) +
  geom_sf(data = streams, 
          color = 'deepskyblue4') +
  scale_fill_hypso_c() +
  labs(title = expression(Потенциальные~водотоки~(A~`>`~100~км^2)), 
       fill = 'Высота, [м]')
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-23-1.png)
]

---

## Водосборные бассейны

Генрируются для указанных точек по растру направлений стока

``` r
setwd('/Users/tsamsonov/GitHub/r-geo-course/data')
wbt_watershed(
 d8_pntr = 'dir.tif',
 pour_pts = 'pour_pts.shp',
 output = 'wsh.tif'
)
wbt_raster_to_vector_polygons(
 input = 'wsh.tif', 
 output = 'watersheds.shp'
)

watersheds = read_sf('watersheds.shp') |> 
  st_set_crs(crs(dem, proj=TRUE))

pour_pts = read_sf('pour_pts.shp') |> 
  st_set_crs(crs(dem, proj=TRUE))
```

---

## Водосборные бассейны

.pull-left[

``` r
ggplot() +
  geom_spatraster(data = dem) +
  scale_fill_hypso_c() +
  labs(title = 'Бассейны', 
       fill = 'Высота, [м]') +
  ggnewscale::new_scale_fill() +
  geom_sf(data = watersheds, 
          mapping = aes(fill = factor(FID)), 
          linewidth = 1, color = 'black', 
          alpha = 0.5, show.legend = FALSE) +
  geom_sf(data = streams, 
          color = 'deepskyblue4') +
  geom_sf(data = pour_pts, size = 2, 
          fill = 'white', shape = 21,
          stroke = 1.5)
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-25-1.png)
]

---

## Анализ расстояний

Доступность станций метро

``` r
roads = read_sf("../r-geo-course/data/roads.gpkg") # Дороги
poi = read_sf("../r-geo-course/data/poi_point.gpkg") # Точки интереса
rayons = read_sf("../r-geo-course/data/boundary_polygon.gpkg") # Границы районов
stations = read_sf("../r-geo-course/data/metro_stations.gpkg") # Станции метро
water = read_sf("../r-geo-course/data/water_polygon.gpkg") # Водные объекты

dist_grid = stations |> 
  rast(resolution = 25) |> 
  rasterize(stations, y = _) |> 
  gridDist(NA)
```

---

## Анализ расстояний

.pull-left[
.code-small[

``` r
ggplot() +
  geom_spatraster_contour_filled(
    data = dist_grid
  ) +
  geom_spatraster_contour(
    data = dist_grid, 
    linewidth = 0.25, 
    color = 'black') +
  geom_sf(data = water, 
          linewidth = 0.1) +
  geom_sf(data = roads, 
          linewidth = 0.1) +
  geom_sf(data = stations, 
          color = 'darkviolet') +
  labs(title = 'Расстояние до ближайшей станции метро',
       fill = 'м') +
  theme_bw()
```
]
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-27-1.png)
]

---

## Зональные операции

``` r
files = list.files('/Volumes/Data/Spatial/WorldClim/2.1/tavg',
                   full.names = TRUE)
temp = rast(files)

data(land, package = 'tmap')
terraland = rast(land['cover'])

pal = c("#003200", "#3C9600", "#006E00", "#556E19", "#00C800", 
        "#8CBE8C", "#467864", "#B4E664", "#9BC832", "#EBFF64", 
        "#F06432", "#9132E6", "#E664E6", "#9B82E6", "#B4FEF0", 
        "#646464", "#C8C8C8", "#FF0000", "#FFFFFF", "#5ADCDC")
ggplot() +
  geom_spatraster(data = terraland) +
  scale_fill_manual(values = pal, 
                    guide = guide_legend(ncol = 3), 
                    name = NULL) + 
  theme(legend.position = 'bottom')
```

---

## Зональные операции

---

## Зональные операции

Необходимо привести растры к общей сетке

``` r
(cover = project(terraland, temp, method = 'near'))
## class       : SpatRaster 
## dimensions  : 1080, 2160, 1  (nrow, ncol, nlyr)
## resolution  : 0.1666667, 0.1666667  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## categories  : label 
## name        :                      cover 
## min value   : Broadleaf Evergreen Forest 
## max value   :               Water bodies
```

---

## Зональные операции

.pull-left[
.code-small[

``` r
cover_north = crop(
  cover, ext(-180, 180, 0, 90)
)
temp_north = crop(
  temp, ext(-180, 180, 0, 90)
)
stats = zonal(temp_north, cover_north, 
              mean, na.rm = TRUE)

zonal_stats = stats |> 
  rename(cover = 1) |> 
  pivot_longer(
    -cover, 
    names_to = 'month', 
    values_to = 'tavg', 
    names_prefix = 'wc2.1_10m_tavg_', 
    names_transform = as.integer
  ) |> 
  mutate(month = ordered(month, 1:12, 
                         month.abb))
```
]
]

.pull-right[
.code-vsmall[

``` r
head(zonal_stats, 20)
## # A tibble: 20 × 3
##    cover                      month   tavg
##    <chr>                      <ord>  <dbl>
##  1 Broadleaf Evergreen Forest Jan   22.5  
##  2 Broadleaf Evergreen Forest Feb   23.1  
##  3 Broadleaf Evergreen Forest Mar   24.0  
##  4 Broadleaf Evergreen Forest Apr   24.7  
##  5 Broadleaf Evergreen Forest May   24.9  
##  6 Broadleaf Evergreen Forest Jun   24.7  
##  7 Broadleaf Evergreen Forest Jul   24.5  
##  8 Broadleaf Evergreen Forest Aug   24.6  
##  9 Broadleaf Evergreen Forest Sep   24.6  
## 10 Broadleaf Evergreen Forest Oct   24.3  
## 11 Broadleaf Evergreen Forest Nov   23.5  
## 12 Broadleaf Evergreen Forest Dec   22.7  
## 13 Broadleaf Deciduous Forest Jan   -7.22 
## 14 Broadleaf Deciduous Forest Feb   -5.11 
## 15 Broadleaf Deciduous Forest Mar    0.246
## 16 Broadleaf Deciduous Forest Apr    7.19 
## 17 Broadleaf Deciduous Forest May   13.1  
## 18 Broadleaf Deciduous Forest Jun   17.6  
## 19 Broadleaf Deciduous Forest Jul   19.6  
## 20 Broadleaf Deciduous Forest Aug   18.4
```
]
]

---

## Зональные операции

.left-40[
Климатограммы:

``` r
ggplot(zonal_stats) +
  geom_line(
    aes(x = month, 
        y = tavg, 
        color = cover, 
        group = cover), 
    size = 1
  ) +
  scale_color_manual(
    values = pal
  )
```
]

.right-60[
<img src="14_Raster_files/figure-html/unnamed-chunk-32-1.png" height="500px" />
]

---

## Глобальные операции

.pull-left[

``` r
global(temp, max, na.rm = T) 
##                        max
## wc2.1_10m_tavg_01 34.00950
## wc2.1_10m_tavg_02 32.82425
## wc2.1_10m_tavg_03 32.90950
## wc2.1_10m_tavg_04 34.19375
## wc2.1_10m_tavg_05 36.25325
## wc2.1_10m_tavg_06 38.35550
## wc2.1_10m_tavg_07 39.54950
## wc2.1_10m_tavg_08 38.43275
## wc2.1_10m_tavg_09 35.79000
## wc2.1_10m_tavg_10 32.65125
## wc2.1_10m_tavg_11 32.78800
## wc2.1_10m_tavg_12 32.82525
```
]

.pull-right[

``` r
global(temp, min, na.rm = T) 
##                         min
## wc2.1_10m_tavg_01 -45.88400
## wc2.1_10m_tavg_02 -44.80000
## wc2.1_10m_tavg_03 -57.92575
## wc2.1_10m_tavg_04 -64.19250
## wc2.1_10m_tavg_05 -64.81150
## wc2.1_10m_tavg_06 -64.35825
## wc2.1_10m_tavg_07 -68.46075
## wc2.1_10m_tavg_08 -66.52250
## wc2.1_10m_tavg_09 -64.56325
## wc2.1_10m_tavg_10 -55.90000
## wc2.1_10m_tavg_11 -43.43475
## wc2.1_10m_tavg_12 -45.32700
```
]

---

## Извлечение данных

Данные можно извлекать посредством `extract` по объектам различной размерности. Например, в точке:

.pull-left[
.code-vsmall[

``` r
pnt = tibble::tibble(
  x = -45, y = 70,
  z = terra::extract(ice.depth, 
                     bind_cols(x,y))[, 2]
)

ggplot() +
  geom_spatraster(data = ice.depth) +
  scale_fill_gradient(low = 'white', 
                      high = 'navyblue', 
                      na.value = "transparent") +
  geom_point(aes(x, y), pnt, 
             size = 3, color = 'red') +
  geom_text(aes(x, y, label = z), 
            pnt, fontface = 'bold', size = 5,
            vjust = 'bottom', hjust = 'left',
            nudge_x = 1, nudge_y = 1) +
  coord_sf(xlim = c(-80, -10), ylim = c(58, 85)) +
  labs(title = 'Мощность покровного оледенения',
       fill = '[м]', x = NULL, y = NULL)
```
]
]

.pull-right[
<img src="14_Raster_files/figure-html/unnamed-chunk-35-1.png" height="450px" />
]

---

## Извлечение данных

Одна из наиболее распространенных задач по извлечению растровых данных — это построение профиля вдоль заданной линии. Воспользуемся интерактивным редактором для проведения линии профиля

``` r
mp = mapview(temp$tmean6)
profile = mapedit::drawFeatures(mp)
```

---

## Извлечение данных

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

``` r
temprof = terra::extract(temp, profile, cells = TRUE, xy = TRUE)

tempdf = temprof %>%
  mutate(dist = 0.001 * c(0, dplyr::select(., x, y) |>  
                             geosphere::distGeo() |> 
                             cumsum() |> 
                             head(-1))) |>  
  select(-c(x, y, cell, ID)) |> 
  pivot_longer(-dist, 
               names_to = 'month', values_to = 'tavg', 
               names_prefix = 'wc2.1_10m_tavg_', 
               names_transform = as.integer)
```

---

## Извлечение данных

.pull-left[

``` r
ggplot(tempdf, 
       aes(x = dist, y = tavg)) +
  geom_line() +
  geom_smooth(span = 0.1) +
  annotate('text', x = 0, 
           y = 10, label = 'A') +
  annotate('text', x = max(tempdf$dist), 
           y = 10, label = 'B') +
  ggtitle('Профиль среднемесячной 
          температуры июня по линии A—B') +
  facet_wrap(~month)
```
]

.pull-right[
![](14_Raster_files/figure-html/unnamed-chunk-40-1.png)
]