09-SpatialData.qmd

# Пространственные данные {#spatial}

```{r setup-spatial, echo = FALSE, purl = FALSE, include=FALSE}
knitr::opts_knit$set(global.par = TRUE)
knitr::opts_chunk$set(warning = FALSE, message = FALSE, collapse = TRUE, out.width = '100%')
knitr::knit_hooks$set(crop = knitr::hook_pdfcrop)
```

## Предварительные требования {#spatial_prerequisites}

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

Необходимые для работы пакеты:

```{r}
library(sf)
library(stars)
library(mapview)
library(dplyr)
library(readr)
```

## Модели пространственных данных {#spatial_models}

**Пространственные данные** (spatial data) --- это данные о пространственных объектах и их наборах. В свою очередь, пространственный объект определяется как *цифровая модель материального или абстрактного объекта реального или виртуального мира с указанием его идентификатора, координатных и атрибутивных данных* [^09-spatialdata-1].

[^09-spatialdata-1]: ГОСТ Р 52438-2005 \<\<Географические информационные системы. Термины и определения\>\>. В стандарте поясняется, что объектом может быть неподвижный или движущийся простой или сложный объект, явление, событие, процесс и ситуация. Моделируемый объект может относиться к территории, акватории, недрам и воздушному пространству Земли, околоземному космическому пространству, другим космическим телам и небесной сфере. В широком смысле под пространственным объектом в геоинформатике понимается как сам объект, так и адекватная ему цифровая модель

Если говорить по сути, то пространственные данные можно определить как *данные о географических объектах или явлениях, фиксирующие их местоположение и/или распределение в системе координат, привязанной к телу Земли или любого другого небесного тела*. Таким образом, отличительной особенностью пространственных данных перед непространственными является координатное описание местоположения. На профессиональном жаргоне пространственные данные также часто называют **геоданными**. Следует помнить, что этот термин не является научным, и его не следует использовать в публикациях и квалификационных работах.

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

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

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

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

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

## Векторные данные {#vector_data_r}

### Simple Features {#simple_features}

**Simple Features** (официально *Simple Features Access*) --- это стандарт [OGC 06-103](http://www.opengeospatial.org/standards/sfa), разработанный Open Geospatial Consortium (OGC) и реализованный также в виде международного стандарта [ISO 19125](https://www.iso.org/standard/40114.html), который определяет общую модель хранения и доступа к векторным объектам (точка, линия, многоугольник, мульти точечные, мультилинии и т. д.), в географических информационных системах.

Геометрическое представление пространственных объектов базируется на следующих принципах:

-   Все геометрии состоят из точек.
-   Точки являются координатами в 2-, 3- или 4-мерном пространстве.
-   Все точки в геометрии имеют одинаковую размерность.

В дополнение к координатам $X$ и $Y$ имеются два дополнительных дополнительных параметра:

-   координата $Z$, обозначающая высоту
-   координата $M$, обозначающая некоторую меру, связанную с точкой, а не с признаком в целом (в этом случае это будет атрибут объекта). Измерение $M$ может быть использовано, например, для представления времени или линейных координат (для маршрутов).

Координаты простой геометрии всегда содержат компоненты $X$ и $Y$, поэтому все разнообразие возможных представлений определяется наличием или отсутствием дополнительных измерений $Z$ и $M$. Таким образом, получаем **четыре** варианта геометрии:

-   двумерные точки $XY$
-   трехмерные точки $XYZ$
-   трехмерные точки $XYM$
-   четырехмерные точки $XYZM$

В случае использования широт и долгот $X$ соответствует долготе, $Y$ соответствует широте.

Всего стандарт **Simple Features** включает в себя 17 типов геометрий. Из них наиболее употребительными являются следующие 7:

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

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

```{r, echo = FALSE}
p = st_point(c(0.5,0.5))

pc = rbind(c(0.5,0.5), c(1, 3), c(2, 1), c(0.2, 2), c(2, 3), c(1.5, 1.5))
mp = st_multipoint(pc)

s1 = rbind(c(0,1),c(0.5,1.5),c(1.2,1.2),c(2,1.3),c(3,2))
ls = st_linestring(s1)

s1 = rbind(c(0.5,1.5),c(1.2,1.2),c(2,1.3))
s2 = rbind(c(0,1.5),c(0.5,2.0),c(1.2,1.7))
s3 = rbind(c(2,1.8),c(3,2.5))
mls = st_multilinestring(list(s1,s2,s3))

p1 = rbind(c(0.5,0.5), c(2,0), c(3,2), c(1.5,4), c(0,3), c(0.5,0.5))
p2 = rbind(c(1,1), c(0.8,2), c(2,2.2), c(1.4,1.1), c(1,1))
pol = st_polygon(list(p1, p2))
p3 = rbind(c(3,3.3), c(3.5, 3.1), c(4,3), c(4,3.7), c(3.7, 3.96), c(3.2,4), c(3,3.3))
p4 = rbind(c(3.2,3.4), c(3.8,3.2), c(3.8,3.7), c(3.3,3.8), c(3.2,3.4))
p5 = rbind(c(3,1.2), c(2.5,0.2), c(3.5,0.2), c(3.5,1.2), c(3,1.2))
p6 = rbind(c(0,1), c(0.1,0.8), c(0.2,0.5), c(0.1,0.3), c(0, 0.7), c(0,1))
mpol = st_multipolygon(list(list(p1,p2), list(p3,p4), list(p5), list(p6)))

gc = st_geometrycollection(list(pol, mp + c(3, -1), mls + c(2.5,2)))
```

```{r, echo=FALSE,fig.width=16, fig.height=8, fig.show='hold'}
par(mfrow = c(2,3))
plot(p, cex = 5, pch = 20, main = 'POINT', cex.main=3)
plot(ls, lwd = 3, main = 'LINESTRING', cex.main=3)
plot(pol, lwd = 3, col = 'grey', main = 'POLYGON', cex.main=3)
plot(mp, cex = 5, pch = 20, main = 'MULTIPOINT', cex.main=3)
plot(mls, lwd = 3, main = 'MULTILINESTRING', cex.main=3)
plot(mpol, lwd = 3, col = 'grey', main = 'MULTIPOLYGON', cex.main=3)
# plot(gc, lwd = 3, col = 'grey', border = 'black', main = 'GEOMETRYCOLLECTION', cex.main=3)
par(mfrow = c(1,1))
```

Оставшиеся виды геометрий *Simple Features* включают: `CIRCULARSTRING`, `COMPOUNDCURVE`, `CURVEPOLYGON`, `MULTICURVE`, `MULTISURFACE`, `CURVE`, `SURFACE`, `POLYHEDRALSURFACE`, `TIN`, `TRIANGLE`.

Существует два официально закрепленных формата представления SF: *Well-Known Text (WKT)* и *Well-Known Binary (WKB)*, которые необходимы для чтения таких данных человеком и машиной соответственно.

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

```{r, echo=FALSE, collapse=T}
cat(st_as_text(p))
cat(st_as_text(ls))
cat(st_as_text(pol))
cat(st_as_text(mp))
cat(st_as_text(mls))
cat(st_as_text(mpol))
cat(st_as_text(gc))
```

**Well-Known Binary (WKB)** --- бинарный формат хранения координат. Именно этот формат фактически используется в базах данных, поскольку он обеспечивает высокую скорость чтения и записи данных (в отличие от текстового). Однако внешний вид данных в формате WKB мало о чем говорит человеку, поскольку он предназначен для чтения компьютером. Например, вышеприведенная строка `LINESTRING` будет выглядеть так:

```{r, echo=FALSE}
cat(st_as_binary(ls))
```

### Базовые библиотеки {#vector_data_packages}

В R существует высоко развитая инфраструктура для работы с векторными данными, которая обеспечивается пакетом [**sf**](https://cran.r-project.org/web/packages/sf/index.html).

Пакет **sf** базируется на библиотеках [PROJ](https://proj.org), [GDAL](https://gdal.org), [GEOS](https://trac.osgeo.org/geos/) и [S2](https://s2geometry.io), которые устанавливаются вместе с ним. Их назначение кратко описано на следующем рисунке:

```{r, echo=FALSE, fig.cap='Архитектура программных библиотек для работы с пространственными данными в R'}
knitr::include_graphics('images/sf_architecture_new.svg')
```

Со многими функциями **sf** мы познакомимся в последующих разделах нашего курса. Некоторые из них (такие как `arrange`, `filter`, `mutate` из пакета **dplyr**), должны быть уже знакомы вам по предыдущим лекциям. Можно обратить внимание на то, что практически все функции начинаются с префикса `st_`, что означает **"spatiotemporal"**. Данные префиксы были выбраны для унификации с аналогичными названиями функций, используемых в широко распространенной СУБД PostgreSQL для оперирования объектами *Simple Features*.

### Чтение {#vector_data_reading}

Существует большое количество форматов хранения пространственных данных. Но в общем и целом их можно разделить на две категории: файловые форматы (наиболее привычные пользователям) и хранение данных в СУБД --- системах управления базами данных. Благодаря библиотеке GDAL пакет **sf** имеет возможность читать и записывать [более 90 различных форматов векторных даных](http://www.gdal.org/ogr_formats.html).

Исторически наиболее распространенным форматом был (и остается) [ESRI Shapefile](https://www.esri.com/library/whitepapers/pdfs/shapefile.pdf). Данный формат, однако не отвечает современным техническим требованиям с точки зрения гибкости, соответствия стандартам и возможностям хранения разнообразных типов геометрий (напомним, что в стандарте **Simple Features** их 17, а с учетом четырех вариантов размерности точек получается целых 68 ). Современный формат, который обеспечивает полную поддержку стандарта **Simple Features** (и не только) --- это [GeoPackage](http://www.geopackage.org). Именно его мы и будем использовать в нашем практикуме.

Для чтения данных средствами sf необходимо использовать функцию `st_read()`:

```{r, collapse = TRUE}
countries = st_read('data/ne/countries.gpkg')
```

Лог функции сообщил нам следующую информацию:

-   Набор данных представляет собой коллекцию из 183 пространственных объектов с 72 атрибутами
-   Тип геометрии `MULTIPOLYGON`
-   Размерность геометрии $XY$
-   Ограничивающий прямоугольник (разброс координат) по осям $X$ и $Y$ имеет диапазон $[-180, 180] \times [-90, 83.64513]$
-   Проекция (CRS --- coordinate reference system) имеет название *WGS 84*.

Подгрузим также для работы данные по другим типам объектов:

```{r}
oceans = st_read('data/ne/oceans.gpkg')
rivers = st_read('data/ne/rivers.gpkg')
lakes = st_read('data/ne/lakes.gpkg')
cities = st_read('data/ne/cities.gpkg')
```

### Внутренняя структура {#sf_structure}

Традиционно во всех ГИС-приложениях и базах пространственных данных множество пространственных объектов представляется в виде таблицы атрибутов, где каждая строка соответствует объекту, а каждый столбец --- атрибуту объекта. С каждой строкой таблицы должна быть ассоциирована информация о геометрии объекта, которая, в зависимости от формата данных, может либо храниться непосредственно в таблице (в специальном столбце), либо быть вынесена в отдельную структуру данных, которая связана с таблицей атрибутов посредством ключа[^09-spatialdata-2].

[^09-spatialdata-2]: Например, в широко распространенном формате **Esri Shapefile** атрибутивная таблица хранится в файле `*.dbf` формата *DBASE*, геометрия хранится в отдельном файле `*.shp`, а связь между ними осуществляется через файл `*.shx`. Разбиение формата хранения на несколько файлов --- это одна из уязвимостей шейп-файлов: при отсутствии хотя бы одного из этих файлов данные прочесть стандартными средствами (без дополнительного хакинга) будет нельзя.

В **R** используется первый подход, в котором информация о геометрии хранится в специальном столбце таблицы. Каждая ячейка этого столбца соответствует геометрическому объекту *Simple Features*. Представление геометрических объектов реализовано стандартными средствами, такими как списки, матрицы и векторы. Эти структуры данных упорядоченным образом хранят координаты объектов и естественным образом соответствуют способу организации данных, который регламентируется стандартом *Simple Features*. Поскольку геометрический столбец хранит не обычные переменные, а структуры данных, он реализуется в виде так называемого *списка-колонки (list-column)*, каждый элемент которой соответствует отдельному объекту.

Исходя из этих соображений, представление пространственных объектов реализовано в **R** в виде иерархии из трех классов объектов:

1.  `sf` (simple features) --- объект класса `data.frame`, представляющий множество пространственных объектов со списком-колонкой для хранения геометрии
2.  `sfc` (simple features geometry column) --- список-колонка в объекте `sf`, представляющий множество геометрий пространственных объектов
3.  `sfg` (simple feature geometry) --- геометрия пространственного объекта внутри списка `sfc`

В соответствии с перечисленными спецификациями происходит работа с пространственными объектами. То что, объекты типа Simple Features реализованы в виде самых обычных фреймов данных, означает что *любая операция, применимая к фрейму данных, будет также применима к объекту типа* `sf`. Это очень важная особенность объектов типа sf, которой сильно не хватало в экосистеме исторического пакета `sp`.

Посмотрим, как все это реализовано, на конкретном примере:

```{r, collapse = TRUE}
class(countries)
```

Данная форма записи говорит о том, что прочитанный слой имеет класс *sf*, который, в свою очередь, является расширением класса *data.frame*.

А теперь посмотрим на последние колонки в первых строках таблицы:

```{r}
head(countries[tail(colnames(countries))])
```

Видно, что геометрия пространственных объектов хранится в заключительном столбце с названием `geometry`. Данный столбец можно быстро извлечь, применив функцию `st_geometry()`. Полученный объект будет иметь тип **sfc** (Simple Feature Geometry Column)

```{r, collapse = TRUE}
outlines = st_geometry(countries)
class(outlines)
```

Полученный вывод говорит нам о том, что наши объекты имеют класс `sfc_MULTIPOLYGON`, который является расширением класса `sfc` (simple feature geometry column).

Теперь если просмотреть начало данных, то мы увидим, что это больше не фрейм данных, а аннотированный список:

```{r, collapse = TRUE}
head(outlines)
```

Далее можно опуститься на базовый уровень геометрии, получив доступ к отдельному объекту. Поскольку объект класса `sfc` представляет собой список, любой элемент можно извлечь по его порядковому номеру. Класс полученного объекта будет:

```{r, collapse = TRUE}
class(outlines[[8]])
```

Исходя из полученной информации можно сделать вывод, что геометрия 8-го объекта таблицы `countries` имеет класс `sfg`, реализованный в виде мультиполигонов (`MULTIPOLYGON`) с плоскими координатами (`XY`)

Наконец, чтобы добраться до координат в чистом виде, необходимо развернуть иерархию списков, из которых состоит объект `sfg`. Количество уровней вложенности всегда зависит от конкретного объекта, их может быть достаточно много, особенно если объекты представлены мультиполигонами (несколько компонент связности), каждый из которых также состоит из полигонов с дырками. В нашем случае все достаточно просто, так как в слое `countries` дырок в полигонах нет, а 8-й по счету полигон состоит из одной-единственной геометрии, координаты которой в виде матрицы можно извлечь как:

```{r, collapse = TRUE}
outlines[[8]][[1]]
```

### Визуализация {#sf_plotting}

#### Базовая графическая система {#sf_plotting_basic}

Если попытаться применить функцию `plot()` к геометрии объекта, она попытается нарисовать тематические карты по всем имеющимся атрибутам (но остановится, если их более 9):

```{r, collapse = TRUE, error=FALSE, warning=FALSE}
plot(countries)
```

Если задача стоит нарисовать границы объектов, то нужно отображать объект **sfc**:

```{r, collapse = TRUE, error=TRUE}
plot(outlines, col = 'red')
```

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

```{r, collapse = TRUE}
plot(countries['sovereignt'], key.pos = NULL) # Здесь легенда не нужна
```

Для отображения координатной сетки надо указать параметр `graticule = TRUE`, а подписей координат --- `axes = TRUE`:

```{r}
plot(countries['gdp_md_est'], graticule = TRUE, axes = TRUE)
```

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

```{r}
cities_large = cities |> 
  filter(scalerank == 0, 
         ! name %in% c('Washington, D.C.', 'Paris', 'Riyadh', 'Rome', 'São Paulo', 'Kolkata'))

plot(st_geometry(countries), lwd = 0.5, border = 'gray')
plot(oceans, col = 'steelblue1', border = 'steelblue', add = TRUE)
plot(lakes, col = 'steelblue1', border = 'steelblue', add = TRUE)
plot(rivers, col = 'steelblue', add = TRUE)
plot(cities_large, col = 'black', pch = 19, cex = 0.25, add = TRUE)
text(cities_large$longitude, cities_large$latitude, 
     label = cities_large$name, cex = 0.5, pos = 2, offset = 0.25)
```

> **Внимание**: чтобы слои совместились на карте, они должна иметь одинаковую систему координат.

Ясно, что на полученных нами картах можно много что улучшить, однако это мы отложим до следующей главы, где подробно разбирается построение тематических карт в **R**.

> **Внимание**: чтобы слои данных можно было совместно анализировать и наносить на одну карту, они должны иметь одну и ту же координатную систему (проекцию).

#### Интерактивные карты {#spatial_interactive}

R предоставляет возможности для интерактивного просмотра пространственных данных средствами библиотек веб-картографирования. В данном разделе мы кратко познакомимся с возможностями пакета [**mapview**](https://r-spatial.github.io/mapview/), который использует возможности библиотеки [Leaflet](https://leafletjs.com/). Функции данного пакета не предназначены для создания тематических карт высокого качества и рассчитаны на выполнение исследовательского анализа данных.

Чтобы отобразить векторный или растровый слой средствами mapview, достаточно вызвать одноименную функцию данного пакета:

```{r, eval = FALSE}
mapview(countries)
```

```{r, echo = FALSE}
knitr::include_graphics('images/mapview1.png')
```

Чтобы отобразить определенный показатель, можно использовать параметр `zcol`, а палитру передать в параметр `col.regions`:

```{r, eval = FALSE}
nconts = length(unique(countries$continent))
mapview(countries, zcol = 'continent', 
        col.regions = RColorBrewer::brewer.pal(nconts, 'Set1'))
```

```{r, echo = FALSE}
knitr::include_graphics('images/mapview2.png')
```

Чтобы скомбинировать несколько слоев, необходимо сложить несколько вызовов `mapview()`:

```{r, eval=FALSE}
{ mapview(countries, zcol = 'continent', 
        col.regions = RColorBrewer::brewer.pal(nconts, 'Set1')) + 
  mapview(cities_large, col.regions = 'black', label = 'name', cex = 3) } |> 
  leafem::addStaticLabels(cities_large, label = cities_large$name,
                          offset = c(0.1, 0),
                          style = list("color" = "black", "font-weight" = "bold"))
```

```{r, echo = FALSE}
knitr::include_graphics('images/mapview3.png')
```

### Атрибутивные операции {#sf_attrs}

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

```{r, collapse = TRUE, message=FALSE, warning=FALSE}
italy = countries |> 
  filter(sovereignt == 'Italy')
plot(st_geometry(italy))
```

Следующий пример иллюстрирует как выбрать страны с населением более 100 млн человек:

```{r, collapse = TRUE, message=FALSE, warning=FALSE}
largest = countries |> 
  select(pop_est) |> 
  filter(pop_est > 100000000)
plot(outlines, col = 'lightgrey')
plot(largest, col = 'red', add = TRUE)
```

Обратите внимание на то, что при вызове функции `select()` столбец `geometry` не был указан в числе выбираемых переменных. Тем не менее, то, что мы смогли построить карту по результатам выборки, говорит о том, что данный столбец был сохранен. *Функции **dplyr** определены для объектов `sf` таким образом, чтобы всегда сохранять геометрический столбец.*

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

Например, мы можем агрегировать валовой региональный продукт по континентам:

```{r, collapse=TRUE, warning=F}
continents = countries %>% # этот пайп из пакета magrittr для подстановки в точку
  filter(., st_is_valid(.)) |>
  group_by(continent) |>  
  summarise(gdp = sum(gdp_md_est))
plot(continents['gdp'])
```

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

### Создание пространственных объектов {#sf_creation}

Пространственные объекты в R можно собирать "вручную", если есть такая необходимость. Например, вам известны координаты границ участков полевого обследования, полученные посредством GPS, а вам необходимо превратить их в полигоны, чтобы выполнить анализ и картографирование. Придется из координат собрать полигоны программным путем. Процесс создания пространственных объектов осуществляется в последовательности их иерархического соподчинения: **sfg** \> **sfc** \> **sf**.

#### Геометрические объекты (sfg) {#sf_sfg}

Для создания геометрических объектов в пакете sf существует ряд функций с говорящими названиями:

| Функция                   | Тип пространственного объекта |
|---------------------------|-------------------------------|
| `st_point()`              | *POINT*                       |
| `st_linestring()`         | *LINESTRING*                  |
| `st_polygon()`            | *POLYGON*                     |
| `st_multipoint()`         | *MULTIPOINT*                  |
| `st_multilinestring()`    | *MULTILINESTRING*             |
| `st_multipolygon()`       | *MULTIPOLYGON*                |
| `st_geometrycollection()` | *GEOMETRYCOLLECTION*          |

В зависимости от типа создаваемого объекта, данные функции принимают координаты, организованные в виде одной из трех структур данных:

-   Вектор координат (*POINT*)
-   Матрица координат (*MULTIPOINT* или *LINESTRING*), в которой строки соответствуют точкам, столбцы --- координатам
-   Список (для всех остальных типов)

Проще всего создаются отдельные **точки** (*POINT*):

```{r}
st_point(c(0, 2)) # XY POINT
st_point(c(0, 2, -1)) # XYZ POINT
st_point(c(0, 2, 5), dim = 'XYM') # XYM POINT
st_point(c(0, 2, -1, 5)) # XYZM POINT
```

Дополнительный параметр `dim=` служит для уточнения типа геометрии точек и по сути нужен только тогда, когда необходимо создать редко используемые точки типа *XYM*. во всех остальных случаях (*XY*, *XYZ*, *XYZM*) размерность геометрии распознается по умолчанию.

При создании **мультиточек** (*MULTIPOINT*) и **линий** (*LINESTRING*) необходимо подавать на вход функции уже матрицу координат:

```{r}
coords = matrix(c(
  0, 2,
  1, 3,
  3, 1,
  5, 0
), ncol = 2, byrow = TRUE)

mp = st_multipoint(coords) # XY MULTIPOINT
print(mp)

ls = st_linestring(coords) # XY LINESTRING
print(ls)
```

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

```{r}
plot(ls)
plot(mp, col = 'red', pch = 19, add = TRUE)
```

Создание трех-(*XYZ*, *XYM*) и четырехмерных (*ZYXM*) мультиточек и линий выполняется аналогично, но матрица должна содержать не 2, а, соответственно 3 или 4 столбца, и при необходимости параметр `dim = 'XYM'`.

Создание **полигонов** (*POLYGON*), **мультиполигонов** (*MULTIPOLYGON*) и **мультилиний** (*MULTILINESTRING*) требует уже создания списков из матриц.

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

-   Первая матрица координат в списке отвечает за контур полигона
-   Все остальные матрицы координат отвечают за дыры в полигоне
-   Координаты первой и последней точки в каждой матрице должны совпадать

Если дыр в полигоне нет, его список будет содержать только одну матрицу. Рассмотрим оба примера построения **полигонов**:

```{r}
coords = matrix(c( # Координаты главного полигона
  1, 0,
  0, 2,
  2, 3,
  4, 2,
  3, 0.5,
  1, 0
), ncol = 2, byrow = TRUE)

pol = st_polygon(list(coords)) # Простой полигон
print(pol)

plot(pol, col = 'lightblue')

hole = matrix(c( # Координаты дыры
  2, 1,
  3, 1.5,
  3, 2,
  2, 2,
  1.5, 1.5,
  2, 1
), ncol = 2, byrow = TRUE)

pol2 = st_polygon(list(coords, hole)) # Полигон с дырой
print(pol2)

plot(pol2, col = 'lightblue')
```

Мультиполигоны (*MULTIPOLYGON*) и мультилинии (*MULTILINESTRING*) требуются тогда, когда один и тот же географический объект состоит из нескольких геометрических объектов. Простейший пример --- островные государства. Чтобы представить страну, занимающую архипелаг (Багамские острова, Индонезия, Япония и т.д.) как один пространственный объект, необходимо создать мультиполигон. Все компоненты мультиполигона будут иметь общий набор атрибутов (непространственных характеристик). Мультилинии используются реже мультиполигонов и необходимы для представления линейных объектов, разорванных в пространстве. Примером такого объекта может быть любая река или канал, которые разорваны в тех местах, где они протекают через озеро или водохранилище, представленное полигональным объектом.

В мультиполигонах добавляется еще один уровень списка, то есть искомые матрицы координат будут располагаться как минимум на втором уровне вложенности:

```{r}
coords1 = matrix(c(
  0.5, 0,
  0, 1,
  1, 1.5,
  2, 1,
  1.5, 0.25,
  0.5, 0
), ncol = 2, byrow = TRUE)

coords2 = matrix(c(
  3, 1,
  2.5, 2,
  3.5, 2.5,
  4, 2,
  4, 1.25,
  3, 1
), ncol = 2, byrow = TRUE)

mpol = st_multipolygon(list(list(coords1), list(coords2)))
print(mpol)

plot(pol, col = 'grey') # Обычный полигон (серый)
plot(mpol, col = 'pink', add = TRUE) # Мультиполигон (розовый)
```

Как насчет острова на озере? Если остров и суша, окружающая озеро, составляют единое целое (например, подлежат учету как единый массив леса), их можно собрать как мультиполигон. В этом случае первая компонента мультиполигона будет представлять собой полигон с дыркой, а вторая компонента --- остров. Порядок компонент в данном случае роли не играет:

```{r}
coords4 = matrix(c(
  2.2, 1.2,
  2.8, 1.5,
  2.8, 1.8,
  2.2, 1.8,
  2.0, 1.6,
  2.2, 1.2
), ncol = 2, byrow = TRUE)

island = st_polygon(list(coords4))

mpol2 = st_multipolygon(list(pol2, island))

print(mpol2)

plot(mpol2, col = 'darkolivegreen4')
```

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

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

```{r}
coords1 = matrix(c(
  -3, 0,
  -1, 2,
  0, 2
), ncol = 2, byrow = TRUE)

coords2 = matrix(c(
  4, 2,
  5, 3,
  6, 5
), ncol = 2, byrow = TRUE)

mline = st_multilinestring(list(coords1, coords2))
print(mline)

plot(mline, lwd = 3, col = 'blue')
plot(pol2, col = 'lightblue', add = TRUE)
```

Наконец, еще один вид геометрии --- это геометрическая коллекция (GEOMETRYCOLLECTION), который позволяет хранить вместе любые виды геометрий. Эта возможность используется достаточно редко, тем не менее, рассмотреть ее нужно. Геометрическая коллекция собирается из списка объектов с простыми типами геометрии (мы создали их ранее):

```{r}
col = st_geometrycollection(list(ls, mp, mline, pol2))
print(col)
plot(col)
```

#### Списки геометрических объектов (sfc) {#sf_sfc}

Списки геометрических объектов (класс `sfc`) используются в таблицах пространственных объектов в качестве столбца, который хранит геометрию объектов. Создание таких списков осуществляется функцией `st_sfc()`, которой достаточно передать в качестве перечня параметров объекты типа `sfg`. Рассмотрим создание списка геометрий на примере точечных объектов (для остальных типов объектов порядок действий не меняется):

```{r}
moscow.sfg = st_point(c(37.615, 55.752))
irkutsk.sfg = st_point(c(104.296, 52.298))
petro.sfg = st_point(c(158.651, 53.044))

cities.sfc = st_sfc(moscow.sfg, irkutsk.sfg, petro.sfg)
print(cities.sfc)
```

При создании списка геометрий для него может быть определена система координат (это можно сделать и позднее при создании таблицы пространственных объектов). Для этого используем уже знакомую нам функцию `st_crs()`:

```{r}
st_crs(cities.sfc) = st_crs(4326) # WGS84
print(cities.sfc)
```

> Для списка геометрий может быть определена только одна система координат

Можно посмотреть, куда легли наши точки:

```{r}
plot(cities.sfc, pch = 19)
countries |> 
  filter(sovereignt == 'Russia') |> 
  st_geometry() |>
  plot(add = TRUE)
```

#### Пространственные объекты (sf) {#sf_sf}

Пространственные объекты (класс `sf`) организуются в виде фрейма данных, один из столбцов которого имеет класс `sfc`. Для этого следует сначала создать обычный фрейм данных с атрибутами, а затем соединить его со списком геометрий посредством функции `st_sf`:

```{r}
city.attr = data.frame(
  name = c('Москва', 'Иркутск', 'Петропавловск-Камчатский'),
  established = c(1147, 1661, 1740),
  population = c(12500, 620, 180)
)

cites.sf = st_sf(city.attr, geometry = cities.sfc)
print(cites.sf)
```

#### Точки по координатам {#sf_geom_points}

Достаточно распространенной является следующая задача: имеются координаты точек в табличной форме, необходимо создать на их основе набор пространственных объектов. Для решения этой задачи можно воспользоваться функцией `st_as_sf()`. Рассмотрим задачу на примере файла координат станций из базы метеорологических данных [**ВНИИГМИ-МЦД**](http://meteo.ru/):

```{r}
(stations = read_fwf('data/vniigmi/stations.txt', 
                    col_positions = fwf_widths(diff(c(1, 7, 42, 47, 53, 59, 67, 71)), 
                                               col_names = c('id', 'name', 'lat', 't1', 'lon', 't2', 'z')),
                    locale = locale(encoding = 'CP1251')))
```

Теперь создадим пространственные точки на основе этой таблицы, взяв координаты из столбцов *lat* и *lon* соответственно и указав код системы координат:

```{r}
sf_stations = st_as_sf(stations, coords = c("lon", "lat"), crs = 4326)

plot(st_geometry(sf_stations), pch = 19, col = 'red', cex = 0.25)
plot(st_geometry(countries), border = 'grey', add = TRUE)
box()
```

#### Преобразование типов геометрии {#sf_cast}

Для преобразования типов геометрии существует функция `st_cast()`. Функция принимает объекты классов `sfg`, `sfc` или `sf`, а также название типа геометрии, к которому необходимо привести входные объекты. Довольно часто возникает задача конвертации площадного объекта в линейный и обратно, а также задача получения координат вершин линейного или площадного объекта в виде точек. Примеры преобразований:

```{r}
italy.borders = st_cast(italy, 'MULTILINESTRING')
class(st_geometry(italy.borders))

italy.regions = st_cast(italy.borders, 'MULTIPOLYGON')
class(st_geometry(italy.regions))

italy.points = st_cast(italy.borders, 'POINT')
class(st_geometry(italy.points))

plot(st_geometry(italy.regions), lwd = 0.5)
plot(italy.points, pch = 20, add = TRUE)
```

#### Полигонизация и разбиение линий {#sf_polygonize}

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

1.  Каждая линия является простой (не имеет самопересечений)
2.  Линии касаются только своими начальными и конечными точками
3.  Линии образуют замкнутую последовательность (т.е. выйдя из любой конечной точки и двигаясь вдоль множества линий, можно вернуться в ту же точку.)

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

Рассмотрим операции полигонизации и добавления узлов на простом примере трех пересекающихся отрезков:

```{r}
# Создадим три линии
coords1 = rbind(c(0, 0), c(0, 6))
line1 = st_linestring(coords1)

coords2 = rbind(c(-1,1), c(5,1))
line2 = st_linestring(coords2)

coords3 = rbind(c(-1,5), c(4,0))
line3 = st_linestring(coords3)

# Создадим мультилинию
mls = st_multilinestring(list(line1, line2, line3))
plot(mls)

# Посмотрим на ее точки
points = st_cast(mls, 'MULTIPOINT')
plot(points, pch = 20, add = TRUE)
```

Из рисунка видно, что линии образуют треугольную замкнутую область. Также рисунок показывает, что у компонент мультилинии нет вершин в точках пересечения. Мы можем попытаться найти замкнутые области и превратить их в полигоны, используя `st_polygonize()`:

```{r}
st_polygonize(mls)
```

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

```{r}
mls2 = st_node(mls)
poly2 = st_polygonize(mls2)
points2 = st_cast(mls2, 'MULTIPOINT')

plot(mls2)
plot(poly2, col = 'grey', add = TRUE)
plot(points2, pch = 20, add = TRUE)
```

Таким образом, после разбиения линий на куски в точках пересечения стала возможной операция полигонизации.

### Геометрические атрибуты {#sf_geom_attrs}

К описательным характеристикам геометрии относятся ограничивающий прямоугольник, периметр (для линий и полигонов), площадь (для полигонов), центроид и список координат, которые можно получить с помощью функций `st_bbox()`, `st_length()`, `st_area()`, `st_centroid()` и `st_coordinates()` соответственно. Функции корректно работают для простых объектов, мультиобъектов, списков геометрий и пространственных объектов. Применительно к полигону Италии эти параметры будут учитывать части геометрии, занимаемые островами:

```{r}
st_bbox(italy)        # Координаты органичивающего прямоугольника
st_area(italy)        # Площадь
st_length(italy)      # Периметр
st_centroid(italy) |> st_geometry()    # Центроид (может быть не внутри для невыпуклых фигур)
st_point_on_surface(italy) |> st_geometry() # Точка гарантированно внутри, но не обязательно в центре
st_coordinates(italy) |> head() # Список координат
```

Обратите внимание на то, что площадь и периметр выводятся с указанием единиц измерений! Это возможно благодаря тому, что объекты типа `sf` поддерживают единицы измерений на основе пакета [units](https://cran.r-project.org/web/packages/units/index.html).

> Если данные находятся в плоской прямоугольной системе координат, то единицы измерения как правило указываются в параметрах проекции --- следовательно, они могут быть использованы при вычислении геометрических параметров объектов. Если же данные хранятся в широтах и долготах, то вычисление геометрических параметров осуществляется пакетом *sf* по формулам сферической тригонометрии через пакет [geosphere](https://cran.r-project.org/web/packages/geosphere/index.html). Это позволяет выводить результат в плоских единицах измерения.

Ограничивающий прямоугольник можно быстро преобразовать в полигон и нанести на карту, применив функцию `st_as_sfc()`:

```{r}
box = st_as_sfc(st_bbox(italy)) # Ограничивающий прямоугольник

plot(italy |> st_geometry(), 
     col = 'lightgrey')
plot(box, 
     border = 'orangered', 
     add = TRUE)
plot(st_centroid(italy), 
     col = 'darkgreen', 
     pch = 19,
     add = TRUE)
plot(st_point_on_surface(italy), 
     col = 'steelblue4', 
     pch = 19,
     add = TRUE)
```

Как видно, в данном случае центроид и характерная точка расположились относительно рядом. Однако так бывает далеко не всегда. Выполним аналогичные вычисления для Индонезии:

```{r}
indonesia = countries |> filter(sovereignt == 'Indonesia')

box = st_as_sfc(st_bbox(indonesia))

plot(indonesia |> st_geometry(), 
     col = 'lightgrey')
plot(box, 
     border = 'red', 
     add = TRUE)
plot(st_centroid(indonesia), 
     col = 'darkgreen', 
     pch = 19,
     add = TRUE)
plot(st_point_on_surface(indonesia), 
     col = 'steelblue4', 
     pch = 19,
     add = TRUE)
```

Как видно, в данном случае центроид мультиполигона оказался за пределами какой-либо из его полигональных компонент, в то время как характерная точка находится внутри одного из полигонов. Таким образом, если необходимо получить точку, находящуюся гарантированно в пределах исходного множества, следует использовать `st_point_on_surface()`. При этом следует помнить, что характерная точка, в отличие от центроида, может не располагаться в визуальном центре тяжести множества объектов, и выбор между этими способами описания геометрии остается за разработчиком.

### Экспорт {#sf_export}

Для экспорта векторных пространственных данных можно воспользоваться функцией `st_write()`, которая определит формат выходного файла по указанному вами расширению:

```{r, eval = FALSE}
st_write(cites.sf, 'data/mycities.shp') # Шейп-файл
```

## Растровые данные {#raster_data_r}

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

### Теоретические сведения {#raster_data}

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

В отличие от векторных данных, которые требуют указания координат для каждой вершины, регулярно-ячеистый характер растровой модели позволяет вычислять координаты пикселов на основе их индексов. Поэтому фактически растровые данные хранятся в виде линейно упорядоченного списка значений *(raster values)* и описания геометрии растра *(raster geometry)*.

**Геометрия растра** определяет, где именно располагаются в пространстве пикселы растра и может быть описана путем указания следующих компонент[^09-spatialdata-3]:

[^09-spatialdata-3]: Названия перечисленных компонент геометрии растра укоренились благодаря распространенности стандарта [Esri ASCII Grid](https://en.wikipedia.org/wiki/Esri_grid)

| Параметр    | Назначение                                       |
|-------------|--------------------------------------------------|
| `NCOLS`     | Количество столбцов                              |
| `NROWS`     | Количество строк                                 |
| `XLLCENTER` | Координата $X$ центра левой нижней ячейки растра |
| `YLLCENTER` | Координата $Y$ центра левой нижней ячейки растра |
| `CELLSIZE`  | Размер ячейки                                    |

Иногда вместо параметров `XLLCENTER`/`YLLCENTER` указываются `XLLCORNER`/`YLLCORNER`, которые кодируют координаты левого нижнего угла, а не центра левой нижней ячейки растра. Выбор одного из двух этих вариантов определяет тип *регистрации растра*, а их значения указывают, в какое именно место необходимо "посадить" растр, чтобы его ячейки заняли соответствующие им области в системе координат. Если геометрия растра характеризуется *анизотропией*, то вместо одного значения `CELLSIZE` могут быть указаны разные размеры ячеек по осям координат `CELLSIZEX` и `CELLSIZEY`.

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

В настоящее время для работы с растровыми данными в R используются два пакета: [**stars**](https://r-spatial.github.io/stars/) и [**terra**](https://cran.r-project.org/web/packages/terra/index.html). **terra** является наследником пакета [raster](https://cran.r-project.org/web/packages/terra/index.html), который исторически был основным средством работы с растровыми данными и обладает широким спектром функций растрового анализа. **stars** --- относительно новый, разработан с целью поддержки многомерных данных и более тесного взаимодействия с пакетом `sf`. В целом можно сказать, что пакеты **terra** и **stars** частично пересекаются по функциональности, но скорее дополняют друг друга, нежели дублируют.

В этой и ближайших лекциях мы будем работать с растрами в формате **stars**, поскольку он концептуально близок к пакету **sf**.

### Чтение {#raster_read}

Для чтения растров любой размерности можно использовать функцию `read_stars()`:

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

img = read_stars('data/world/BlueMarbleJuly.tif') # Цветной космический снимок (RGB)
img
```

### Внутренняя структура {#stars_inner}

Для работы с данными типа stars необходимо понимать их внутреннюю структуру. Для начала можно взглянуть на нее посредством стандартной функции `str()`:

```{r}
str(img)
```

Видно, что данный трёхканальный растр представляет собой список из единственного элемента с названием `BlueMarbleJuly.tif` --- это имя было присвоено автоматически при чтении растра. Каждый такой элемент соответствует *переменной* данных. В данном случае переменная одна --- это интенсивность цвета. Хранится она в виде трехмерного массива (`array`) размерностью $720 \times 360 \times 3$:

```{r}
str(img[[1]])
img[[1]][100, 200, 2]
```

Каждой оси этого массива соответствует измерение (`dimension`), которое определяет параметры отображения индексов массива на соответствующую систему координат (пространственную, временную, спектральную и т.д.). Например, чтобы понять, что ячейка растра с индексами `[36, 18, ]` имеет географические координаты (широту и долготу) `(0, 0)`, нужно знать направления осей растра, размер ячейки и координаты одной из угловых ячеек растра. Необходимая информация находится в атрибуте `dimensions` объекта `stars`, т.е. является общей для всех переменных. При печати параметры измерений выводятся в удобном табличном виде:

```{r}
attr(img, 'dimensions')
```

Этот атрибут представляет собой список, длина которого равна количеству измерений в массиве данных переменной. Обычно измерения имеют имена, в данном случае это `x`, `y` и `band`. Описание каждого измерения выполнено по единому шаблону, который включает следующие параметры:

-   `from`: начальный индекс (будет меняться при обрезке растра при постоянной точке отсчета индексов);
-   `to`: конечный индекс (будет меняться при обрезке растра при постоянной точке отсчета индексов);
-   `offset`: координата первого пиксела (точки отсчета);
-   `delta`: размер ячейки;
-   `refsys`: координатная (референцная) система: для систем счета координат, времени, высот и других измерений будет своя;
-   `point`: логическое значение, которое указывает, следует ли интерпретировать элементы растра по этой оси как измеренные в точке (мгновенные) или агрегированные по площади (за временной период);
-   `values`: значения координат ячеек по данной оси
    -   `NULL` (используется в большинстве случаев, т.к. координаты могут быть вычислены на основе `from`, `delta` и индекса пиксела),
    -   вектор координат (используется для представления ректилинейных растров в переменным размром пиксела),
    -   объект класса `intervals` (список из двух векторов --- начал и концов интервалов), or
    -   матрица координат такой же размерности, что и пространственные измерения растра. В случае стекущего примера будет иметь размер $720 \times 360$. Используется для представления *криволинейных* растров.

Например, посмотрим параметры измерения `x` растра:

```{r}
attr(img, 'dimensions')[['x']]
```

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

```{r}
img |> attr('dimensions') |> attr('raster') |> str()
```

Видно, что атрибут `raster` содержит 3 элемента:

-   `dimensions`: названия измерений, которые являются пространственными
-   `affine`: параметры аффинного преобразования, которое будет применяться к пространственным измерениям перед их отображением или применением в операциях пространственного анализа
-   `curvilinear`: логическое значение, которое устанавливает, является ли растр криволинейным (в этом случае в параметре `values` пространственных измерений должна быть матрица координат)

### Визуализация {#raster_viz}

#### Статичные карты {#raster_viz_static}

Для визуализации одноканальных растров используется функция `plot()`. В простейшем виде ей достаточно просто передать визуализируемый растр:

```{r}
par(mfrow = c(1,1))
plot(dem)
```

Поскольку растры часто используют в классифицированном виде, вы можете сформировать вектор граничных значений классов, вектор цветов классов, и передать их в параметры `breaks` и `col` функции `plot()` соответственно. Если параметр `breaks` не определять, то весь диапазон значений растра будет разбит на равные интервалы соответственно количеству цветов. Если не определять параметр `col`, то будет применена стандартная палитра `terrain.colors`. Вы также можете использовать одну из готовых палитр цветов или создать ее вручную (см. посвященную графической подсистеме R):

```{r}
brks = c(-12000, -5000, -2500, -1000, -200, 0, 200, 500, 1000, 2000, 4000, 8000)
clrs = c(
  "steelblue4",
  "steelblue3",
  "steelblue2",
  "steelblue1",
  "lightskyblue1",
  "darkseagreen",
  "lightgoldenrod1",
  "darkgoldenrod1",
  "darkorange",
  "coral2",
  "firebrick3")

plot(dem, breaks = brks, col = clrs)

plot(dem, col = colorRampPalette(c("black", "white"))(255))

plot(dem, col = rainbow(10))
```

Для синтезирования цветного изображения на основе многоканального растра необходимо объект `stars` предварительно подать в функцию `st_rgb()`:

```{r}
plot(st_rgb(img))
```

Поскольку при визуализации космических снимков часто используют различные варианты синтеза каналов (чтобы лучше дешифрировать те или иные категории объектов), функция `st_rgb()` предоставляет такую возможность. Достаточно перечислить последовательность каналов растрового стека (по умолчанию эти каналы будут подставлены в каналы R, G и B соответственно):

```{r}
st_rgb(img[,,,c(1, 2, 3)]) |> plot()
st_rgb(img[,,,c(1, 3, 2)]) |> plot()
st_rgb(img[,,,c(2, 1, 3)]) |> plot()
st_rgb(img[,,,c(2, 3, 1)]) |> plot()
st_rgb(img[,,,c(3, 1, 2)]) |> plot()
st_rgb(img[,,,c(3, 2, 1)]) |> plot()
```

Вы можете совмещать на картах несколько растровых и векторных слоев точно так же как и при совмещении векторных данных (указав параметр `add = TRUE` при вызове функции `plot()`):

```{r}
plot(st_rgb(img), reset = FALSE)
plot(outlines, border = rgb(1,1,1,0.5), lwd = 0.5, add = TRUE)
```

#### Интерактивные карты {#raster_viz_inter}

Объекты типа stars могут быть визуализированы аналогично векторным на интерактивных картах `mapview`:

```{r, eval = F}
mapview(dem, at = brks, col = clrs)
```

```{r, echo = FALSE}
knitr::include_graphics('images/mapview4.png')
```

### Обрезка {#raster_crop}

Одна из распространенных задач при работе с растрами --- это обрезка, то есть удаление растровых данных, находящихся за пределами указанной территории. Чаще всего обрезку делают либо ограничивающим прямоугольником, либо полигональным объектом. Рассмотрим оба варианта:

```{r}
# Обрезка по ограничивающему прямоугольнику
box = st_bbox(c(xmin = -80, xmax = -10, ymax = 85, ymin = 58), crs = st_crs(4326))
dem_greenland = dem[box]

dem_greenland

plot(dem_greenland)
```

Аналогичным образом можено обрезать растр контуром выбранной страны:

```{r}
country = countries |> 
  filter(name == 'Afghanistan')
dem_country = dem[country]
dem_country
plot(dem_country)
```

### Индексирование

Ортогональная структура объектов типа stars позволяет выполнять по ним различные срезы, отсекая ненужные данные. Для этого используется привычный по работе с векторами оператор квадратной скобки `[`, который работает следующим образом:

-   первый аргумент выбирает атрибут
-   второй и последующий аргументы выбирают измерения.

Таким образом, при работе с растрами, которые содержат один атрибут, вам необходимо указать 4 индекса: `[var, x, y, band]`, где `var` - это название или порядковый номер атрибута, а `x, y, band` --- порядковые номера двух координатных и одного семантического измерения.

Например:

```{r}
# выбрать 1 канал
ch1 = img[,,,1] 
ch1
plot(ch1)

# выбрать диапазон ячеек растра
frag = img[, 320:470, 100:255, ] 
frag
plot(st_rgb(frag))
```

### Манипуляции

Объекты типа `stars` поддерживают манипуляции, аналогичные тем, что могут применяться к векторным данным. Посмотрим это на примере данных по высоте земной поверхности с учетом и без покровного оледенения:

```{r}
etopo = read_stars(c('data/etopo1_bed.tif', 'data/etopo1_ice.tif'))
etopo
```

Для начала переименуем переменные:

```{r}
etopo = etopo |> setNames(c('bed', 'ice'))
etopo
```

После этого посчитаем, например, толщину покровного оледеления как разность `ice` и `bed` через *мутирование*:

```{r}
etopo = etopo |>
  mutate(depth = ice - bed)

plot(etopo['depth'], 
     col = cm.colors(5),
     breaks = c(0, 500, 1000, 2000, 3000, 4000),
     main = 'Мощность покровного оледенения',
     reset = FALSE)
plot(oceans, col = 'steelblue', add = TRUE)
```

*Фильтрация* происходит по измерениям, но применяется не к индексам ячеек, а к соответствующим величинам измерений:

```{r}
greenland = etopo |> 
  filter(x > -80, x < -10, y > 58, y < 85)

plot(greenland)
```

*Выбор* переменной позволяет оставить только ее:

```{r}
icedepth = etopo |> 
  select(depth)

icedepth
```

### Экспорт {#raster_export}

Чтобы экспортировать (сохранить в файл) любой растр, можно воспользоваться функцией `write_stars()`, указав имя выходного файла:

```{r, eval = FALSE}
write_stars(greenland, 'data/world/greenland.tif')
```

## Пространственная привязка {#spatref}

### Компоненты пространственной привязки {#spatref_components}

**Пространственная привязка** (*spatial reference* или *georeference*) --- важнейшая составляющая пространственных данных, которая говорит нам о том, как правильно интерпретировать координаты объектов. Пространственная привязка в простейшем случае включает несколько фундаментальных компонент:

1.  *Эллипсоид вращения* --- тело, по отношению к которому вычисляются геодезические координаты точек (широты и долготы)
2.  *Исходные геодезические даты (датум)* --- параметры положения эллипсоида в теле Земли
3.  *Географическая система координат* --- включает датум, положение начального меридиана и единицы измерения широт и долгот
4.  *Проекция* --- математический способ перехода от географических координат на эллипсоиде к плоским прямоугольным координатам карты.
5.  *Плоская прямоугольная система координат* --- включает проекцию, ее параметры и единицы измерения координат.

Если точки имеют также координаты $Z$, то для их правильной интерпретации необходимы дополнительные компоненты пространственной привязки:

1.  *Система счета высот* (геодезические, нормальные, ортометрические) - определяют содержательный смысл и порядок вычисления высот и глубин (координата Z)
2.  *Модель геоида, квазигеоида или эллипсоида* --- определяет поверхность, относительно которой вычисляются высоты точек.
3.  *Вертикальная система координат* --- фактическая реализация системы счета высот относительно конкретной поверхности относимости с заданным положением нулевого уровня. Например, в России это *Балтийская система нормальных высот* с нулем в г. Кронштадт.

Аналогичным образом требуется введение системы счета дополнительных координат $M$, если они используются в представлении координат.

### Форматы описания пространственной привязки {#spatref_formats}

Существует три распространенных способа задания (хранения) пространственной привязки:

-   *PROJ.4 String* --- представление в виде строки.
-   *WKT (Well-Known Text)* --- представление в виде иерархического списка. Это *наиболее полный* формат описания пространственной привязки, который рекомендуется к использованию для избежания неоднозначностей.
-   *EPSG (European Petroleum Survey Group)* --- представление в виде числового кода.

Для поиска проекций в перечисленных форматах представления удобно воспользоваться порталом [spatialreference.org](spatialreference.org).

**PROJ.4 String** --- строковый формат представления информации о пространственной привязки, используемый в библиотеке [**PROJ**](http://proj.org). Данная библиотека лежит в основе координатных систем пространственных данных, используется в **R**, **Python**, **QGIS** и прочих средах. Основные параметры строки:

```         
+datum     Datum name (see `proj -ld`)
+ellps     Ellipsoid name (see `proj -le`)
+lat_0     Latitude of origin
+lat_1     Latitude of first standard parallel
+lat_2     Latitude of second standard parallel
+lat_ts    Latitude of true scale
+lon_0     Central meridian
+proj      Projection name (see `proj -l`)
+units     meters, US survey feet, etc.
+vunits    vertical units.
+x_0       False easting
+y_0       False northing
+zone      UTM zone
```

Примеры записи координат в формате PROJ.4:

-   Географические координаты в *WGS84* (без проекции):

```{r, echo = FALSE}
cat(st_crs(4326)[['proj4string']])
```

-   Координаты в проекции *Web Mercator* (проекция Google Maps, Яндекс.Карт и т.д.):

```{r, echo = FALSE}
cat(st_crs(3857)[['proj4string']])
```

-   Координаты в *конической равнопромежуточной проекции*:

```{r, echo = FALSE}
cat(st_crs(54027)[['proj4string']])
```

-   Координаты в проекции *UTM, зона 37*:

```{r, echo = FALSE}
cat(st_crs(32637)[['proj4string']])
```

**WKT** предполагает представление вышеуказанных компонент пространственной привязки к виде иерархического списка. Например, так будет выглядеть информация о *полярной стереографической проекции для карт России*:

```         
PROJCS["WGS 84 / EPSG Russia Polar Stereographic",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",90],
    PARAMETER["central_meridian",105],
    PARAMETER["scale_factor",0.994],
    PARAMETER["false_easting",2000000],
    PARAMETER["false_northing",2000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","5940"]]
```

**EPSG (European Petroleum Survey Group)** --- европейская рабочая группа нефтегазовой области, которая ведет реестр систем координат с уникальными цифровыми кодами вида `EPSG:xxxxxx`. Коды EPSG оказались настолько удобны, что используются повсеместно для быстрой инициализации проекций со стандартными параметрами. Например, вышеприведенные проекции имеют следующие коды EPSG:

-   *WGS84*: `EPSG:4326`
-   *Web Mercator*: `EPSG:3857`
-   *UTM*: `EPSG:326..` , например для UTM 37N: `EPSG:32637`

### Преобразование координат {#spatref_transform}

Преобразование координат включает три различных операции:

1.  **Трансформирование** --- пересчет географических координат с одного датума на другой

2.  **Проецирование** --- переход от географических координат к плоским прямоугольным

3.  **Обратное проецирование** --- переход от плоских координат к географическим.

Например, чтобы пересчитать координаты *UTM* в проекцию *Гаусса-Крюгера*, необходимо:

1.  Обратно проецировать координаты в географические *WGS84*
2.  Трансформировать географические координаты c *WGS84* в *ГСК-2011*
3.  Проецировать координаты *ГСК-2011* в проекцию *Гаусса-Крюгера*

*Несоответствие датумов часто является причиной того, что данные из разных наборов плохо совмещаются друг с другом*

### Работа с пространственной привязкой в R {#sf_crs}

Работа с пространственной привязкой данных в R состоит в основном из четырех операций:

-   чтение информации о системе координат
-   создание информации о системе координат
-   замена информации о системе координат
-   изменение системы координат (проецирование)

Первые три операции (чтение, создание, замена) осуществляются функцией `st_crs()`. Чтобы прочитать информацию о проекции, достаточно передать в качестве параметра объект типа `sf`:

```{r}
st_crs(countries)    # Координатная система
```

Эта же функция позволяет создать новую координатную систему, путем передачи ей кода *EPSG* или строки *PROJ.4*:

```{r}
st_crs(3857) # Проекция Меркатора для карт мира
st_crs('+proj=robin') # Проекция Робинсона для карт мира

# Проекция UTM, зона 37.
st_crs('+proj=utm +zone=37 +datum=WGS84 +units=m')
```

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

```{r}
st_crs(countries) = NA
st_crs(countries) 

st_crs(countries) = st_crs(4326)
st_crs(countries)
```

> **Внимание**: замена координатной системы не осуществляет перепроецирования данных и не меняет координаты точек. Она лишь влияет на то, как эти координаты будут интерпретироваться. Если вместо проецирования выполнить замену информации о координатной системе, данные будут позиционироваться в неправильном месте.

Для проецирования данных в другую систему координат следует использовать функцию `st_tranform(x, crs)`. Данная функция принимает в качестве параметров класс объектов *sf* и координатную систему, в которую необходимо проецировать данные.

```{r, fig.width=14, fig.height=16, crop = TRUE}
# Проекция Меркатора
countries_merc = st_transform(countries, 3857)
par(mar = c(2,16,2,16))
plot(st_geometry(countries_merc), 
     col = 'lightgray',
     lwd = 0.5,
     graticule = TRUE, 
     axes = TRUE)
```

```{r, fig.width=14, fig.height=10, crop = TRUE}
# Проекция Мольвейде (используем dplyr)
countries_moll = countries |> st_transform('+proj=moll')
plot(st_geometry(countries_moll), 
     col = 'lightgray',
     lwd = 0.5,
     graticule = TRUE, 
     axes = TRUE)
```

```{r, fig.width=14, fig.height=16, crop = TRUE}
# Зарубежная Европа в Конической равнопромежуточной проекции. 
# Задаем только необходимые параметры проекции
countries_eqdc = countries |> 
  filter(continent == 'Europe' & sovereignt != 'Russia') |> 
  st_transform('+proj=eqdc +lon_0=10 +lat_1=30 +lat_2=60 +datum=WGS84 +units=m')

plot(st_geometry(countries_eqdc), 
     col = 'lightgray',
     lwd = 0.5,
     graticule = TRUE, 
     axes = TRUE)
```

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

```{r}
st_crs(dem) # читаем систему координат
st_crs(dem) = NA # очищаем систему координат
st_crs(dem)
st_crs(dem) = st_crs(4326) # создаем систему координат
st_crs(dem)
```

Для проецирования растра в новую систему координат необходимо использовать функцию `st_warp()` Приведем несколько примеров проецирования.

Проекция *Меркатора*:

```{r, fig.width = 14, dpi = 300, crop = TRUE}
img_merc = st_warp(img, crs = 3857)

plot(st_rgb(img_merc), main = NULL, reset = FALSE)
plot(st_geometry(countries_merc), 
     border = rgb(1,1,1,0.5), lwd = 0.25, add = TRUE)
```

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

```{r, fig.width = 14, dpi = 300, crop = TRUE}
img_moll = st_warp(img, crs = st_crs('+proj=moll'), use_gdal = TRUE)
plot(st_rgb(img_moll, probs = c(0.01, 0.99),
       stretch = "percent"), main = NULL, reset = FALSE)
plot(st_geometry(countries_moll), 
     border = rgb(1,1,1,0.5), lwd = 0.5, add = TRUE)
```

Проекция *коническая равнопромежуточная*:

```{r, fig.width = 14, out.width = '50%', dpi = 300, crop = TRUE}
prj = '+proj=eqdc +lon_0=10 +lat_1=30 +lat_2=60 +datum=WGS84 +units=m'

img_eqdc = st_warp(
  img, 
  crs = st_crs(prj)
)

img_eqdc_euro = img_eqdc[st_bbox(countries_eqdc)]

plot(st_rgb(img_eqdc_euro), main = NULL, reset = FALSE)
plot(st_geometry(countries_eqdc), 
     border = rgb(1,1,1,0.5), lwd = 0.5, add = TRUE)
```

Более подробно вопросы выбора проекций и построения сеток координат рассматриваются в следующей главе.

## Краткий обзор {#temporal_review}

Для просмотра презентации щелкните на ней один раз левой кнопкой мыши и листайте, используя кнопки на клавиатуре:

```{r, echo=FALSE}
knitr::include_url('https://tsamsonov.github.io/r-geo-course-slides/09_SpatialData.html#1', height = '390px')
```

> Презентацию можно открыть в отдельном окне или вкладке браузере. Для этого щелкните по ней правой кнопкой мыши и выберите соответствующую команду.

## Контрольные вопросы и упражнения {#questions_tasks_spatial}

### Вопросы {#questions_spatial}

1.  Что такое пространственные данные и какие модели пространственных данных существуют?
2.  Назовите номер стандарта ISO, в котором описана модель *Simple Features*.
3.  Перечислите основные принципы представления объектов в рамках стандарта *Simple Features*.
4.  Какие размерности координат допустимы в объектах *Simple Features*?
5.  Перечислите основные 7 типов геометрий. Сколько всего их описано в стандарте *Simple Features*?
6.  Как называются основные два формата представления объектов *Simple Features*?
7.  Перечислите основные компоненты пространственной привязки.
8.  Перечислите основные форматы описания пространственной привязки.
9.  Дайте расшифровку основных параметров строки *PROJ.4*.
10. Какой номер *EPSG* имеет географическая система координат *WGS84*?
11. В чем отличие трансформирования координат и проецирования?
12. Какие три программных библиотеки составляют основу функциональности пакета **sf**? Каково их назначение?
13. В чем отличие объектов типа **sp** от **sf**?
14. Что означает префикс `st_`, используемый в названиях функций пакета **sf**?
15. Какая функция используется для чтения данных средствами пакета **sf**?
16. Перечислите три класса, слагающих иерархию представления пространственных объектов, реализуемую пакетом **sf**.
17. Какой тип данных имеет колонка с геометрией объекта **sf**?
18. Какая функция позволяет извлечь геометрическую колонку из объекта **sf**?
19. С помощью какой структуры данных фактически реализован класс объектов **sfg**?
20. Сколько карт будет построено функцией `plot()` применительно к объекту **sf**?
21. Как с помощью функции `plot()` нарисовать только геометрию объектов, не отображая атрибутивные характеристики?
22. Какой параметр функции `plot()` отвечает за отображение/не отображение градусной сетки координат?
23. Каким способом можно узнать и задать систему координат объекта **sf**?
24. Какая функция позволяет осуществить проецирование данных?
25. Можно ли применять к объектам типа **sf** стандартные манипуляции **dplyr**?
26. Что произойдет с геометрией пространственных объектов при выполнении агрегирования данных по группам значений заданных атрибутов?
27. Перечислите функции, с помощью которых создаются объекты типа **sfg**, и структуры данных с координатами, которые должны быть поданы на вход этих функций.
28. Назовите три правила, которым подчиняется формат представления координат вершин *полигональных* объектов.
29. Может ли обычный полигон **sf** содержать дырку, или же для этого требуется создание мультиполигона?
30. Как можно быстро собрать слой точечных объектов по их координатам, не собирая объекты вручную?
31. Какая функция позволяет осуществлять преобразование типа геометрии **sf**?
32. Перечислите требования, которым должно удовлетворять множество линейных объектов для того, чтобы к нему была применима операция полигонизации?
33. Назовите функции **sf**, реализующие операцию добавления вершин в точках пересечения линий и операцию полигонизации линий.
34. Перечислите названия функций **sf**, позволяющих получать ограничивающий прямоугольник, периметр, площадь, центроид, характерную точку и координаты объекта.
35. С помощью какой функции осуществляется запись (экспорт) **sf** в файлы пространственных данных?
36. Назовите основные параметры, определяющие геометрию растра.
37. Какие пакеты отвечают за поддержку растровых данных в **R**?
38. Как можно прочитать одноканальный и многоканальный растры в **R**?
39. Какие функции можно использовать для визуализации одноканальных и многоканальных растров?
40. Можно ли совмещать растровые и векторные слои на одном изображении? Если да, то как эта возможность реализуется?
41. Каким образом можно узнать и задать пространственную привязку растрового набора данных?
42. Какая функция отвечает за проецирование растровых данных? Перечислите ее параметры и их назначение.
43. Объясните систему индексирования объектов типа **stars**.
44. Как вычислить новую переменную объекта типа **stars**?
45. Как осуществить экспорт растра в файл?

### Упражнения {#tasks_spatial}

1.  Преобразуйте точки землетрясений из набора данных *quakes* в пространственные объекты и отобразите их сначала средствами стандартной графической подсистемы, а затем на интерактивной карте средствами пакета **mapview**. Передайте магнитуду землетрясения в параметр `zcol` функции `mapview()`, чтобы дифференцировать точки цветом по этому параметру.

2.  Таблица *storms* из пакета **dplyr** содержит данные трекинга тропических циклонов c 1975 по 2015 год. Выберите любой циклон и постройте для него линию трека прохождения и точки прохождения. Отобразите эти данные средствами стандартной графической подсистемы, а затем на интерактивной карте средствами **mapview**. Напишите программу таким образом, чтобы можно было выбирать имя циклона и программа отображала его трек на интерактивной карте.

3.  Скачайте [базу данных](https://github.com/tsamsonov/r-geo-course/blob/master/data/Satino.gpkg) и [цифровую модель рельефа](https://github.com/tsamsonov/r-geo-course/blob/master/data/Satino_DEM.zip) на территорию Сатинского полигона МГУ. Изучите содержимое базы данных и постройте на основе этих данных общегеографическую карту средствами стандартной графической подсистемы.

|                                                                                                                                                                                                                    |
|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| *Самсонов Т.Е.* **Визуализация и анализ географических данных на языке R.** М.: Географический факультет МГУ, `r lubridate::year(Sys.Date())`. DOI: [10.5281/zenodo.901911](https://doi.org/10.5281/zenodo.901911) |