Глава 14 Интерполяция геополей

Программный код главы

В лекции рассмотрены детерминистические методы восстановления непрерывных поверхностей по данным в точках. Детерминистические методы интерполяции производят интерполяцию значений на основе заданной аналитической зависимости между значением в точке и значениями в пунктах с исходными данными. Это отличает их от методов геостатистических, где эта зависимость находится статистическим путем. Геостатистические методы будут рассмотрены далее.

14.1 Введение

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

Методы которые производят интерполяцию на основе заданной аналитической зависимости, называют детерминистическими. Параметры этой зависимости могут быть как априори заданы пользователем, так и определяться автоматически одним из методов оптимизации — в частности, по методу наименьших квадратов.

Например, мы можем сказать, что между соседними точками показатель меняется линейным образом (здесь нужно еще указать, что мы понимаем под соседством). Такие методы достаточно просты в использовании и интерпретации. В то же время, они не учитывают статистических особенностей поведения величины между точками, которое определяется ее автокорреляционными свойствами. Методы интерполяции, которые учитывают пространственную автокорреляцию, называют геостатистическими. Они более сложны в использовании, но потенциально могут дать более достоверные результаты.

В этом модуле мы познакомимся со следующими детерминистическими методами интерполяции.

  • Метод ближайшего соседа (nearest neighbor)
  • Метод естественного соседа (nearest neighbor)
  • Метод интерполяции на основе триангуляции
  • Метод обратно взвешенных расстояний (ОВР)
  • Метод радиальных базисных функций (РБФ)
  • Метод иерархических базисных сплайнов (ИБС)

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

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

library(sf)
library(stars)
library(tmap)
library(raster)
library(plotly)
library(mapview)
library(tidyverse)
library(ggrepel)
library(dismo) # библиотека species distribution modelling
library(akima) # библиотека для интерполяции на основе триангуляции
library(gstat) # библиотека для геостатистической интерполяции, построения трендов и IDW
library(deldir) # библиотека для построения триангуляции Делоне и диаграммы Вороного
library(fields) # радиальные базисные функции (сплайны минимальной кривизны)
library(MBA) # иерархические базисные сплайны

# Убираем экспоненциальное представление больших чисел
options(scipen=999)

# Читаем слои картографической основы
cities = st_read("data/Italy_Cities.gpkg") %>% # Города
  bind_cols(st_coordinates(.) %>% as_tibble())
## Reading layer `Italy_Cities' from data source 
##   `/Users/tsamsonov/GitHub/r-geo-course/data/Italy_Cities.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 8 features and 37 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 368910.4 ymin: 4930119 xmax: 686026 ymax: 5115936
## Projected CRS: WGS 84 / UTM zone 32N
rivers = st_read("data/Italy_Rivers.gpkg") # Реки
## Reading layer `Italy_Rivers' from data source 
##   `/Users/tsamsonov/GitHub/r-geo-course/data/Italy_Rivers.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 23 features and 10 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 332178.5 ymin: 4922880 xmax: 758033.3 ymax: 5121416
## Projected CRS: WGS 84 / UTM zone 32N
lakes = st_read("data/Italy_Lakes.gpkg")   # Озера
## Reading layer `Italy_Lakes' from data source 
##   `/Users/tsamsonov/GitHub/r-geo-course/data/Italy_Lakes.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 6 features and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 460686.2 ymin: 4938750 xmax: 757443.7 ymax: 5113557
## Projected CRS: WGS 84 / UTM zone 32N

# Читаем ЦМР — цифровую модель рельефа на регулярной сетке
dem = read_stars("data/gtopo.tif")

# Читаем данные об осадках
pts = read_table2("data/Rainfall.dat") %>% 
  st_as_sf(coords = c('x', 'y'), 
           crs = st_crs(cities),
           remove = FALSE)

# Координаты пригодятся нам в дальнейшем
coords = st_coordinates(pts)

Построение карты

# Цветовая шкала для рельефа
dem_colors = colorRampPalette(c("darkolivegreen4", "lightyellow", "orange", "firebrick"))

# Шкала высот для рельефа
dem_levels = c(0, 50, 100, 200, 500, 1000, 1500, 
               2000, 2500, 3000, 3500, 4000, 5000)

dem_ncolors = length(dem_levels) - 1

dem_contours = st_contour(dem, breaks = dem_levels, contour_lines = TRUE)

old = theme_set(theme_bw())

# Смотрим как выглядит результат
ggplot() +
  geom_stars(data = cut(dem, breaks = dem_levels)) +
  coord_sf(crs = st_crs(dem)) +
  scale_fill_manual(name = 'м',
                    values = dem_colors(dem_ncolors),
                    labels = dem_levels,
                    na.translate = FALSE,
                    guide = guide_legend(label.vjust = -0.3, reverse = TRUE)) +
  geom_sf(data = dem_contours, color = 'black', size = 0.2) +
  geom_sf(data = rivers, color = 'midnightblue', size = 0.2) +
  geom_sf(data = lakes, color = 'midnightblue', fill = 'lightblue', size = 0.2) +
  geom_sf(data = pts, color = 'black', size = 0.5) +
  geom_sf(data = cities, shape = 21, fill = 'white') + 
  geom_text_repel(data = cities, mapping = aes(x = X, y = Y, label = name), box.padding = 0.2)

14.2 Построение сетки

Любопытным свойством пакетов R, отвечающих за интерполяцию данных, является их индифферентность относительно того, в каких точках эта интерполяция будет производиться. Это может быть как регулярная растровая сетка, так и множество точек в совершенно произвольных конфигурациях. Подобная гибкость делает процесс интерполяции данных чуть более сложным, чем в ГИС-пакетах, однако способствует полному и глубокому пониманию происходящего. Вы своими руками задаете пункты, в которых следует интерполировать значения.

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

# ПОСТРОЕНИЕ СЕТКИ ДЛЯ ИНТЕРПОЛЯЦИИ
# получим ограничивающий прямоугольник вокруг точек:
box = st_bbox(pts)
envelope = box[c(1,3,2,4)]

px_grid = st_as_stars(box, dx = 10000, dy = 10000)

ggplot() + 
  geom_sf(data = pts, color = 'red') +
  geom_sf(data = st_as_sf(px_grid), size = 0.5, fill = NA)

Получившееся представление можно назвать сеточной моделью. В пределах каждой ячейки величина будет считаться постоянной. Ее значение будет интерполировано в центре пиксела. Можно видеть, что интерполяционная сетка слегка выходит за пределы исходного охвата множества точек. Это связано с тем, что размеры прямоугольника, ограничивающего множество точек, не кратны выбранному разрешению растра (\(10 000\) м). Проблема, однако, не столько критична, если вы выбираете достаточно подробное (малое) разрешение растра, что мы и сделаем. Зададим его равным \(1 000\) м:

# создадим детальную растровую сетку
px_grid = st_as_stars(box, dx = 1000, dy = 1000)

# извлечем координаты точек в соответствующие столбцы, они нам пригодятся:
coords_grid = st_coordinates(px_grid)

14.3 Интерполяционные методы

# Цветовая шкала для осадков
rain_colors = colorRampPalette(c("white", "dodgerblue", "dodgerblue4"))

# Шкала количества осадков и соответствющее число цветов
rain_levels = seq(0, 80, by=10)
rain_ncolors = length(rain_levels)-1

rain_legend = scale_fill_manual(name = 'мм',
                                values = rain_colors(rain_ncolors),
                                guide = guide_legend(label.vjust = -0.3, reverse = TRUE, title.position = "bottom"),
                                labels = rain_levels,
                                na.value = 'white',
                                drop = FALSE)
rain_mapping = aes(fill = cut(rain_24, breaks = rain_levels))

14.3.1 Метод ближайшего соседа (nearest neighbour)

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

# МЕТОД БЛИЖАЙШЕГО СОСЕДА (NEAREST NEIGHBOR)

# Диаграмма Вороного
voronoi_sf = voronoi(coords, envelope) %>% 
  st_as_sf() %>% 
  st_set_crs(st_crs(pts)) %>% 
  st_join(pts)

# Триангуляция Делоне
edges = pts %>% 
  st_union() %>% 
  st_triangulate()

# Визуализация
ggplot() +
  geom_sf(data = voronoi_sf,
          mapping = rain_mapping,
          size = 0.2) +
  rain_legend +
  geom_sf(data = pts, color = 'red', size = 0.5) +
  geom_sf(data = edges, color = 'red', size = 0.1, fill = NA)

Если есть задача конвертировать это в растр, то надо сформировать новый растр и “перенести” на него информацию с полигонов:

# Создаем растр
rnn = st_rasterize(voronoi_sf["rain_24"], px_grid)

# Визуализируем:
ggplot() +
  geom_stars(data = cut(rnn, breaks = rain_levels)) +
  rain_legend +
  coord_sf(crs = st_crs(pts))

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

  • x - вектор координат ячеек по оси \(Х\)
  • y - вектор координат ячеек по оси \(Y\)
  • z - матрицу значений, имеющую размеры \(length(x) \times length(y)\)

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

rain_colors3d = colorRamp(c("white", "dodgerblue", "dodgerblue4"))

x = coords_grid[,'x'] %>% unique() # Получим координаты столбцов
y = coords_grid[,'y'] %>% unique() # Получим координаты строк

p = plot_ly(x = x, 
            y = y, 
            z = rnn$rain_24, 
            type = "surface",
            colors = rain_colors3d)
layout(p, scene = list(aspectratio = list(x = 1, y = 1, z = 0.3)
))