class: center, middle, inverse, title-slide .title[ # Векторный анализ ] .subtitle[ ## Визуализация и анализ географических данных на языке R ] .author[ ### Тимофей Самсонов ] .institute[ ### МГУ имени Ломоносова, Географический факультет ] .date[ ### 2023-11-28 ] --- ## Пакеты Для выполнения кода данной лекции вам понадобятся следующие пакеты: ```r library(sf) library(dplyr) library(classInt) library(mapsf) # Удобное построение тематических карт средствами plot() library(circular) # статистика направлений library(NPCirc) # статистика направлений library(pracma) options(scipen = 999) ``` --- ## Пакет mapsf <iframe src="https://riatelab.github.io/mapsf/articles/web_only/img/mapsf_cheatsheet.png" width="100%" height="400px" data-external="1"></iframe> https://raw.githubusercontent.com/riatelab/mapsf/master/vignettes/web_only/img/mapsf_cheatsheet.pdf --- ## Пространственные запросы __Пространственные запросы__ связаны с поиском объектов (географических локаций), удовлетворяющих условию, заданному на множестве пространственных отношений. Они бывают трех типов: - _Метрические_ (расстояния) - _Топологические_ (взаимное размещение) - _Дирекционные_ (направления) Примеры пространственных запросов знакомы любому географу: - Найти все объекты не далее 100 метров от дороги (метрические отношения) - Найти все объекты внутри границы заповедника (топологические отношения) - Найти все объекты в северо-восточной четверти относительно точки (дирекционные отношения) --- ## Метрические отношения .code-small[ ```r # Чтение данных roads = read_sf("../r-geo-course/data/roads.gpkg") # Дороги rayons = read_sf("../r-geo-course/data/boundary_polygon.gpkg") # Границы районов stations = read_sf("../r-geo-course/data/metro_stations.gpkg") # Станции метро water = read_sf("../r-geo-course/data/water_polygon.gpkg") # Водные объекты poi = read_sf("../r-geo-course/data/poi_point.gpkg") # Точки интереса poi_food = poi |> # Только общепит select(NAME, AMENITY) |> filter(AMENITY %in% c("restaurant", "bar", "cafe", "pub", "fast_food")) # Прочитаем текущие параметры компоновки def = par(no.readonly = TRUE) # Уберем поля, чтобы карта занимала весь экран par(mar = c(0,0,0,0)) # Получим ограничивающий прямоугольник слоя дорог в качестве общего охвата карты frame = roads |> st_bbox() |> st_as_sfc() |> st_geometry() ``` ] --- ## Метрические отношения .pull-left[ __Задача__: подсчитать сколько заведений общепита на каждой улице. ```r # Визуализируем входные данные plot(frame) plot(water, col = "lightskyblue1", border = "lightskyblue3", add = TRUE) plot(roads, col = "gray70", add = TRUE) plot(poi_food, col = "deepskyblue4", pch = 20, cex = 0.5, add = TRUE) ``` ] .pull-right[ <!-- --> ] --- ## Метрические отношения Для этого вычислим матрицу попарных расстояний между дорогами и точками и найдем ближайшую улицу к каждой точке: ```r dist_matrix = st_distance(roads, poi_food) print(dist_matrix[1:5,1:5]) ## Units: [m] ## [,1] [,2] [,3] [,4] [,5] ## [1,] 4962.292 3420.6849 3240.2066 2145.2044 4748.686 ## [2,] 2247.737 705.2923 524.9411 570.3986 2035.341 ## [3,] 2213.236 670.7904 490.4167 605.1606 2000.759 ## [4,] 2197.874 655.4285 475.0629 620.4411 1985.242 ## [5,] 3910.957 2368.5560 2188.1348 1092.5472 3698.246 id = apply(dist_matrix, 2, \(X) order(X)[1]) head(id) ## [1] 13 19 12 29 39 23 ``` --- ## Метрические отношения .pull-left[ Результаты необходимо проагрегировать: ```r stats = table(id) |> as_tibble() |> mutate(id = as.integer(id)) print(stats, n = 5) ## # A tibble: 657 × 2 ## id n ## <int> <int> ## 1 6 1 ## 2 7 1 ## 3 10 2 ## 4 12 4 ## 5 13 2 ## # ℹ 652 more rows ``` ] .pull-right[ И далее присоединить к линиям дорог: ```r roads_stats = roads |> mutate(id = row_number()) |> left_join(stats, by = 'id') |> mutate(n = ifelse(is.na(n), 0, n)) |> arrange(n) ``` ```r roads_stats |> select(NAME, n) |> tail(5) |> st_drop_geometry() ## # A tibble: 5 × 2 ## NAME n ## * <chr> <dbl> ## 1 Пятницкая улица 17 ## 2 Мясницкая улица 20 ## 3 Цветной бульвар 22 ## 4 улица Новый Арбат 24 ## 5 улица Арбат 43 ``` ] --- ## Метрические отношения .pull-left[ Визуализируем результаты с помощью __mapsf__: ```r mf_init(frame) mf_map(water, col = 'steelblue3', add = TRUE) mf_map(roads_stats, var = 'n', type = 'choro', breaks = 'fisher', nbreaks = 7, pal = 'Reds', lwd = 4, col_na = 'grey', add = TRUE) mf_map(poi_food, cex = 0.5, col = "deepskyblue4", add = TRUE) ``` ] .pull-right[ <!-- --> ] --- ## Топологические отношения <img src="13_Vector_files/figure-html/unnamed-chunk-10-1.png" height="450px" /> --- ## Топологические отношения Функция | Топологическое отношение ------------------------------|-------------------------------------------------------------------- `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` совпадают) --- ## Топологические отношения Выборка по пересечению: .pull-left[ .code-small[ ```r countries = st_read( '../r-geo-course/data/ne/countries.gpkg', quiet = T ) outlines = st_geometry(countries) cities = st_read( '../r-geo-course/data/ne/cities.gpkg', quiet = T ) # Наносим исходную конфигурацию plot(outlines, lwd = 0.5) plot(cities, col = 'black', pch = 20, cex = 0.5, add = TRUE) ``` ] ] .pull-right[ <!-- --> ] --- ## Топологические отношения Выборка по пересечению: .pull-left[ .code-small[ ```r # отключаем s2 чтобы не было ошибок # при выборке sf::sf_use_s2(FALSE) largest = countries |> select(pop_est) |> filter(pop_est > 100000000) # Отбираем точки внутри стран # с максимальным ВВП *sel = cities[largest, ] # Смотрим что получилось plot(outlines, lwd = 0.5) plot(largest, col = 'gray', add = TRUE) plot(sel, pch = 20, col = 'black', add = TRUE) ``` ] ] .pull-right[ <!-- --> ] --- ## Топологические отношения Выборка по касанию границ: .pull-left[ .code-small[ ```r cz = countries |> filter(sovereignt == 'Czechia') *neighbors = countries[cz, op=st_touches] plot(st_geometry(neighbors), col = 'lightgray', lwd = 0.5) plot(cz, col = 'darkgray', add = TRUE) ``` ] ] .pull-right[ <!-- --> ] --- ## Топологические отношения .pull-left[ Количественное агрегирование по районам: .code-small[ ```r poi_food = poi_food |> mutate(count = 1) rayons_poi = aggregate( poi_food['count'], rayons, sum ) %>% mutate(density = 1000000 * count / st_area(.)) mf_init(frame) mf_map(rayons_poi, var = 'density', lwd = 4, type = 'choro', breaks = 'pretty', nbreaks = 7, pal = 'Greens', add = TRUE) mf_map(water, col = 'steelblue3', add = TRUE) mf_map(roads, col = 'black', add = T) mf_map(poi_food, cex = 0.5, col = "deepskyblue4", add = TRUE) ``` ] ] .pull-right[ <!-- --> ] --- ## Абсолютные зоны окружения (буферные) .pull-left[ .code-small[ ```r # Выберем станцию метро и построим буферную зону krop = stations |> filter(NAME == "Кропоткинская") zone = st_buffer(krop, dist = 300) sel_poi = poi_food[zone, ] mf_init(frame) mf_map(water, col = 'steelblue3', add = TRUE) mf_map(zone, col = 'sienna1', lwd = 4, add = T) mf_map(roads, col = 'black', add = T) mf_map(poi_food, cex = 0.5, col = "deepskyblue4", add = TRUE) mf_map(sel_poi, cex = 1.5, col = "darkred", add = TRUE) ``` ] ] .pull-right[ <!-- --> ] --- ## Конкурентные зоны окружения (диаграмма Вороного) .pull-left[ ```r zones = stations |> as('Spatial') |> dismo::voronoi() |> st_as_sf() |> st_crop(frame) plot(zones %>% st_geometry()) plot(stations, add = TRUE, pch = 19, col = 'black') ``` ] .pull-right[ <!-- --> ] --- ## Конкурентные зоны окружения (диаграмма Вороного) .pull-left[ Агрегирование по ячейкам диаграммы Вороного: .code-small[ ```r zones_poi = aggregate(poi_food['count'], zones, sum) mf_init(frame) mf_map(water, col = 'steelblue3', add = TRUE) mf_map(zones_poi, lwd = 3, col = NULL, add = T) mf_map(zones_poi, var = 'count', col = 'olivedrab', type = 'prop', add = T) mf_map(roads, col = 'black', add = T) mf_map(poi_food, cex = 0.5, col = "deepskyblue4", add = TRUE) ``` ] ] .pull-right[ <!-- --> ] --- ## Дирекционные данные Среднемесячные значения метеопараметров в пограничном слое атмосферы по полярным аэрологическим обсерваториям России (ВНИИГМИ-МЦД) ```r (obs = readxl::read_excel('../r-geo-course/data/bound/scheme.xlsx', 2)) ## # A tibble: 13 × 4 ## Индекс Название Широта Долгота ## <dbl> <chr> <dbl> <dbl> ## 1 20674 Остров Диксон 73.5 80.4 ## 2 21824 Тикси 71.4 129. ## 3 22113 Мурманск 68.6 33.1 ## 4 22217 Кандалакша 67.1 32.2 ## 5 22271 Шойна 67.5 44.1 ## 6 23078 Норильск 69.2 88.2 ## 7 23205 Нарьян-Мар 67.4 53.1 ## 8 23330 Салехард 66.3 66.4 ## 9 24125 Оленек 68.3 112. ## 10 24266 Верхоянск 67.6 133. ## 11 24343 Жиганск 66.5 123. ## 12 89512 Новолазаревская -70.8 11.8 ## 13 89592 Мирный -66.6 19.7 ```