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

Сразу скажу, что при помощи гео я решал очень узкий набор задач:

• Оценить существующие локации с т.з. плотности населения, конкуренции, объема рынка. Найти новые точки для открытия или переезда бизнеса;

• Использовать признаки близости покупателя к бизнесу/конкурентам в клиентской аналитике для предсказаний оттока и откликов на рассылки/оффлайн рекламу;

Сегодня расскажу поподробнее про оценку локаций. Все работы я производил на языке R.

Проблема классическая: ищем место для открытия нового магазина (или хотим закрыть старый, переехать). За сторонние сервисы платить не хотим (или почти не хотим). Есть аналитик, он сможет.

Весь процесс можно разбить на несколько простых этапов:

1. Сбор данных по плотности населения;

2. Сбор данных по бизнесу/ конкурентам;

3. Расчет расстояний между бизнесом, конкурентами и населением;

4. Получение метрики оценки качества точки, сравнение и выбор наиболее интересных;

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

Сбор данных по плотности населения

На этом шаге производим сбор домов с координатами и количеством домохозяйств в городе. Изначально я работал с обширной географией. Т.е. просто проанализировать один город и успокоится не получится. Поэтому я настроился анализировать сразу области и республики. Также на обширной географии были интересны населенные пункты от 5 - 10 тыс. человек, где доминирует частный сектор. А частный сектор классически представлен в ГИС хуже многоквартирных домов. В целом нужен был единый подход к поиску детализированной информации.

Приступая к решению задачи, мне повезло натолкнуться на доклад «Пиво против кофе» с указанием отличного источника на кол-во домохозяйств в конкретных многоквартирных домах. Естественно, реестр содержит не все дома, а частный сектор отсутствует напрочь. В итоге, реестр я  объединил с реестром ФИАС. Всем домам из ФИАС без кол-ва домохозяйств я присваивал по одному домохозяйству. Да, школа на районе тоже получит одно домохозяйство, но погрешность будет небольшая. Зато теперь есть информация по частному сектору. У нас страна большая и есть города, где частный сектор – наше всё. Я исключал огороды, дачи и гаражи, т.к. ценности для меня они не представляли. Также я исключал населенные пункты до 300 домохозяйств.  Но в целом можно учитывать и их. Итак, у нас есть набор адресов для геокодирования.

Далее для каждого адреса нужно получить широту и долготу. И тут поджидает первая серьезная проблема. Если пользоваться Open Street Map, то геокодирование в Москве или Питере +- приемлемо, но чем дальше от центра, тем всё печальнее. В итоге, с ходу OSM закрывает около 30-35% домохозяйств, что, согласитесь, маловато. Например, при геокодировании Свердловской области нужно получить данные по 568 тыс. адресов. Это 2 млн. домохозяйств. В 2021 проблема хорошо решалась через API here, который позволял в день выкачивать 250 тыс. адресов. Но в марте 2022 компания ушла из России. Пришлось искать альтернативу. Сейчас из всех геокодеров отлично показывает себя dadata, но кол-во бесплатных запросов в день – 10 тыс., для платных версий – начиная от 50 тыс. Если обрабатывать только адреса с числом домохозяйств от 2-х и более, качество геокодирования поднимается до 85-90% для крупных городов. Если собирать частный сектор, получим еще 5-7%. Финальный штрих – ручное геокодирование адресов от 50 и более домохозяйств. Немного трудозатратно, около 40 минут на регион, но качество геокодирования становится более чем приемлемым, за исключением отдельных городов, где видимо уже ничего не поможет.

Сбор данных по бизнесу / конкурентам

Более простой этап с т.з. подготовки реестров или сбора данных. По аналогии с домами, получаем координаты наших магазинов и магазинов конкурентов.

Есть проблема скорее концептуальная: кого считать конкурентом, а кого нет? Нужно ли учитывать весь рынок или достаточно взять максимально похожих по форматам? Я рассматриваю 2 слоя по отдельности: прямые конкуренты; весь рынок по формату. Более сложные, но (кажется) правильные решения – построение близости бизнеса в рамках моделей Рейли или Хаффа. В них можно учитывать силу бренда, ассортимент, площадь, привлекательность и расстояния до конкурентов. Но пока не приходилось применять их на практике (если у читателей есть опыт, поделитесь пожалуйста).

Расчет расстояний между бизнесом, конкурентами и населением

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

На R я быстро нашел библиотеки osmdata для получения координат дорог, а также sf для объектов Simple Features, в которых изначально происходит сбор географии. Отличным подспорьем были следующие источники тут, тут и тут.

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

Сперва нам нужно получить граф дорог для дальнейшей обработки. Но перед этим опять есть концептуальный вопрос: какую географию мы хотим анализировать: город по его территориальной границе или учитывать еще и близлежащие населенные пункты? И здесь я стараюсь опираться на понимание «торговых зон». Своими словами, торговая зона – радиус вблизи магазина, или бизнеса, где концентрируется основная масса его клиентов. Таких зон м.б. несколько. Очень упрощенно: я использую первичную и вторичную торговую зону как расстояния, в рамках которого генерируются 50% и 80% всех заказов (естественно выкупленных). Более научный и последовательный подход к определениям торговых зон вы можете посмотреть тут (да и в целом прекрасная статья по геопространственному анализу).

Получается, если мы используем концепцию торговых зон, нам нужно анализировать не только сам город, но и близлежащие территории. Поэтому для сбора нужных нам дорог в определенной географической зоне мы будем брать географию города с запасом. Пусть в результате отдельного исследования мы установили, что 80% оборота пиццерия получает от клиентов, проживающих в радиусе 2 км. Тогда к границам города прибавим по 2 км. Границы города я задаю квадратом (по граничным адресам в городе в каждой стороне), передавая 4 координаты в список (bbox) с добавлением 2 км с каждой стороны. Получим квадратную рамку, где будем искать дороги. Далее, bbox передается в запрос OSM и получаем дороги. Есть один нюанс: при геокодинге домов время от времени получаем не верное определение адреса, поэтому я беру не самые крайние значения, а 1% и 99% процентиль по широте и долготе. Таким образом рамка не распухнет в широту или долготу. В то же время, добавление 2 км к границам по перцентилям компенсирует 1% выкинутой с каждой стороны географии.

Чтобы было понятнее, я собрал уже геокодированные адреса в файлы формата parquet, отдельно по 2-м городам в рамках bbox. Предварительно я уже выкинул крайние значения адресов с каждой стороны. Исходные данные можно получить отсюда.

library(data.table)
library(magrittr)
library(dplyr)
library(tidyr)
library(sf)
library(osmdata)
library(tidygraph)
library(sfnetworks)
library(ggplot2)

# Загрузка координат с указанием кол-ва домохозяйств
houses <- arrow::read_parquet("./Каменск-Уральский_homes.parquet") %>%
  data.table()

# определим рамку по координатам для загрузки графа дорог
bbox_city <- st_bbox(
  c(xmin = round(min(as.numeric(houses$lon_home)) - 2000/111000,3),
    xmax = round(max(as.numeric(houses$lon_home)) + 2000/111000,3),
    ymin = round(min(as.numeric(houses$lat_home)) - 2000/111000,3),
    ymax = round(max(as.numeric(houses$lat_home)) + 2000/111000,3)
  ),
  crs = st_crs(4326)
)

# Определим типы дорог для парсинга в OSM 
# https://wiki.openstreetmap.org/wiki/Key:highway
road_type_for_parsing = c(
  "primary",
  "motorway",
  "unclassified",
  "trunk",
  "tertiary",
  "secondary",
  "residential",
  "living_street",
  "motorway_link",
  "trunk_link",
  "primary_link",
  "secondary_link",
  "tertiary_link",
  "footway"
)

# Дороги внутри bbox (в виде объекта sf)
road_osm <- opq(bbox_city) %>%
  add_osm_feature(key = 'highway',
                  value = road_type_for_parsing) %>% 
  osmdata_sf() %>%
  .$osm_lines %>% # забираем только отрезки (lines)
  st_as_sf(, crs = 4326) %>%
  select(osm_id,highway)

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

ggplot(road_osm) +
  geom_sf() +
  ggtitle("Граф дорог г.Каменск-Уральский")

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

Далее я жадно иду через многие методы приведения сети в порядок, заимствуя из разных источников нужные мне приемы.

Отдельно хочу обратить внимание на библиотеку sfnetworks. Она буквально создана для задач: взять граф дорог, причесать его, обогатить данными по домам, магазинам, вести расчеты расстояний на лету между разными объектами. До неё было очень много самописных функций, возни с матчингом графов и т.д. А тут комплексное решение в стиле tidy data.

# отсекаем лишние дороги (вне связной компоненты)
roads_groups <- sf::st_touches(road_osm) %>%
  igraph::graph.adjlist(.) %>%
  igraph::components(.)

roads_table_order <- table(roads_groups$membership) %>%
  .[order(., decreasing = TRUE)]

road_osm_not_comp <- road_osm[
  roads_groups$membership == names(roads_table_order[1]), ]

rm(roads_table_order)


# Чистка графа по источнику:
# https://r-spatial.org/r/2019/09/26/spatial-networks.html

# Функция чистки графа дорог
sf_to_tidygraph = function(x, directed = TRUE) {
  # browser()
  # https://r-spatial.org/r/2019/09/26/spatial-networks.html
  # разбивка графа на отрезки, извлечение всех узлов
  # Give edges a unique index
  edges <- x %>%
    mutate(edgeID = c(1:n()))
  
  # Extract beginning and end node for each edge
  nodes <- edges %>%
    st_coordinates() %>%
    as_tibble() %>%
    rename(edgeID = L1) %>%
    group_by(edgeID) %>%
    slice(c(1, n())) %>%
    ungroup() %>%
    mutate(start_end = rep(c('start', 'end'), times = n()/2)) %>%
    mutate(xy = paste(.$X, .$Y)) %>%
    mutate(xy = factor(xy, levels = unique(xy))) %>%
    group_by(xy) %>%
    mutate(nodeID = cur_group_id()) %>%
    ungroup() %>%
    select(-xy)
  
  source_nodes <- nodes %>%
    filter(start_end == 'start') %>%
    pull(nodeID)
  
  target_nodes <- nodes %>%
    filter(start_end == 'end') %>%
    pull(nodeID)
  
  # Specify for each edge in which node it starts and in which node it ends
  edges <- edges %>%
    mutate(from = source_nodes, to = target_nodes)
  
  # Remove duplicate nodes and convert to sf object with correct CRS
  nodes <- nodes %>%
    distinct(nodeID, .keep_all = TRUE) %>%
    select(-c(edgeID, start_end)) %>%
    st_as_sf(coords = c('X', 'Y')) %>%
    st_set_crs(st_crs(edges))
  
  # Convert to tbl_graph object
  tbl_graph(nodes = nodes, edges = as_tibble(edges), directed = directed)
  
}

graph <- sf_to_tidygraph(road_osm_not_comp, directed = FALSE)

# Преобразование графа через sfnetworks
# https://cran.r-project.org/web/packages/sfnetworks/vignettes/sfn03_join_filter.html
graph <- sfnetworks::as_sfnetwork(graph) %>%
  st_transform(crs = 4326)

graph = convert(graph, to_spatial_subdivision)

# полная чистка от изолянтов
graph <- graph %>%
  activate("edges") %>%
  filter(!edge_is_multiple()) %>%
  activate("nodes") %>%
  filter(!node_is_isolated())

# вырежем все за пределами рамки bbox
graph <-  st_filter(graph , bbox_city  %>% 
                      st_as_sfc() %>%  st_cast("POLYGON")  %>%
                      st_sfc(crs = 4326))

# могут образоваться изолированные маршруты. Поправим
graph <-   graph %>%
  activate("nodes") %>%
  filter(group_components() == 1)

На выходе видим уже более чистую объединенную дорожную сеть.

ggplot() +
  geom_sf(data = st_as_sf(graph, "edges")) +
  ggtitle("Граф дорог г.Каменск-Уральский (очищенный)")

Следующий шаг – от каждого дома мы ищем касательную на граф дорог и добавляем её в граф. Касательная становятся частью графа. Мы фиксируем связку – касательная – дом. Про общую концепцию можно посмотреть тут.

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

# Определим координаты домов
homes_df <- houses %>%
  .[order(lat_home,lon_home)]

# Для домов находим участки сети для блендинга

homes_sf <- homes_df %>%
  unique(.,by = c('lat_home','lon_home')) %>%
  st_as_sf(., coords = c("lon_home","lat_home"),crs = 4326)


blended <- graph %>%
  activate("nodes") %>%
  st_network_blend(., homes_sf)

# получаем id домов и координаты
blended_sf <- blended %>%
  activate("nodes") %>%
  st_as_sf()


homes_graph_point <- data.table(
  home_point_index = blended_sf$home_point_index,
  st_coordinates(blended_sf)) %>%
  .[!is.na(home_point_index)]


# Объединим с основным графом дорог
blended_graph <-  st_network_join(graph,
                                  blended %>% 
                                    activate('nodes') %>%
                                    select(c(geometry))
                                 )

Также в отдельный data.table  сохраняем расстояния между домом и касательной, чтобы при расчетах расстояний учитывать и их (из-за высокой близости к графу я не заморачивался и считал по Евклиду).  Я отсекаю дома, за пределами 1 км от дорожной сети. Они скорее относятся к соседним территориям или строениям, не влияющим на общий характер анализа. Расстояния между домом и графом я сохраняю, чтобы каждый раз не пересчитывать.

# -----(Код объединения расчета расстояний).---------------
# быстрая отсечка домов, далеких от графа дорог
# считаем расстояние по евклиду от дома до ближайшего ребра графа дорог
cjdt <- function(a,b){
  cj = CJ(1:nrow(a),1:nrow(b))
  cbind(a[cj[[1]],],b[cj[[2]],])
}

# Считаем попарно расстояния от домов до точек на графе
blended_point <- blended_graph %>% 
  activate(nodes) %>%
  st_coordinates() %>%
  data.table()

# берем 1 градус как 111 км (множитель 111000)
dist_matrix_home_to_road <- cjdt(blended_point[,.(lat_road = Y, lon_road = X)],
                    homes_df[,.(lat_home,lon_home,home_point_index)]) %>%
  .[, dist := sqrt( (lat_home - lat_road)^2 + (lon_home - lon_road)^2) * 111000] %>%
  .[dist<1000] %>% # Отсечка далеких от графа домов
  .[order(dist)] %>% 
  unique(., by = c('home_point_index'))

Аналогично, в дальнейшем посчитаем расстояния от дорог до пиццерий. Т.е. расстояние, которое будет между бизнесом и домохозяйствами = расстояние от бизнеса до дорожного графа + расстояние по дорожному графу + расстояние от графа до дома.

И еще маленькая эвристика: мы же ищем, где открыться. Как производить поиск? Есть подход от тех, кто в нашей сети пиццерий занимается арендными отношениями: «Вот вам 3 адреса, где мы нашли помещение. Оцените!». Попробуем все же оценить всю картину целиком, чтобы точечные оценки нас не пугали: на карту нанесем плотную сетку из точек по всему bbox (напр. через каждые 200 метров). Тогда мы сможем использовать накинутые точки как ориентиры для поиска, при расчете сценария открытия. Проделываем ту же операцию, что и с домами, только по ограниченному числу точек за раз, чтобы не выбить оперативную память (напр. по 5000).

# Рассчитаем координаты точек (по аналогии с домами)
potencial_points <- expand.grid(
  lon_point = seq(from = bbox_city[1], to = bbox_city[3], by = 200/111000),
  lat_point = seq(from = bbox_city[2], to = bbox_city[4], by = 200/111000)) %>%
  data.table() %>%
  .[,point_ind := seq(.N)]%>%
  .[order(lat_point,lon_point)]%>%
  .[,iteration := seq(.N) %/% 5000]
# итерации нужны в дальнейшем, чтобы не перегружать оперативную память для крупных городов


# Для потенциальных точек по итерациям находим участки сети для блендинга
graph_point <- data.table() #для хранения связи индекс потенциальная точка - точка на графе

for(iter in c(0:max(potencial_points$iteration))){
  potencial_points_sf <- potencial_points[iteration == iter] %>%
    unique(.,by = c('lat_point','lon_point')) %>%
    st_as_sf(., coords = c("lon_point","lat_point"),crs = 4326)
  
  min_lat <- min(potencial_points[iteration == 0]$lat_point) - 2000/111000
  max_lat <- max(potencial_points[iteration == 0]$lat_point) + 2000/111000
  min_lon <- min(potencial_points[iteration == 0]$lon_point) - 2000/111000
  max_lon <- max(potencial_points[iteration == 0]$lon_point) + 2000/111000
  
  blended <- graph %>%
    activate("nodes") %>%
    filter(node_Y() >= min_lat & node_Y() <= max_lat &
           node_X() >= min_lon & node_X() <= max_lon) %>%
    st_network_blend(., potencial_points_sf)
  
  # получаем id потенциальных точек и координаты
  blended_sf <- blended %>%
    activate("nodes") %>%
    st_as_sf()
  
  
  blended_df <- data.table(
    point_ind = blended_sf$point_ind,
    st_coordinates(blended_sf)
  ) %>%
    .[!is.na(point_ind)]
  
  graph_point <- rbindlist(list(
    graph_point,blended_df
  ))
  
  # Объединим с основным графом
  blended_graph <-  st_network_join(blended_graph,
                                    blended %>% activate('nodes') %>%
                                      select(c(geometry))
  )
}

# быстрая отсечка потенциальных точек, далеких от графа дорог
# Считаем попарно расстояния от домов до точек на графе
pre_dist_matrix_point <- data.table()
potencial_points[,iteration := seq(.N) %/% 100]

blended_point <- blended_graph %>% 
  activate(nodes) %>%
  st_coordinates() %>%
  data.table()

dist_matrix_point_to_road <- data.table()

for(i in c(0:max(potencial_points$iteration))){
  time_matrix <- cjdt(blended_point[,.(lat_road = Y, lon_road = X)],
                      potencial_points[iteration == i ,.(lat_point,lon_point,point_ind)]) %>%
    .[, dist := sqrt( (lat_point - lat_road)^2 + (lon_point - lon_road)^2)*111000 ]%>%
    .[dist<1000] %>% 
    .[order(dist)] %>% 
    unique(., by = c('point_ind'))
  
  dist_matrix_point_to_road <-  rbindlist(
    list(
      dist_matrix_point_to_road,
      time_matrix
    )
  )
}

# Применим морфер для отображения не конечных точек. Подробнее:
# https://cran.r-project.org/web/packages/sfnetworks/vignettes/sfn05_morphers.html
blended_graph = convert(blended_graph, to_spatial_subdivision)

Итак, на этом этапе у нас есть объединенные в один граф: дороги, дома, точки для поиска локации. Это некоторые стандартные вещи, которые мы можем не пересчитывать для каждого бизнеса.

Теперь для каждой потенциальной точки найдем список домов, которые лежат недалеко от неё (в пределах торговой зоны). Пример сильно не оптимизированного кода:

# Торговая зона (максимально допустимое расстояние до бизнеса)
max_distance_to_point <- 2000

# Добавим кол-во домохозяйств в dt по расстоянию домов от графа дорог
dist_matrix_home_to_road[homes_df[,.(quarters_count = sum(quarters_count)),
                                 by = c('home_point_index')],
             `:=`(quarters_count = i.quarters_count),
             on = c("home_point_index")]

# Глобально посчитаем по каждому узлу графа расстояние 
# до каждого узла с касательной потенциальной точки
list_coord_road <- st_as_sf(x =  dist_matrix_point_to_road,                         
                              coords = c("lon_road", "lat_road"),crs = 4326)


# в граф потенциальной точкой открытия
dist_table <-  st_network_cost(blended_graph,
                               from = list_coord_road,
                               direction = "all") %>%
  t() %>%
  data.table() %>%
  .[,point_ind := seq(.N)] %>%
  melt(.,id.vars = c('point_ind') ) %>%
  .[,distance_road := as.numeric(value)] %>%
  .[distance_road <= max_distance_to_point] %>%
  .[,variable := as.integer(gsub('V','',variable))] 

# Добавим координаты касательной к дорожным расстояниям
road_coord <- blended_graph %>%
  activate('nodes') %>%
  st_geometry() %>%
  st_coordinates() %>%
  data.table() %>%
  .[,road_point_ind := seq(.N)]

dist_table[road_coord,
           `:=`(road_point_lat = i.Y,
                road_point_lon = i.X),
           on = .(point_ind = road_point_ind)]

dist_table[,c('value','point_ind') := NULL]

# для каждого дома (по координатам) подготовим набор потенциальных точек
homes_point_dist <- merge.data.table(x = dist_matrix_home_to_road,
                                     y = dist_table,
                                     by.x = c('lat_road','lon_road'),
                                     by.y = c('road_point_lat','road_point_lon'),
                                     all.x = T,
                                     all.y = F,
                                     allow.cartesian = T) %>%
  .[!is.na(distance_road)] %>%
  .[(distance_road + dist) <= max_distance_to_point] %>%
  .[,.(home_road_point_lat = lat_road,
       home_road_point_lon = lon_road,
       home_lat = lat_home,
       home_lon = lon_home,
       home_point_index,
       distance_home_to_road = round(dist,0),
       quarters_count,
       mag_ind = variable,
       distance_road = round(distance_road,0)
  )]

# Добавим конкретику по потенциальной точке 
# (для удобства дальнейшего представления данных)
homes_point_dist[
  dist_matrix_point_to_road[,
                .(business_type = 'new_point',
                  mag_name = '',
                  mag_type = 'new_point',
                  mag_ind = seq(.N),
                  mag_road_point_lat = lat_road,
                  mag_road_point_lon = lon_road,
                  mag_lon = lon_point,
                  mag_lat = lat_point,
                  mag_dist_to_road = round(dist,0),
                  mag_address = '')],
  `:=`(business_type = i.business_type,
       mag_name = i.mag_name,
       mag_type = i.mag_type,
       mag_lon = i.mag_lon,
       mag_lat = i.mag_lat,
       mag_road_point_lat = i.mag_road_point_lat,
       mag_road_point_lon = i.mag_road_point_lon,
       mag_dist_to_road = i.mag_dist_to_road,
       mag_address = i.mag_address
  ),
  on = .(mag_ind = mag_ind)
]

# отфильтруем все что в пределах торговой зоны с учетом всех типов расстояний
# от дома до графа, расстояния по графу, от графа до нового магазина
homes_point_dist = homes_point_dist[(distance_road + 
                                    distance_home_to_road +
                                    mag_dist_to_road) <= max_distance_to_point]

На старте у нас 6 тыс. домов, 5 тыс. точек для оценки. В итоге отфильтрованных пар потенциальных точек – координат домов – всего 1,4 млн. Все остальное – вне торговых зон.

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

library(tmap)

# получим только потенциальные локации с кол-вом домохозяйств в торговой зоне от 500
point_quarters_density <-  homes_point_dist %>%
  .[,.(quarters = sum(quarters_count)),by = c("mag_lon","mag_lat")] %>%
  .[quarters > 500] %>%
  st_as_sf(., coords = c("mag_lon","mag_lat"),crs = 4326)

tmap_mode('view')
tm_shape(point_quarters_density) +
  tm_dots(col = "quarters",size = 0.06) +
  tmap_options(basemaps = 'OpenStreetMap')

Четко видны 2 основных района, где в радиусе 2 км ожидается высокая плотность населения. Пока картина не особо полезна. Видим два крупных района застройки. Да, где-то в радиусе будет более 30 тыс. жителей. Но и конкуренты распределены не равномерно. Пора их добавить.

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

# Добавление существующего бизнеса
magazine_dt <- fread('./Каменск-Уральский_пиццерии.csv',
                     integer64 = 'character') %>%
  .[,.(name_business = `Наименование`,
       lon = `Lon`,
       lat = `Lat`)]

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

# Проиндексируем магазины
magazine_dt[,magazine_ind := seq(.N)]

# найдем ближайший дом, используем как точку вхождения
#  считаем расстояние по евклиду 
dist_matrix <- cjdt(magazine_dt[,.(magazine_ind,lat_b = lat, lon_b = lon)],
                    blended_graph %>%
                      activate(nodes) %>%
                      st_coordinates() %>%
                      data.table() %>%
                      .[,.(lat_road = Y,lon_road = X)]) %>%
  .[, dist := sqrt( (lat_road - lat_b)^2 + (lon_road - lon_b)^2)*111000]%>%
  .[order(dist)] %>% 
  unique(., by = c('magazine_ind'))

# Добавим расстояния до графа
magazine_dt[dist_matrix,
            `:=`(lat_road = i.lat_road,
                 lon_road = i.lon_road,
                 dist_to_road = i.dist),
            on = .(lat = lat_b, lon = lon_b)]

И по аналогии с сеткой из точек ищем для каждого магазина набор домов в торговой зоне.

# Для каждого бизнеса посчитаем дорожное расстояние до точек домов
stores_list_coord <- st_as_sf(x = magazine_dt,                         
                              coords = c("lon_road", "lat_road"),crs = 4326)

# расчет расстояния магазинов до домов
dist_table <-  st_network_cost(blended_graph, 
                               from = stores_list_coord,direction = "all") %>%
  t() %>%
  data.table() %>%
  .[,point_ind := seq(.N)] %>%
  melt(.,id.vars = c('point_ind') ) %>%
  .[,distance_road := as.numeric(value)] %>%
  .[distance_road <= max_distance_to_point] %>%
  .[,variable := as.integer(gsub('V','',variable))]

# Добавим координаты точек вхождения
road_coord <- blended_graph %>%
  activate('nodes') %>%
  st_geometry() %>%
  st_coordinates() %>%
  data.table() %>%
  .[,road_point_ind := seq(.N)]

dist_table[road_coord,
           `:=`(road_point_lat = i.Y,
                road_point_lon = i.X),
           on = .(point_ind = road_point_ind)]

dist_table[,c('value','point_ind') := NULL]

# для каждого дома (по координатам) найдем все магазины в допустимой торговой зоне
homes_mag_dist <- merge.data.table(x = dist_matrix_home_to_road,
                                   y = dist_table,
                                   by.x = c('lat_road','lon_road'),
                                   by.y = c('road_point_lat','road_point_lon'),
                                   all.x = T,
                                   all.y = F,
                                   allow.cartesian = T) %>%
  .[!is.na(distance_road)] %>%
  .[(distance_road + dist) <= max_distance_to_point] %>%
  .[,.(home_road_point_lat = lat_road,
       home_road_point_lon = lon_road,
       home_lat = lat_home,
       home_lon = lon_home,
       distance_home_to_road = dist,
       quarters_count,
       mag_ind = variable,
       distance_road
  )]

# Добавим информацию по бизнесу
homes_mag_dist[
  magazine_dt[,.(
    mag_name = name_business,
    mag_ind = seq(.N),
    mag_road_point_lat = lat_road,
    mag_road_point_lon = lon_road,
    mag_lon = lon,
    mag_lat = lat,
    mag_dist_to_road = dist_to_road)],
  `:=`(
    mag_name = i.mag_name,
    mag_lon = i.mag_lon,
    mag_lat = i.mag_lat,
    mag_road_point_lat = i.mag_road_point_lat,
    mag_road_point_lon = i.mag_road_point_lon,
    mag_dist_to_road = i.mag_dist_to_road
  ),
  on = .(mag_ind = mag_ind)
]

# отфильтруем все что в пределах торговой зоны с учетом всех типов расстояний
# от дома до графа, расстояния по графу, от графа до магазина
homes_mag_dist <- homes_mag_dist[(distance_road + 
                                  distance_home_to_road +
                                  mag_dist_to_road) <= max_distance_to_point]

Получение метрики оценки качества точки, сравнение и выбор наиболее интересных локаций

Теперь у нас есть по каждому домохозяйству набор пиццерий в торговой зоне. Будем считать конкуренцию следующим образом: по каждому домохозяйству находим среднее число конкурентов в торговой зоне от домохозяйства (т.е. клиент будет выбирать то, что близко к нему). Для каждой пиццерии считаем средневзвешенное число пиццерий по домохозяйствам в торговой зоне от расчетной пиццерии. Чтобы было понятнее, простой пример логики расчета. Пусть есть 2 пиццерии [piz1, piz2] и 3 домохозяйства в радиусе от которых есть [ [piz1], [piz1, piz2], [piz1, piz2] ]. Тогда для домохозяйств число конкурентов: [ 1, 2, 2]. Для пиццерий показатель конкуренции будет [ (1+2+2)/3, (2 + 2)/2 ] = [1.66, 2]. Т.е. пиццерия 2 в более для нее конкурентной среде. А чтобы понять, сколько приходится в среднем на 1 бизнес домохозяйств, просто поделим охват на средневзвешенную конкуренцию [1.81, 1]

# добавим кол-во магазинов на каждое домохозяйство
homes_mag_dist[,n_business := .N, 
               by = c('home_lat','home_lon')]

# для поддержания уникальности магазина создадим уникальное имя
homes_mag_dist[,mag_uniq := paste0(mag_name, ' ', mag_lon, ' ',mag_lat)]

# Получим агрегаты по магазинам
business_aggr <- homes_mag_dist[,.(quarters_count = sum(quarters_count),
                                   avg_magazine_in_trading_area = mean(n_business)),
                                by = c('mag_name','mag_lat','mag_lon')] %>%
  .[, avg_quarters_on_magazine := round(quarters_count 
                                        / avg_magazine_in_trading_area,2)]%>%
  .[order(-quarters_count)] %>%
  .[,.(`Пиццерия` = mag_name,
       `Широта` = mag_lat,
       `Долгота` = mag_lon,
       `Кол-во домох-в в торг.зоне` = quarters_count,
       `Средневзвешенное по домох-вам число конкурентов в торг. зоне` = round(avg_magazine_in_trading_area,2),
       `Охват с учетом конкуренции` = avg_quarters_on_magazine
       )]

Посмотрим уровень конкуренции наших данных:

Интересно, что Додо не боятся конкуренции) Также видно что одна из пиццерий единолично занимает территорию. Остальные серьёзно делят рынок в другой части города.

mag_quarters_density <-  homes_mag_dist %>%
  .[,.(quarters = sum(quarters_count)),by = c("mag_lon","mag_lat")] %>%
  .[quarters > 500] %>%
  st_as_sf(., coords = c("mag_lon","mag_lat"),crs = 4326)

buffers <- st_buffer(mag_quarters_density, 2000)

tmap_mode('view')
tm_shape(buffers)+
  tm_polygons(alpha = 0.5,
              col = "green")+
  tm_shape(mag_quarters_density) +
  tm_dots(col = "quarters",size = 0.06) +
  tmap_options(basemaps = 'OpenStreetMap')

Кто-то скажет, что тут и так все очевидно, зайди на карту да посмотри. Ок, в данном случае да. Но если пиццерий в городе 20, а застройка (как тут) не равномерна и есть интересные нишы? И мы в целом можем сравнивать нишы в разных городах, перейдя на объем рынка (об этом чуть позже) и выигрывать за счет того, что видим картину более масштабно.

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

# Установим порог на не интересные сценарии (в домохозяйствах)
min_people_count <- 1000
# Добавление сценариев открытия
new_points <- homes_point_dist  %>% 
  .[,mag_uniq := 'new_point'] %>% 
  .[,mag_name := 'new_point']

new_points[,quarters_count_point := sum(quarters_count), by = c('mag_ind')]
new_points <- new_points[quarters_count_point > min_people_count]
new_points[,quarters_count_point:= NULL]

new_points_list <- unique(new_points$mag_ind)

# В параллели рассчитываем сценарии
library(foreach)
library(doParallel)
registerDoParallel(cl <- makeCluster(3))
new_point_calculate <- foreach(
  i = new_points_list,
  .packages = c("data.table", "magrittr"),
  .combine=function(x,y)rbindlist(list(x,y)),
  .inorder = F,
  .export = c('new_points','homes_mag_dist')
) %dopar% {
  
  rbindlist(list(
    new_points[mag_ind == i] ,
    homes_mag_dist
  ),use.names = T, fill = T) %>%
    .[,n_business := .N, 
      by = c('home_lat','home_lon')] %>%
    .[,.(quarters_count = sum(quarters_count),
         avg_magazine_in_trading_area = round(mean(n_business),3)),
      by = c('mag_name','mag_uniq','mag_lat','mag_lon')] %>%
    .[, avg_quarters_on_magazine := round(quarters_count 
                                          / avg_magazine_in_trading_area,
                                          2)] %>%
    .[order(-quarters_count)] %>%
    .[, point_calculate_lon := new_points[mag_ind == i]$mag_lon[1]]  %>%
    .[, point_calculate_lat := new_points[mag_ind == i]$mag_lat[1]]
  
}
stopCluster(cl)

# Посмотрим только на сценарии от определенного порога
new_points_full <- new_point_calculate[mag_name == 'new_point' &
                                       avg_quarters_on_magazine >= 10000]
new_points_full_sf <- st_as_sf(new_points_full, 
                               coords = c("mag_lon","mag_lat"),crs = 4326)  

Отобразим на карте:

tmap_mode('view')
tm_shape(new_points_full_sf) +
  tm_dots(col = "avg_quarters_on_magazine",size = 0.02) +
  tmap_options(basemaps = 'OpenStreetMap')

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

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

Итак, вернемся к изначальному вымышленному примеру. У нас есть 2 города, где мы можем разместиться. Я не буду дублировать код по Новоуральску. Приложу только исходные данные по домохозяйствам и координатам. В Новоуральске ситуация с пиццериями хуже, я нашел всего 3 релевантных. После анализа ситуация следующая. Топ 100 сценариев по Каменск-Уральский: лучший сценарий – 14695 домохозяйства с учетом конкурентов, худший сценарий -  12002. По Новоуральску: 8924 и 6520 соответственно. С одной стороны, можно было бы и закончить на этом. Но есть еще разница в доходах и доля рабочего населения. Я копнул gorodrabot и выяснил, что на 24.03.2024 средняя зп в Новоуральске = 68 058 руб., а в Каменск-Уральский = 59 823 руб. Если заложить разницу в ЗП, это не изменит порядок сценариев. Хотя на моей практике часто сценарии из глубинки пробивали сценарии из столиц областей и республик.

Ну что же, рекомендуем Каменск-Уральский и конкретные торговые зоны + можем подсветить, с кем конкурировать. Я бы постеснялся вставать рядом с федеральной сеткой с сильным брендом, доставкой и бизнес процессами. И это тоже можно отобразить в модели, придав им больший вес.

В заключение

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

• гео требует много ресурсов, поэтому крупному бизнесу проще купить готовые решения, выйдет дешевле. Среднему бизнесу сложнее. Поэтому гео часто заменяется более простой методологией либо реализуется на своей стороне;

• часто субъективная оценка заказчика мешает корректно оценивать нишы. В моей практике был случай, когда на место, которое рекомендовал, отвечали резким отказом именно из-за субъективной оценки городской среды. Потом это место заняли федералы. Но этот аспект относится скорее к аналитике в целом;

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

Я замечал, что гео пространственной аналитикой занимаются многие федеральные сети, четко понимая, какую долю рынка на конкретной территории они будут занимать. И для меня вывод очевиден – этим стоит заниматься, это отличное направление аналитики для всего, что связано с оффлайн, розницей, общепитом, доставкой. Это дорогой, но полезный кирпичик в построении стратегии принятия решений. А бонусом – масса спиноффов в виде приложений к CRM, понимания конкуренции, анализа ряда оффлайн каналов привлечения трафика.