Карен Валитов, “Пространственные операции, применяемые при подготовки данных для анализа”

План митапа

1. Работа с геометриями
2. Система координат для городов
3. Геометрические операции
4. Пространственные операции
5. Пространственные соотношения
6. Пространственные преобразования
7. Сетки и симуляция карт плотности показателя

Секция Интро

Загрузим библиотеки:

pkgs <- c("sf", "tidyverse", "geojsonio", "geojsonsf", "readr", "mapview", "leaflet", "sfnetworks", "tidygraph", "data.table", "units", "openxlsx")

if (length(setdiff(pkgs, rownames(installed.packages()))) > 0) {
  
  install.packages(setdiff(pkgs, rownames(installed.packages())), quiet = T)
}

lapply(pkgs, library, character.only = T)

Загрузим датасеты. Для чтения векторных данных из основных форматов работы с пространственными данными (geojson, geopackage, shapefiles etc) есть функция ‘st_read’. Она загрузит ваш файл и представит его в виде объекта sf (который как мы помним также dataframe)

districts <- st_read('data/districts.geojson')

## Reading layer `districts' from data source 
##   `/builds/upravitelev/r_meetups/content/post/2022-05-21-eu_spb_meetup/data/districts.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 146 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 36.80313 ymin: 55.14222 xmax: 37.96747 ymax: 56.02125
## Geodetic CRS:  WGS 84

districts <- st_make_valid(districts)

houses <- st_read('data/buildings.shp')

## Reading layer `buildings' from data source 
##   `/builds/upravitelev/r_meetups/content/post/2022-05-21-eu_spb_meetup/data/buildings.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1035 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 37.58994 ymin: 55.78266 xmax: 37.63449 ymax: 55.81012
## Geodetic CRS:  WGS 84

houses <- st_make_valid(houses)

social_objects <- st_read('data/education_poi.shp')

## Reading layer `education_poi' from data source 
##   `/builds/upravitelev/r_meetups/content/post/2022-05-21-eu_spb_meetup/data/education_poi.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 15 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 37.60057 ymin: 55.78663 xmax: 37.62163 ymax: 55.80323
## Geodetic CRS:  WGS 84

districts <- st_make_valid(districts)

streets <- st_read('data/highway.shp')

## Reading layer `highway' from data source 
##   `/builds/upravitelev/r_meetups/content/post/2022-05-21-eu_spb_meetup/data/highway.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 20290 features and 13 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 362176.6 ymin: 6115267 xmax: 434921.9 ymax: 6209619
## Projected CRS: WGS 84 / UTM zone 37N

streets <- select(streets, OSM_ID, NAME, HIGHWAY)
districts <- st_make_valid(districts)

Первый - данные по районам города Москвы. Посмотрим кратко на эти данные

mapview(districts, zcol='name', legend=F)

Второй - данные расположения домов в районе Марьина Роща, где указано, жилой дом или нет.

mapview(houses, legend=F)

Третий сет - данные расположения школ и детских садов в районе Марьина Роща.

mapview(social_objects, zcol='type_ed', legend=F)

Посмотрим на все наборы вместе

ggplot()+
  geom_sf(data=districts,fill="#b8dae0")+
  geom_sf(data=houses,color="#f5f53d")+
  geom_sf(data=social_objects,color="red")+
  geom_sf(data=streets)+
  theme_minimal()

Комментарий Карена: Эти данные совершенно обычны для урбанистов, их можно найти в разных местах: OSM, Портал Открытых Данных Москвы, Реформа ЖКХ и другие

1. Топологические операции

Изменение типа геометрий

Вспоминаем про точки - линии - полигоны. Можем всегда разломать одни в объекты в другие, но важно помнить про структуру самих геометрий. Можно ли точки разложить на линии?

districts_center <- filter(districts, name %in% c("Арбат", "Басманный", "Замоскворечье", "Красносельский", "Мещанский", "Пресненский", "Таганский", "Тверской", "Хамовники", "Якиманка"))

polygon <- st_cast(districts_center,"POLYGON")

ggplot()+
  geom_sf(data=polygon)+
  ggtitle('Москва POLYGON')

polylines <- st_cast(polygon,"LINESTRING")

ggplot()+
  geom_sf(data=polylines)+
  ggtitle('Москва LINESTRING')

multipoint <-st_cast(polylines,"MULTIPOINT")

ggplot()+
  geom_sf(data=multipoint)+
 ggtitle('Москва MULTIPOINT')

# Ступенчатая логика - от более сложной геометрии к простой. Не советуем перескакивать через "уровни".

Moscow_simp <- list()

for (i in c(100, 300, 500)) {
  simplif_list <- list(st_simplify(st_transform(districts_center, 32637), dTolerance =  i))
  names(simplif_list) <- i
  Moscow_simp <- append(Moscow_simp,
                        list(simplif_list))
}

plot(Moscow_simp[[1]][["100"]][1], main = paste("Simplification =", names(Moscow_simp[[1]]), sep = " "))

plot(Moscow_simp[[2]][["300"]][1], main = paste("Simplification =", names(Moscow_simp[[2]]), sep = " "))

plot(Moscow_simp[[3]][["500"]][1], main = paste("Simplification =", names(Moscow_simp[[3]]), sep = " "))

2. Выбор проекции для городских проектов

Как подобрать рабочую проекцию?

Немного снова вспомним про проекции? При работе в городах можно использовать локальную проекцию UTM. Что это значит? Вспомним про проекцию Меркатор и как ее построили.

А теперь мы повернем цилиндр на бок и начнем его вращать!

Итого получили 60 зон (по 6 градусов), где нет почти нет никаких искажений. Подходит для городов, поскольку они обычно лежат в одной зоне (редко на границе двух, но не критично).

3. Геометрические вычисления и операции

Работая с геометриями, мы можем делать обычные геометрические вычисления. Это могут быть площади, длины и другие. Важно, что эти значения зависят от системы координат. Библиотека SF перешла на использования библиотеки S2 от Google и теперь в WGS84 ведет геодезические вычисления. У этого есть плюс и минусы.

districts <- st_transform(districts, crs = 32637)
districts <- mutate(districts, area=st_area(districts))
districts <- st_transform(districts, crs = 4326)
districts <- mutate(districts, area_s2=st_area(districts))

knitr::kable(head(districts))

centroid <- st_centroid(districts_center)

st_coordinates(centroid)

##              X        Y
##  [1,] 37.57483 55.72911
##  [2,] 37.58998 55.75116
##  [3,] 37.63375 55.73363
##  [4,] 37.62874 55.78002
##  [5,] 37.60729 55.73057
##  [6,] 37.66681 55.74079
##  [7,] 37.55942 55.76089
##  [8,] 37.65453 55.77787
##  [9,] 37.67080 55.76668
## [10,] 37.60663 55.77024

4. Пространственные операции

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

social_objects <- mutate(social_objects, dist = ifelse(type_ed == "школа", 500, 300))

social_objects <- st_transform(social_objects, crs = 32637)
social_buffer <- st_buffer(social_objects, social_objects$dist)
social_buffer <- st_transform(social_buffer, crs = 4326)

mapview(social_buffer, zcol='type_ed', legend=F)

Мы также можем резать объекты по геометриям друг друга, объединять их и многое другое. По ссылке представлены несколько операций, которые являются основными в рамках работы с геометриями (https://gisgeography.com/geoprocessing-tools/)

Наши дороги выходят за район исследования, давайте обрежем их

districts <- st_transform(districts, crs = 32637)
streets_mr <- st_intersection(streets, filter(districts, name == "Марьина Роща"))
districts <- st_transform(districts, crs = 4326)
streets_mr <- st_transform(streets_mr, crs = 4326)


mapview(streets_mr, zcol='NAME', legend=F)

5. Пространственные предикаты

Каждая пара слоев (и объектов к них) как-то пространственно соотносится друг с другом. Для ответа на вопрос как, есть набор пространственных предикатов. Каждый задает некоторый вопрос паре и в результате дает на ответ 0 или 1.

Внутри sf такие операции реализованы через следующие команды:

Посмотрим просто на карте, как соотносятся буферы социальных объектов и дома района

 ggplot()+
  geom_sf(data = filter(districts, name == "Марьина Роща"), fill="#b8dae0")+
  geom_sf(data=social_buffer,color="red")+
  geom_sf(data=houses,color="orange", fill = "green")+
  geom_sf(data=streets_mr)+
  theme_minimal()

Можем перейти к центроидам и проверить уже их соотношение с буферами.

houses_centroid <- st_centroid(houses)

ggplot()+
  geom_sf(data = filter(districts, name == "Марьина Роща"), fill="#b8dae0")+
  geom_sf(data=social_buffer,color="red")+
  geom_sf(data= houses_centroid, color="orange")+
  geom_sf(data=streets_mr)+
  theme_minimal()

# 6. Пространственное присоединение

Обычный джоин работает по ключевому полю или по какому-то условия. Пространственный предикат может стать таким условием.

7. Сетки и симуляция карт плотности показателя

Для расчета карт плотности можно использовать алгоритмы растрового анализа. Это - стандарт индустрии, которые много где используется. Но для решения это задачи можно также использовать сетки (растр по сети тоже сетка, но она статичная, не векторная).

grid <- st_as_sf(st_make_grid(st_transform(districts_center, crs = 32637), cellsize = 500)) 
grid <- mutate(grid, id=rownames(grid))

grid <- st_intersection(grid, st_union(st_transform(districts_center, crs = 32637)))
streets_grid <- st_intersection(streets, grid)
streets_grid <- select(streets_grid, OSM_ID, NAME, HIGHWAY)

streets_grid <- streets_grid |>
  mutate(len = drop_units(st_length(streets_grid)))
grid <- grid |> 
  mutate(area = drop_units(st_area(grid)))

# Расчет плотности сети
grid_count <- st_join(grid, streets_grid, join = st_contains)
grid_sum <- grid_count |> 
  group_by(id) |> 
  summarise(tot_len = sum(len),
            area = area[1]) |> 
  mutate(dens = tot_len / area)

knitr::kable(head(grid_sum))

library(ggplot2)
ggplot()+
  geom_sf(data=grid_sum, aes(fill=dens))+
  scale_fill_viridis_c(option = "magma",alpha = 0.6,direction = -1)+
  geom_sf(data=districts_center,fill=NA,color="#222277")+
  coord_sf(crs = 32637)+
  theme_void()+
  guides(fill = guide_bins())