13-VectorAnalysis.qmd

# Векторный анализ {#vector}

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

## Предварительные условия {#mapping3d_prerequisites}

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

```{r}
library(readxl)
library(sf)
library(tidyverse)
library(classInt)
library(mapsf) # Удобное построение тематических карт средствами plot()

library(circular) # статистика направлений
library(NPCirc) # статистика направлений
library(pracma)
```

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

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

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

## Пространственные запросы {#spatial_queries}

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

-   Найти все объекты внутри административного района (топологические отношения)
-   Найти все объекты не далее 100 метров от дороги (метрические отношения)
-   Найти все объекты, расположенные к северу от точки (дирекционные отношения)

Пространственные запросы могут объединять несколько условий. Можно найти объекты, удовлетворяющие одновременно *всем* (*логическое И*) вышеперечисленным условиям: внутри района, не далее 100 м от дороги и к северу от выбранной точки; или *хотя бы одному* (*логическое ИЛИ*) из вышеперечисленных условий. Результат выполнения такого комплексного запроса будет являться, соответственно, пересечением множеств объектов, полученных каждым из запросов, или их объединением.

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

### Контекстные и целевые объекты {#context_objects}

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

-   **контекстные** --- объекты, относительно которых будет оцениваться размещение других объектов, то есть, определяющие *контекст анализа*
-   **целевые** --- объекты, размещение которых анализируется по отношению к контекстным объектам, что является *целью* анализа

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

### Метрические отношение {#sf_spat_metric}

### Топологические отношения {#sf_spat_topology}

Поиск объектов по местоположению базируется на проверке топологических отношений между объектами. Топологические отношения описывают взаимное расположение объектов. Различные варианты топологических отношений для площадных объектов представлены на следующем рисунке, где серым цветом показаны пересечения *внутренних областей* объектов $A$ и $B$, синим цветом --- пересечения *границ* объектов $A$ и $B$:

```{r, echo = FALSE}
def = par(no.readonly = TRUE)

par(mfrow=c(2,4))
par(mar = c(2,1,5,1))

a = list(
  c(0,0,0,20,20,20,20,0,0,0),
  c(0,0,0,20,20,20,20,0,0,0),
  c(0,0,0,20,20,20,20,0,0,0),
  c(0,0,0,20,20,20,20,0,0,0),
  c(0,0,0,20,20,20,20,0,0,0),
  c(0,5,0,15,10,15,10,5,0,5),
  c(0,0,0,20,20,20,20,0,0,0),
  c(2,5,2,15,12,15,12,5,2,5)
)
  
b = list(
  c(30,0,30,20,50,20,50,0,30,0),
  c(10,10,10,30,30,30,30,10,10,10),
  c(0,0,0,20,20,20,20,0,0,0),
  c(20,5,20,15,30,15,30,5,20,5),
  c(10,5,10,15,20,15,20,5,10,5),
  c(0,0,0,20,20,20,20,0,0,0),
  c(8,5,8,15,18,15,18,5,8,5),
  c(0,0,0,20,20,20,20,0,0,0)
)

labels = c(
  'Не пересекается \n(A disjoint B)',
  'Перекрывает \n(A overlaps B)',
  'Совпадает \n(A equals B)',
  'Касается \n(A touches B)',
  'Покрывает \n(A covers B)',
  'Покрыт \n(A covered by B)',
  'Содержит \n(A сontains B)',
  'Содержится \n(A within B)'
)

ashift = c(0, -5, -5, 0, -5, 0, -6, 0)
bshift = c(0, 5, 5, 0, 0, 5, 0, 6)

ncases = length(b)

for (i in 1:ncases) {
  aline = st_linestring(matrix(a[[i]], ncol = 2, byrow = TRUE))
  bline = st_linestring(matrix(b[[i]], ncol = 2, byrow = TRUE))
  
  apol = st_cast(aline, 'POLYGON')
  bpol = st_cast(bline, 'POLYGON')
  
  ac = st_centroid(apol)
  bc = st_centroid(bpol)
  
  geom = st_sfc(list(apol, bpol))
  
  poly = st_intersection(apol, bpol)
  line = st_intersection(bline, aline)
  
  plot(geom, main = labels[i])
  plot(poly, col = 'lightgrey', add = TRUE)
  plot(line, lwd = 4, col = 'blue', add = TRUE)
  text(c(ac[1] + ashift[i], bc[1] + bshift[i]), c(ac[2], bc[2]), labels = c('A', 'B'), cex = 2)
}

par(def)
```

Отношение *Пересекает (intersects)* будет истинно для любого случая когда две геометрии имеют хотя бы одну общую точку, то есть во всех случаях кроме *Не пересекает (disjoint)*. Для проверки этих, а также некоторых других отношений, в пакете `sf` существует ряд функций:

| Функция                      | Топологическое отношение                                                                                                                                      |
|------------------------------|---------------------------------------------------------------------------------------------------------------------------------------------------------------|
| `st_intersects(x, y)`        | `x` имеет общие точки с `y`                                                                                                                                   |
| `st_disjoint(x, y)`          | `x` не имеет общих точек с `y`                                                                                                                                |
| `st_touches(x, y)`           | `x` касается `y` (граница `x` имеет общие точки с границей `y` И внутренняя область `x` не имеет имеет общих точек с внутренней областью `y`)                 |
| `st_crosses(x, y)`           | `x` пересекает `y` (граница `x` имеет общие точки с границей `y`, при этом размерность их пересечения меньше размерности хотя бы одного из исходных объектов) |
| `st_within(x, y)`            | `x` внутри `y` (все точки `x` содержатся в `y` И внутренняя область `x` имеет общие точки с внутренней областью `y`)                                          |
| `st_contains(x, y)`          | `x` содержит `y` (все точки `y` содержатся в `x` И внутренняя область `y` имеет общие точки с внутренней областью `x`)                                        |
| `st_contains_properly(x, y)` | `x` содержит `y` полностью (все точки `y` содержатся в `x` И граница `x` не имеет общих точек с границей `y`)                                                 |
| `st_overlaps(x, y)`          | `x` перекрывает `y` (внутренняя область `x` имеет как общие, так и не общие точки с внутренней областью `y`)                                                  |
| `st_equals(x, y)`            | `x` совпадает `y` (множества точек `x` и `y` совпадают)                                                                                                       |
| `st_covers(x, y)`            | `x` покрывает `y` (все точки `y` содержатся в `x`)                                                                                                            |
| `st_covered_by(x, y)`        | `x` покрыт `y` (все точки `x` содержатся в `y`)                                                                                                               |
| `st_equals_exact(x, y)`      | `x` совпадает `y` точно (упорядоченные множества точек `x` и `y` совпадают)                                                                                   |

Между `covered_by` и `within`, а также `covers` и `contains` нет разницы в случае, когда оба объекта являются площадными. Эта разница будет сказываться если хотя бы один из объектов является линией либо точкой. В этом случае `within`, `contains` и `contains_properly` будут давать ложный результат (FALSE), поскольку ни у линий, ни у точек нет внутренней области.

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

```{r}
countries = st_read('data/ne/countries.gpkg')
outlines = st_geometry(countries)
cities = st_read('data/ne/cities.gpkg')
city.pts = st_geometry(cities)

largest = countries |> 
  select(pop_est) |> 
  filter(pop_est > 100000000)

# Наносим исходную конфигурацию
plot(outlines, lwd = 0.5)
plot(cities, col = 'black', pch = 20, cex = 0.5, add = TRUE)

sf::sf_use_s2(FALSE)

# Отбираем точки внутри стран с максимальным ВВП
sel = cities[largest, ]

# Смотрим что получилось
plot(outlines, lwd = 0.5)
plot(largest, col = 'gray', add = TRUE)
plot(sel, pch = 20, col = 'black', add = TRUE)
```

Разумеется, при выполнении пространственных запросов могут возникать и другие пространственные отношения. Например, мы можем выбрать все страны, имеющие общую границу с Чехией. Для этого можно использовать топологическое отношение `st_touches` вместо `st_intersects` --- это будет гарантировать, что сама Чехия в результате не выберется (касающиеся объекты не могут перекрываться). Тип отношения необходимо поставить в параметр `op =` при выполнении фильтрации фрейма данных:

```{r}
cz = filter(countries, sovereignt == 'Czechia')
neighbors = countries[cz, op = st_touches]

plot(st_geometry(neighbors), col = 'lightgray', lwd = 0.5)
plot(cz, col = 'darkgray', add = TRUE)
```

### Зоны окружения объектов {#zones}

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

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

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

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

Если же при построении зон окружения учитывается взаимное размещение объектов, то в данном случае зоны доступности строятся не исходя из порогового значения показателя (хотя оно может использоваться дополнительно), а исходя из того, какой объект является ближайшим. Конкурентные зоны окружения представляют собой [разбиение](https://ru.wikipedia.org/wiki/Разбиение_множества) пространства на неперекрывающиеся участки без дыр, каждый из которых является зоной окружения соответствующего пространственного объекта. При этом любая точка внутри зоны окружения объекта ближе к этому объекту по выбранному признаку (времени или расстоянию), нежели к любому другому объекту. Конкурентные зоны окружения, построенные по расстоянию, можно реализовать средствами [диаграммы Вороного](https://ru.wikipedia.org/wiki/Диаграмма_Вороного).

## Постановка задач и изучение данных {#use_case}

В настоящем модуле мы рассмотрим вышеперечисленные методы на примере анализа размещения пунктов общественного питания --- кафе, ресторанов и т.д. Используя методы пространственного анализа в среде R, мы ответим на следующие вопросы:

-   Какие улицы являются местами наибольшей концентрации заведений общественного питания?
-   Как распределены заведения общественного питания по районам центра Москвы?
-   Какие заведения общественного питания находятся вблизи метро и на берегу реки?
-   В какие заведения общественного питания можно доехать от выбранной точки в течение 5 минут?
-   Каков оптимальный маршрут между вашим местоположением и заведением, в котором вы хотите пообедать?

В качестве источника данных используем [OpenStreetMap](http://www.openstreetmap.org) --- краудсорсинговый интернет-проект по созданию бесплатных и открытых пространственных данных глобального охвата. Данные OpenStreetMap в удобном для использования в ГИС виде доступны на [портале GIS-Lab](http://beryllium.gis-lab.info/project/osmshp/).

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

-   [osrm](https://cran.r-project.org/web/packages/osrm/index.html) --- построение зон доступности, маршрутов и матриц корреспонденции онлайн на основе данных OpenStreetMap и [OSRM API](http://project-osrm.org).
-   [mapsf](https://riatelab.github.io/mapsf/) --- пакет, облегчающий построение тематических карт и легенд средствами стандартной функции `plot()`.

Начнем наше исследование с визуального анализа исходных данных

```{r, message = FALSE, results = "hide", collapse=TRUE}
# Чтение данных
roads = read_sf("data/roads.gpkg") # Дороги
poi = read_sf("data/poi_point.gpkg") # Точки интереса
rayons = read_sf("data/boundary_polygon.gpkg") # Границы районов
stations = read_sf("data/metro_stations.gpkg") # Станции метро
water = read_sf("data/water_polygon.gpkg") # Водные объекты

# Прочитаем текущие параметры компоновки
def = par(no.readonly = TRUE)

# Уберем поля, чтобы карта занимала весь экран
par(mar = c(0,0,0,0))

# Получим ограничивающий прямоугольник слоя дорог в качестве общего охвата карты
frame = roads |> st_bbox() |> st_as_sfc() |> st_geometry()

## ОБЗОР ИСХОДНЫХ ДАННЫХ -------------------------------------

# Визуализируем входные данные

basemap = function(add = FALSE) {
  mf_base(frame, col = NA, add = add)
  mf_base(water, 
       col = "lightskyblue1",
       border = "lightskyblue3",
       add = TRUE)
  mf_base(roads,
       col = "gray70", 
       add = TRUE)
}

basemap()
mf_base(poi, 
     col = "deepskyblue4", 
     pch = 20, 
     cex = 0.5, 
     add = TRUE)
```

Теперь приступим к изучению данных, хранящихся в слое `poi` (от англ. POI --- Point Of Interest). Данный слой содержит все точечные маркеры OSM, которыми были отмечены на карте объекты, представляющие (по мнению создателей данных) интерес для пользователей. В POI включаются самые разнообразные объекты, такие как: объекты сферы услуг (amenity), места для отдыха (leisure), офисные здания (office), магазины и торговые центры (shop), туристические достопримечательности (tourism), спортивные объекты (sport), примечательные инженерные сооружения (man_made). В наших данных информация разнесена по соответствующим полям, каждый объект снабжен уникальным идентификатором:

```{r, echo = FALSE, purl = FALSE, paged.print = TRUE}
poi
```

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

```{r, collapse = TRUE, paged.print = TRUE}
data.frame(table(poi$AMENITY))
```

Для дальнейшего анализа отберем из всего множества объектов сферы услуг заведения, где можно поесть: рестораны, кафе, бары, пабы и заведения быстрого питания (фастфуд). В классификаторе OSM эти заведения имеют тип *restaurant*, *bar*, *cafe*, *pub* и *fast_food*. Для отбора нужных строк и столбцов используем dplyr:

```{r, collapse=TRUE}
poi_food = poi |>  
  select(NAME, AMENITY) |> 
  filter(AMENITY %in% c("restaurant", "bar", "cafe", 
                        "pub", "fast_food"))
head(poi_food)
```

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

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

1.  Вычислить матрицу расстояний между пунктами обслуживания и улицами. Размер матрицы $M \times N$, где $M$ --- количество улиц (строк), $N$ --- количество пунктов (столбцов)
2.  Найти в каждом столбце минимальное расстояние.
3.  Получить идентификатор улицы (номер строки), соответствующий данному расстоянию.
4.  Записать идентификатор в выходной вектор.

Таким образом, мы получим вектор из идентификаторов улиц, при этом каждый идентификатор будет встречаться в этом векторе столько раз, сколько раз данная улица оказалась ближайшей к какому-то объекту.

Вычислим матрицу расстояний с помощью функции `st_distance()` из пакета **sf**:

```{r, collapse=TRUE}

## АНАЛИЗ РАССТОЯНИЙ -------------------------------------

dist_matrix = st_distance(roads, poi_food)

# посмотрим, как выглядит результат на примере первых пяти объектов
print(dist_matrix[1:5,1:5])
```

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

```{r, collapse=TRUE}
ids = apply(dist_matrix, 2, function(X) order(X)[1])
```

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

```{r, collapse=TRUE}
count_stats = ids |> 
  table() |> 
  as_tibble() |> 
  mutate(ids = as.integer(ids))
roads = roads |> mutate(id = row_number())
roads_poi = left_join(roads, 
                   count_stats, 
                   by = c('id' = 'ids'))
```

Посмотрим первые 10 улиц по количеству общепита:

```{r, collapse=TRUE}
# Статистика по улицам в табличном представлении (первые 10)
roads_poi |> 
  select(NAME, n) |> 
  arrange(desc(n)) |> 
  head(10)
```

Для завершения анализа осталось визуализировать результаты:

```{r, collapse=TRUE}
basemap()
mf_map(roads_poi,
       type = 'prop',
       var = c('n'),
       col = 'red')
mf_base(poi_food, 
     col = "deepskyblue4", 
     pch = 20, cex = 1, 
     add = TRUE)
```

## Анализ взаимного положения (топологический) {#topology_analysis}

Пространственные запросы, основанные на топологических отношениях, позволяют находить объекты, находящиеся внутри других объектов, соприкасающиеся с другими объектами, пересекающиеся с ними и так далее. Топологические отношения сохраняются при взаимно-однозначных и непрерывных преобразованиях плоскости.

Отличия от метрических отношений легко пояснить на примере преобразования проекции. Представьте, что карту России в [конической проекции](http://geocnt.geonet.ru/projections/html/viewer.html) с концентрическими параллелями (известную по учебникам и атласам) вы трансформировали в карту России в проекции Меркатора (такую же как на [Google Maps](https://www.google.ru/maps/place/Россия/)). Изогнутые параллели превратились в прямые линии; форма регионов, площади и расстояния между населенными пунктами значительно изменились. Однако Красноярск по-прежнему находится в Красноярском крае, Ярославль --- на реке Волге, Нижний Новгород --- на правом берегу Волги, озеро Белое --- внутри Вологодской области, а Московская область как не граничила с Тамбовской, так и не граничит после трансформации проекции. Это и есть топологические отношения.

Формально топологические отношения в ГИС описываются с помощью [модели девяти пересечений DE-9IM](https://en.wikipedia.org/wiki/DE-9IM), которая была рассмотрена в начале этой лекции.

```{r}
## АНАЛИЗ ВЗАИМНОГО ПОЛОЖЕНИЯ -------------------------------------
poi_food = poi_food |> mutate(count = 1)
rayons_poi = aggregate(poi_food['count'], rayons, sum)
```

```{r, collapse=TRUE}
# Преобразуем результат в относительный показатель
# (единиц на кв.км. площади) и запишем в таблицу районов:
rayons_poi$density = 1000000 * rayons_poi$count / st_area(rayons_poi)
```

Масштабный множитель 1000000 в коде понадобился чтобы перевести площадь, хранящуюся в поле Shape_Area из квадратных метров в квадратные километры. Обратите внимание на то, что в данном случае мы не стали ограничивать фигурными скобками тело анонимной функции (`table(X)[2]`) внутри `apply()`, поскольку выполняемая операция достаточно компактна.

Визуализируем результат:

```{r, collapse=TRUE}
mf_map(rayons_poi,
       var = 'density',
       type = 'choro',
       leg_title = "Заведений\nна 1 кв.км")
basemap(add = TRUE)
mf_base(rayons,
     border = "black",
     lwd = 3, col = NA,
     add = TRUE)
mf_label(rayons,
     var = 'NAME',
     col = 'black',
     overlap = FALSE,
     font = 2,
     cex = 1)
```

Итак, используя топологический пространственный запрос "Содержит", мы смогли агрегировать точечные объекты внутри площадных и построить картограммы плотности распределения пунктов питания по районам центра Москвы.

## Анализ абсолютных зон окружения {#absolute_zones}

Задача данного раздела модуля звучит следующим образом: определить, какие пункты питания находятся в радиусе 300 метров от метро "Кропоткинская". Контекстом анализа в данном случае служит 300-метровая зона окружения станции метро. Поставленную задачу можно решить двумя способами:

-   Рассчитать расстояния от каждого пункта питания до станции метро "Кропоткинская" и выбрать точки, для которых это расстояние меньше или равно 300 метрам.
-   Построить буферную зону радиусом 300 метров и выбрать ею точки, используя топологическое отношение пересечения

Мы будем использовать второй вариант решения. Алгоритм выглядит следующим образом:

1.  Построить буферную зону, используя функцию `st_buffer()` из пакета **sf**.
2.  Выбрать полученной зоной точки пунктов питания, используя стандартный оператор `[]`.
3.  Визуализировать на карте полученные точки и буферную зону.

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

```{r, collapse=TRUE}
## АНАЛИЗ АБСОЛЮТНЫХ ЗОН ОКРУЖЕНИЯ -------------------------------------

stations$label = 'М' # пригодится для подписей

# Функция отвечает за рисование базовой карты
basemap2 = function(add = FALSE){
  
  mf_base(frame, col = NA, add = add)
  mf_base(water, 
       col = "lightskyblue1",
       border = "lightskyblue3",
       add = TRUE)
  mf_base(roads,
       col = "gray70", 
       add = TRUE)
  
  mf_base(poi_food, 
     col = "deepskyblue4", 
     pch = 20, cex = 0.5, 
     add = TRUE)
  
  mf_base(stations, 
     col = "slategray4", 
     pch = 20, 
     cex = 3, 
     add = TRUE)
  
  mf_label(stations,
     var = 'label',
     col = "white",
     cex = 0.6)
}
```

Определив вспомогательные функции, можем приступать к выполнению анализа:

```{r, collapse=TRUE}
# Выберем станцию метро и построим буферную зону
krop = filter(stations, NAME == "Кропоткинская")
zone = st_buffer(krop, dist = 500)

# Применим разработанную функцию для отбора точек
selected_poi = poi_food[zone, ]

# Применим разработанную функцию для рисования картографической основы
basemap2()

# Визуализируем результаты анализа
mf_base(zone,
     col = adjustcolor("sienna3", alpha.f = 0.5),
     border = "sienna3",
     add = TRUE)

mf_base(selected_poi, 
     col = "sienna4", 
     pch = 20, 
     cex = 1, 
     add = TRUE)

mf_base(krop, 
   col = "red", 
   pch = 20, 
   cex = 4, 
   add = TRUE)

mf_label(krop,
   var = 'label',
   col = "white",
   cex = 0.7)
```

```{r, echo = FALSE, purl = FALSE, paged.print = TRUE}
# Найденные объекты в табличном представлении:
selected_poi
```

В качестве примера аналогичного анализа отберем все пункты питания, находящиеся в пределах 100 метров от реки Москвы:

```{r, collapse=TRUE}
river = water |> filter(NAME == "Москва") |> st_union()
zone = st_buffer(river, dist = 100)

selected_poi = poi_food[zone, ]

basemap2()

# Визуализируем результаты анализа
mf_base(zone,
     col = adjustcolor("sienna3", alpha.f = 0.5),
     border = "sienna3",
     add = TRUE)

mf_base(river,
     col = adjustcolor("cyan", alpha.f = 0.5),
     border = "black",
     add = TRUE)

mf_base(selected_poi, 
     col = "red", 
     pch = 20, 
     cex = 1.5, 
     add = TRUE)
```

```{r, echo = FALSE, purl = FALSE, paged.print = TRUE}
# Найденные объекты в табличном представлении:
selected_poi
```

## Анализ конкурентных зон окружения {#conc_zones}

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

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

```{r, collapse=TRUE}

## АНАЛИЗ КОНКУРЕНТНЫХ ЗОН ОКРУЖЕНИЯ -------------------------------------

zones = stations |> 
  st_combine() |> 
  st_voronoi() |> 
  st_collection_extract() |> 
  st_crop(frame)
  
mf_base(zones)
mf_base(stations, add = TRUE, pch = 19, col = 'black')
```

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

```{r, collapse = TRUE}
# Агрегруем данные по каждой зоне
zones_poi = aggregate(poi_food['count'], zones, sum)

# Визуализируем результат

basemap2()
mf_base(zones, col = adjustcolor("white", alpha.f = 0.5),
        add = TRUE)
mf_map(zones_poi, var = 'count', type = 'prop', 
       col = adjustcolor("turquoise3", alpha.f = 0.5),
       border = F,
       leg_title = "Заведений\nпитания")
mf_label(zones_poi, 
     var = 'count',
     col = "turquoise4",
     cex = log(zones_poi$count)/4)
```

## Интерполяция, взвешенная на площадь {#areal_interpolation}

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

-   метод анализа (например, моделирование диффузии) предполагает, что данные распределены по регулярной сетке, в то время как исходная сетка нерегулярна.
-   необходимо обеспечить сравнимость пространственных распределений показателя для разных территорий, в то время как дробность исходного территориального деления существенно меняется в пространстве.

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

--- *экстенсивные* параметры суммируются и делятся при агрегировании/агрегировании территориальных единиц. Например, площадь, покрытая лесом или численность населения --- это экстенсивный параметр. - *интенсивные* параметры осредняются или остаются постоянными при агрегировании/дизагрегировании территориальных единиц. Например, густота древостоя и плотность населения --- интенсивные параметры.

Рассмотрим это метод интерполяции на примере данных по графствам Северной Каролины (показатель --- количество новорожденных в 1974 году). Для расчета векторной регулярной сетки используем функцию `st_make_grid()` из пакета **sf**.

```{r}
# Данные по Северной Каролине
nc = st_read(system.file("shape/nc.shp", package="sf"))

cells = st_make_grid(nc, cellsize = 0.25)

birth = st_interpolate_aw(nc["BIR74"], 
                       cells, 
                       extensive = FALSE)

# исходное распределение
mf_map(nc, var = 'BIR74', type = 'choro')

# пересчет на регулярную сетку
mf_map(birth, var = 'BIR74', type = 'choro')
mf_base(nc, col = NA, border = 'white', add = TRUE)
```

## Дирекционные отношения

### Статистика направлений {#circular_circ}

#### Теория {#circular_circ_theory}

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

Обработкой данных о направлениях занимается особая область математической статистики --- **статистика направлений**, или **круговая (циркулярная)** статистика (Mardia, Jupp, 2000; Pewsey et al., 2013). В круговой статистике каждое направление $\theta \in [0, 2\pi)$ представляется в виде вектора $x = (\cos \theta, sin \theta)$. Все операции производятся над подобными векторами и их координатами. Аналогом нормального распределения для круговой случайной величины является распределение фон Мизеса (von Mises, 1918), которое задается функцией плотности вероятности: $$
f(θ)=\frac{1}{2 \pi I_0(\kappa)} e^{\kappa \cos (\theta - \mu)},
$$

где $\kappa \geq 0$ --- параметр концентрации, $\mu$ --- среднее значение (для $\kappa > 0$) и

$$
I_p(\kappa) = \frac{1}{2π} \int_{0}^{2\pi} \cos (p \theta) e^{\kappa \cos θ} d \theta
$$ есть модифицированная функция Бесселя первого рода и порядка $p$. Из формул видно, что по своему эффекту параметр концентрации противоположен среднеквадратическому отклонению $\sigma$, которое является параметром нормального распределения. Чем больше значение $\kappa$, тем более сконцентрировано распределение относительно среднего значения --- отсюда идет название этого параметра. Распределение фон Мизеса используется для построения ядра при аппроксимации плотности распределения направлений методом [ядерной оценки](https://en.wikipedia.org/wiki/Kernel_density_estimation) (оценки по методу Парзена-Розенблатта).

> В метеорологии значения $\cos \theta$ и $\sin \theta$ определяют соотношение **зональной** и **меридиональной** составляющей скорости \[ветра\] (для получения самих составляющих их надо умножить на скорость ветра).

Для вычисления статистических моментов круговой случайной величины требуется найти средний равнодействующий вектор первого порядка: $$R = (C, S),$$ где

$$C = \frac{1}{n} \sum_{j=1}^{n} \cos \theta_j,\\ 
S = \frac{1}{n} \sum_{j=1}^{n} \sin \theta_j.$$

Данный вектор имеет направление $\bar\theta$, которое является **выборочным средним направлением** исследуемой величины.

**Выборочная средняя равнодействующая длина** $\bar R = \sqrt{C^2 + S^2}$ принимает значения в диапазоне $[0, 1]$ и показывает меру концентрации направлений относительно $\theta$. $\bar R = 1$ означает, что все исходные направления совпадают, $\bar R = 0$ --- что данные равномерно распределены по кругу, либо распределение имеет несколько мод, которые уравновешивают друг друга.

Величина $\bar R$ дает важную информацию для предварительной диагностики картины направлений. Если значение $\bar R$ близко к единице, это означает, что распределение является унимодальным и в качестве основного направления можно принять значение $\bar θ$ [@mardia2000directional].

**Стандартное отклонение направлений** $v$ в радианах может быть найдено как $v=\sqrt{-2 \ln \bar R}$ .

В ряде случаев противоположные направления считаются эквивалентными. Например, нельзя сказать, идет ли улица с юга на север или с севера на юг. Такие данные в теории круговой статистики называются **аксиальными** (Mardia, Jupp, 2000). Для аксиальных данных возможный диапазон значений лежит в интервале $[0, \pi)$. Поскольку методы круговой статистики рассчитаны на круговое замыкание данных, стандартный подход к обработке аксиальных данных предполагает переход от направлений к их удвоенным значениям $\theta' = 2\theta$, обработку полученных значений стандартными методами и отображение полученных значение обратно на интервал $[0, \pi)$. Для среднего, медианы и моды распределения это означает простое деление полученного значения пополам [@pewsey2013circular].

Модальные направления могут быть определены как по гистограмме распределения, так и методом ядерной оценки. Основной вопрос поиска эффективного ядра заключается в параметризации функции $K$. Для распределения фон Мизеса таким параметром является концентрация $\kappa$. Чем больше этот параметр, тем более локализованной будет оценка, тем сильнее будут проявляться в ней существующие моды распределения, но также будут и выделяться новые моды, которые на самом деле не значимы. Малые значения $\kappa$ приведут, наоборот, к «размыванию» плотности распределения в пределах полного круга. Как и в случае с количеством интервалов гистограммы, избыточно малые и большие значения κ нежелательны.

В работе [@Oliveira2012] показано, что оптимальное значение $\kappa$ может быть подобрано также для оценки распределений, являющихся конечной суммой $M$ распределений фон Мизеса, то есть, мультимодальных распределений, имеющих плотность : $$g(\theta)=\sum_{i=1}^{M} \alpha_i \frac{\exp\lbrace{\kappa_i \cos(\theta - \mu_i)\rbrace}}{2 \pi I_0 (\kappa_i)},$$ где $\sum_{i=1}^{M} = 1$.

Поскольку в результате подбора определяется не только параметр концентрации, но и число компонент в сумме распределений [@Oliveira2014], его можно также использовать для определения количества искомых мод, если это необходимо.

Когда подобрана функция ядра и ее параметры, оценка плотности распределения (вычисление функции $\circ f _h (x)$) для круговых данных делается либо для исходных направлений $\theta_j$, либо с равным (достаточно малым) интервалом --- например, через 1 градус [@pewsey2013circular]. После того как произведена оценка, могут быть выбраны направления, в которых функция плотности распределения достигает локального максимума --- первого и второго по величине. Эти направления и будут соответствовать первой и второй моде распределения направлений.

### Практика {#circular_circ_vis}

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

В системе доступны данные по следующим обсерваториям:

```{r}
obs = read_excel('data/bound/scheme.xlsx', 2)
```

```{r, echo = FALSE}
knitr::kable(obs)
```

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

```{r, echo = FALSE}
params = read_excel('data/bound/scheme.xlsx', 1)
```

```{r, echo = FALSE}
knitr::kable(params)
```

Загрузим данные по всем обсерваториям из текстовых файлов в папке *bound*:

```{r}
files = paste('data/bound', list.files('data/bound', "*.txt"), sep = '/')

(tab = lapply(files, function(X) {
    read_table(X, col_names = params$Обозначение)
  }) |> 
  bind_rows() |> 
  left_join(obs, by = c('INDEX' = 'Индекс'))) # присоединим информацию о названиях станций
```

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

```{r}
geo_azimuth = function(dx, dy) {
  a = atan2(dx, dy)
  ifelse(a <= pi/2, pi/2 - a,  5*pi/2 - a)
}

(winds = tab |> 
  mutate(wind = circular(geo_azimuth(MV, MU), template = 'geographics')) %>% 
  select(INDEX, name = Название, GGGG, MM, HH, Z, MU, MV, SS, wind))
```

Выберем данные по высоте 0 метров за 12 часов дня для поселка Тикси, сохранив только составляющие скорости и ее скалярную величину:

```{r}
(tiksi_wind = winds |> filter(name == 'Тикси', HH == 12, Z == 0))
```

Отобразим распределение направлений, розу-диаграмму и плотность распределения. Для построени графиков используем функции `plot.circular()` и `rose.diag` из пакета **circular**. Для аппроксимации плотности распределения направлений воспользуемся функцией `kern.den.circ()` из пакета **NPCirc**. Эта функция использует функцию плотности распределения *фон Мизеса* в качестве ядра и по умолчанию разбивает круг на 250 направлений, по которым производится оценка плотности (при необходимости это значение можно изменить в параметре `len`):

```{r}
plot.circular(tiksi_wind$wind, 
     cex = 0.5, 
     stack = TRUE, 
     sep = 0.035,
     axes = FALSE,
     main = 'Среднемноголетняя роза ветров в Тикси',
     sub = 'Измерения за период с 2007 по 2018 г, высота 0 м')

rose.diag(tiksi_wind$wind, 
          bins = 8, 
          col = 'gray70',
          border = 'gray30',
          prop = 1, 
          add = TRUE, 
          tick = FALSE,
          lwd = 0.5)

kden = kern.den.circ(tiksi_wind$wind)

lines(kden, shrink = 3, # параметр shrink отвечает за масштаб радиус-вектора
      join = F,
      col = 'steelblue')
```

> Параметр `shrink` отвечает за масштаб радиус-вектора на графиках из пакета **circular**. Чем больше его величина, тем сильнее будет сжат график относительно центра круга.

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

```{r}
peak = findpeaks(kden$y, sortstr = T)[1,2] # находим индекс самого высокого пика плотности распределения

(modal = kden$x[peak]) # извлекаем сам угол

# раскладываем на составляющие для отрисовки линии
xp = sin(modal)
yp = cos(modal)

plot.circular(tiksi_wind$wind, 
     cex = 0.5, 
     stack = TRUE, 
     sep = 0.035,
     axes = FALSE,
     main = 'Среднемноголетняя роза ветров в Тикси',
     sub = 'Измерения за период с 2007 по 2018 г, высота 0 м')

rose.diag(tiksi_wind$wind, 
          bins = 8, 
          col = 'gray70',
          border = 'gray30',
          prop = 1, 
          add = TRUE, 
          tick = FALSE,
          lwd = 0.5)

lines(kden, shrink = 3, 
      join = F, col = 'steelblue')

lines(c(0, xp), c(0, yp),
      lwd = 2, col = 'orangered')

text(x = 1.4 * xp, y = 1.4 * yp, 
     col = 'orangered',
     labels = paste0(round(180 * modal / pi, 0), '°')) # приводим к целым градусам
```

Проведем анализ направлений для всех станций. Для этого рассчитаем функции плотности распределения и разместим их в новом фрейме данных с лист-колонкой.

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

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

```{r}
(dens = winds |> 
  filter(HH == 12, Z == 0) |> 
  group_by(name) |> 
  summarise(kden = list(kern.den.circ(wind))) |> 
  mutate(peak = sapply(kden, function(X) {
                  peak = findpeaks(X$y, sortstr = T)[1,2]
                  X$x[peak]
                })
  )
)
```

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

```{r}
# устанавливаем параметры компоновки
par(mar = c(1,1,1,1),
    mfrow = c(1,2))

# строим графики в цикле
for (obs_name in dens$name) {
  
  wind_df = winds |> filter(name == obs_name, HH == 12, Z == 0)
  dens_df = dens |> filter(name == obs_name)
  
  modal = dens_df$peak

  xp = sin(modal)
  yp = cos(modal)
  
  plot.circular(wind_df$wind, 
       shrink = 1.2,
       cex = 0.5, 
       stack = TRUE, 
       sep = 0.035,
       axes = FALSE,
       main = obs_name)
  
  rose.diag(wind_df$wind, 
            bins = 8, 
            col = 'gray70',
            border = 'gray30',
            prop = 1, 
            add = TRUE, 
            tick = FALSE,
            lwd = 0.5)
  
  lines(dens_df$kden[[1]], 
        shrink = 3, join=F,
        col = 'steelblue')
  
  lines(c(0, xp), c(0, yp),
        lwd = 2, col = 'orangered')
  
  text(x = 1.4 * xp, y = 1.4 * yp, 
       col = 'orangered',
       labels = paste0(round(180 * modal / pi, 0), '°')) # приводим к целым градусам
}
```

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

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

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

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

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

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

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

1.  Перечислите три основных вида пространственных отношений, приведите их примеры.
2.  Перечислите 8 вариантов топологических отношений и названий функций **sf**, которые им соответствуют.
3.  Опишите способ, с помощью которого можно выбрать пространственные объекты, пересекающиеся с заданным множеством пространственных объектов.
4.  Каким образом можно заменить тип топологического отношения с пересечения на любой другой при выполнении пространственной выборки?
5.  Чем отличаются контекстные и целевые объекты?
6.  В чем заключается отличие абсолютных и конкурентных зон окружения?
7.  Какая функция пакета sf позволяет вычислять расстояния между объектами? Как с помощью полученного результата определить для каждого объекта из множества $A$ определить ближайший к нему объект из множества $B$?
8.  Опишите последовательность действий, которую необходимо выполнить для подсчета количества точечных объектов по заданной сетке полигонов.
9.  С помощью какой функции можно построить буферную зону вокруг пространственного объекта? Есть ли ограничения на размерность буферизуемого пространственного объекта (точка, линия, полигон)? Можно ли построить буфер вокруг поверхности?
10. Какая геометрическая структура используется для построения конкурентных зон окружения?
11. Что такое OSRM?
12. Какими средствами можно построить зоны транспортной доступности и маршруты в среде R? В какой системе координат должны быть точки, участвующие в сетевом анализе?
13. Опишите возможности и основные функции пакета `cartography`, с помощью которых можно строить тематические карты способами картограмм, картодиаграмм и линейных знаков, а также легенды к ним.
14. В каком виде направления рассматриваются в круговой статистике?
15. Какое распределение является аналогом нормального распределения для круговых данных? Что означают параметры $\kappa$ и $\mu$ в функции этого распределения?
16. Как вычисляется равнодействующий вектор первого порядка и выборочная средняя равнодействующая длина этого вектора для направлений?
17. Чем аксиальные данные отличаются от круговых данных в общем случае? Какие преобразования осуществляются над такими данными для того чтобы применять к ним стандартные методы циркулярной статистики?
18. Какой класс данных (и пакет) можно использовать в R для представления направлений? Как указать, что направления отсчитываются географическим методом, то есть, по часовой стрелки от направления на север?
19. Какую функцию можно использовать для для оценки плотности распределения круговых данных? В каком пакете она находится?
20. Какую функцию можно использовать для выявления модальных направлений по данным функции плотности вероятности?
21. Какие функции позволяют строить диаграммы и розы-диаграммы по круговым данным в среде R?
22. Какой параметр управляет масштабом радиус-вектора на круговых графиках?
23. Что такое лист-колонка в фрейме данных, и какого типа данные можно в ней хранить?

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

1.  Проанализируйте пространственную ассоциацию подтипов почв с типами рельефа [данным](https://github.com/tsamsonov/r-geo-course/blob/master/data/Satino.gpkg) ГИС Сатино. Для этого выполните оверлей между слоями *RelTypes* и *SoilTypes* методом [`st_intersection()`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html). Для каждого подтипа почв рассчитайте долю, которая занята в его площади каждым типом рельефа. Визуализируйте результаты средствами **ggplot2** в виде столбчатой диаграммы, где каждый столбик отвечает за подтип почвы, а его внутреннее разделение соответствует долям типов рельефа. Используя функцию [`cramerV()`](https://www.rdocumentation.org/packages/rcompanion/versions/2.3.25/topics/cramerV) из пакета **rcompanion**, рассчитайте [коэффициент ассоциации Крамера](https://en.wikipedia.org/wiki/Cram%C3%A9r%27s_V), чтобы охарактеризовать силу ассоциации между этими номинальными переменными.

    > **Подсказка**: для вычисления коэффциента Крамера вам необходимо преобразовать данные в широкую форму, где подтипы почв идут по строкам, а типы рельефа --- по столбцам. Полученную таблицу необходимо конвертировать в матрицу и подать на вход функции `cramerV()`.

2.  Одна из гипотез, часто используемых в [геомаркетинге](https://ru.wikipedia.org/wiki/%D0%93%D0%B5%D0%BE%D0%BC%D0%B0%D1%80%D0%BA%D0%B5%D1%82%D0%B8%D0%BD%D0%B3) --- это так называемые *аттракторы потоков* --- пространственные объекты, которые сосредотачивают в своей близости высокую плотность пешеходного трафика. Типичный пример аттрактора --- любая транспортная локация: выход из метро, железнодорожная платформа, автобусная остановка. Владельцы предприятий сферы услуг в теории стремятся размещать свои точки вблизи к аттракторам. Используя данные из настоящей лекции, проведите проверку реалистичности этой теории. Для этого:

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

    Далее постройте график зависимости между радиусом буфера и плотностью объектов интереса. Рассчитайте также коэффициент корреляции между этими величинами.

3.  Субъекты Российской Федерации пронумерованы числами от 1 до 85 в алфавитном порядке. Для решения задач пространственного анализа часто бывает необходимо, чтобы нумерация была пространственная. В этом случае соседние субъекты будут иметь похожие номера. Придумайте методику такой нумерации и реализуйте ее в виде программы на языке R.

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