Введение

Во французском сабреддите я наткнулся на пост с интересной гипотезой:

Чем ближе точка к вокзалу, тем хуже там кебаб.

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

Я решил, что мне нечем заняться, ведь я вылетевший с учёбы выгоревший безработный с новоприобретённым диагнозом «аутизм», поэтому стоит пожертвовать немного своего времени на благую цель — проведение этого неформального «исследования». В пределах следующих трёх рабочих дней мне на почту точно придёт Нобелевская премия мира и куча рабочих оферов.

Методика

Я решил, что лучшим местом для исследования будет Париж. Причины:

1. Исходный пост был французским

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

2. Плотность

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

Данные сети

При помощи OSMnx я скачал и сохранил сеть навигации. Учитывая то, что французский сабреддит сильно центрирован на общественном транспорте, я решил, что будет логично будет учитывать пешеходные расстояния (то есть пешеходные дорожки и тротуары), поэтому задал для network_type OSMnx значение «walk». Учитывая местоположение (и то, что OSMnx автоматически использовал эту CRS (базовую систему координат), если не было указано иное), все данные были спроецированы согласно EPSG:32631 (зона UTM 31N).

import osmnx as ox
from geopandas import GeoDataFrame

#EPSG
PROJECTION = 32631

graph = ox.graph_from_place('Paris, FR', network_type="walk")
graph = ox.project_graph(graph, to_crs=PROJECTION)

ox.save_graphml(graph, filepath="network.graphml")

Дальше нужно разобраться с различными вокзалами/станциями метро. Учитывая особенности сабреддита, я решил, что логично будет использовать и центральные вокзалы дальнего следования, и бесчисленные станции метро. Эта задача тоже достаточно тривиально решалась при помощи OSMnx, нужно было всего лишь выполнить фильтрацию по railway=subway_entrance или railway=train_station_entrance.

stations: GeoDataFrame = ox.features_from_place('Paris, FR', tags = {
    "railway": ["subway_entrance", "train_station_entrance"]
})

# Фильтруем результаты по точкам
station_nodes: GeoDataFrame = stations.loc[stations.geom_type=="Point"]
station_nodes = station_nodes.to_crs(epsg=PROJECTION)

station_nodes.to_file("train_station_entrances.gpkg")

Я тщательно сохранил все выходные данные, чтобы можно было легко изучить их в QGIS. Также я попробовал заставить Python-ноутбуки работать с моей конфигурацией NeoVIM, но все усилия пропали даром.

... И теперь у нас есть первая половина данных. Далее перейдём к заведениям общепита.

Данные заведений

Мне показалось, что разумно будет использовать Google Places API (и его данные отзывов). Отзывы в Google, естественно, далеки от идеальных, большая их написана ботами, не говоря о других проблемах, но это лучшее, что я смог на тот момент придумать. Есть альтернативы, например Yelp, но их API ужасно дороги для нищего меня, а веб-скрейпер у меня настроения писать не было (этот процесс так же угнетает меня, как промтинг LLM). Привлёк меня и бесплатный кредит Google в 200 долларов.

Однако когда я начал исследовать API, то осознал, что в Places API, похоже, нет никаких способов поиска внутри многоугольника, только внутри радиуса точки. Спасибо тебе, мегакорпорация, очень удобно.

Не способствовало дальнейшей разработке и неработающее автозаполнение для библиотеки googlemaps. Python — хороший язык, но его инструментарий слишком часто любит проверять моё терпение. И если уж мы начали жаловаться, скажу, что дэшборд Google Cloud — это, наверно, самый медленный «веб-сайт», с которым я имел неудовольствие взаимодействовать.

Итак... Это значит, что мне придётся выполнять некий поиск по сетке по всему Парижу, надеясь, что это не исчерпает мой кредит. К тому же, возникла пара интересных вопросов:

1. Что такое «кебаб»?

Когда я выполняю поиска «kebab» (дополнительный контекст не требуется)... то как Google определяет, в каких заведениях подают кебаб?

Изучив вопрос, я понял, что всё не так сложно, как мне казалось. В названиях многих заведений просто есть слово «kebab», некоторые помечены как «Mediterranean» («средиземноморское»; родина кебаба — Турция, Персия и Ближний Восток в целом), а у других было довольно много отзывов со словом «kebab». Меня это вполне устраивало.

2. Проблема в мире запросов

Оказалось, что при запросе мест в пределах указанного радиуса это только «смещение», а не чёткая граница, помогающая сузить диапазон сбора данных и снизить количество ненужных запросов. Мне становилось всё очевиднее, что Google не любит, когда люди занимаются что-то подобным.

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

Административные границы Парижа включают в себя пару крупных зелёных областей. На западе это парк, а на востоке — какой-то спортивный институт.

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

Я решил, что нет ничего страшного в том, чтобы оставить сеть и станции, поэтому ничего не менял.

Чтобы максимизировать объём собираемых данных, я выбрал шестиугольную сетку с расстоянием по вертикали 1 км. Это обеспечит нам радиус поиска 500 м * √3 ~= 866 метров. Да, радиусы будут достаточно сильно пересекаться, зато у нас не будет никаких дырок. Не знаю, зачем я тратил так много времени на обеспечение «целостности данных», ведь из-за Google процесс всё равно мог пойти наперекосяк, но мне важна была иллюзия контроля.

Так мы получили 99 точек сэмплирования. Наверно, этого должно хватить?

Впрочем, ладно. Вот, как выглядел мой код на Python, написанный в три часа ночи:

# Уже спроецировано на WGS84
sample_points: GeoDataFrame = GeoDataFrame.from_file("samples.gpkg")
gmaps: googlemaps.Client = googlemaps.Client(key='get-your-own')

output = {}

for point in sample_points.geometry:
    lat, lon = point.y, point.x

    next_page_token = None
    num_fetches = 3

    while num_fetches > 0:
        result = {}

        if next_page_token == None:
            result = gmaps.places(
                "kebab",
                location=(lat, lon),
                radius=866,
            )
        else:
            result = gmaps.places(
                page_token=next_page_token
            )

        next_page_token = result.get("next_page_token")
        print(result["status"], next_page_token)

        for p in result["results"]:
            output[p["place_id"]] = p

        if next_page_token == None:
            break

        num_fetches -= 1

        sleep(2)

json_out = json.dumps(output)

with open("output.json", "w") as file:
    file.write(json_out)

Код сработал довольно неплохо. Изначально я не стал заниматься разбиением на страницы, получив 322 результата. Однако я заметил, что некоторые заведения, найденные мной на Google Maps, в списке результатов отсутствовали.

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

Сырой вывод Google Places API нужно было также усечь по области исследования, спроецированной на локальную зону UTM, а также преобразованной в геопространственный формат:

import pandas as pd

with open("output.json", "r") as file:
    data = json.load(file)
    file.close()
    
    for id in data:
        place = data[id]
        point = place["geometry"]["location"]
        data[id]["lng"] = point["lng"]
        data[id]["lat"] = point["lat"]
        del data[id]["geometry"]

    data = pd.DataFrame.from_dict(data).T
    data.rating = pd.to_numeric(data.rating)
    data.user_ratings_total = pd.to_numeric(data.user_ratings_total)
    data = data[data["user_ratings_total"] > 0]

    # Очистка была добавлена после того, как сделан показанный ниже скриншот
    data = data.drop(columns=[
        "icon",
        "icon_background_color",
        "icon_mask_base_uri",
        "plus_code",
        "reference",
        "photos",
        "opening_hours"
    ])

    gdata = GeoDataFrame(
        data, geometry=geopandas.points_from_xy(data.lng, data.lat),
        crs=4326
    )

    gdata: GeoDataFrame = gdata.to_crs(PROJECTION)

    # Изменённые границы с рисунка 4.
    paris = GeoDataFrame.from_file("mod_bounary.gpkg");

    gdata: GeoDataFrame = gdata.clip(paris)

    gdata.to_file("establishments.gpkg")

Маршруты и расстояние

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

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

Погуглив и покопавшись в документации API я, однако, понял, что для этой цели нам вполне подойдёт sindex библиотеки GeoPandas. Хоть в ней и нет возможности «верни ближайшие N», как в моей любимой Rust-библиотеке r-tree, к которой я так привык, она позволила мне выполнять поиск в пределах определённого радиуса (одного километра мне было вполне достаточно), и после этого двигаться дальше. Результаты запросов не отсортировывались, поэтому мне приходилось сортировать индексы по расстоянию и усекать их по нужному размеру.

Благодаря NetworkX анализ сетей оказался достаточно простым процессом, и спустя пару часов я смог составить такой код:

import networkx as nx
import shapely as shp

establishments: GeoDataFrame = GeoDataFrame.from_file("establishments.gpkg")
entrances: GeoDataFrame = GeoDataFrame.from_file("entrances.gpkg")
graph = ox.load_graphml("network.graphml")
graph = ox.project_graph(graph, to_crs = PROJECTION)

# Гарантируем, что используется та же CRS
if (establishments.crs != entrances.crs != PROJECTION):
    exit(100)

# Вспомогательная функция для получения расстояния между узлом графа и геометрией заведения
def node_geom_dist(node_id: int, geom: shp.Point):
    node = graph.nodes[node_id]
    return math.sqrt((geom.x - node['x']) ** 2 + (geom.y - node['y']) ** 2)

distances: list[float] = []

for (id, establishment) in establishments.iterrows():
    establishment_geom: shp.Point = establishment.geometry
    establishment_node: int = ox.nearest_nodes(graph, establishment_geom.x, establishment_geom.y)
    establishment_dist_to_node: float = node_geom_dist(establishment_node, establishment_geom)
    
    # Пространственный индекс для входов в станции
    index: shp.STRtree = entrances.sindex
    nearest_q = index.query(establishment_geom, predicate="dwithin", distance = 1000)
    nearest_entrances: list[tuple[int, float]] = []

    for i in nearest_q:
        ent = entrances.iloc[i]
        ent_geom: shp.Point = ent.geometry

        dist = ent_geom.distance(establishment.geometry)
        
        nearest_entrances.append((i, dist))
     
    nearest_entrances = sorted(nearest_entrances, key = lambda e: e[1])[:3]
    entrance_geom: list[shp.Point] = [entrances.iloc[i].geometry for (i, _) in nearest_entrances]
    entrance_nodes: list[int] = [ox.nearest_nodes(graph, point.x, point.y) for point in entrance_geom]
    entrance_geom_dist_to_node: list[float] = [node_geom_dist(entrance_nodes[i], entrance_geom[i]) for i in range(len(nearest_entrances))]

    result_paths = [nx.shortest_path(graph, establishment_node, dest_node, weight="length") for dest_node in entrance_nodes]
    result_lengths: list[float] = [nx.path_weight(graph, path, "length") + entrance_geom_dist_to_node[i] + establishment_dist_to_node for (i, path) in enumerate(result_paths)]

    distances.append(min(result_lengths))

establishments["distance"] = distances 
establishments.to_file("establishment_results.gpkg")

Не самый лучший мой код. Общий объём возни со списками, наверно, немного пугает, зато всё работает.

Потестировав получившиеся данные и сети в QGIS (и использовав множество операторов print()), я обрёл уверенность в точности результатов.

Результаты

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

Благодаря могучему коэффициенту корреляции Пирсона в 0,091 можно считать, что это может быть правдой! Если только закрыть глаза на то, что корреляция настолько слаба, что назвать её «статистически незначимой» было бы слишком щедрым заявлением.

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

... что увеличило коэффициент аж до целых 0,098.

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

Коэффициент ухудшился, став ниже на 0,001, так что, думаю, настало время выбросить белый флаг.

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

Вопросы и ответы

- Действительно ли отзывы Google — объективная метрика вкуса кебабов?

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

- Может ли оказывать влияние туризм?

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

- Были ли результаты Google точными?

Да, в определённой мере. Судя по моим наблюдениям, все точки из запроса действительно в том или ином виде подавали кебаб. Было несколько странных выбросов и нюансов, например Pizza Hut, который, вероятно, подаёт кебаб-пиццу, а не один из множества различных видов кебабов.

- Почему не заведения общепита в целом?

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

- А как насчёт количества отзывов?

Разумеется, это вполне могло оказать влияние, но на момент проведения исследований я не совсем понимал, как внедрить эту метрику в анализ.

- Дай нам данные

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

К тому же я показал весь код. Сжигайте сами свои кредиты.

- У тебя там всё нормально вообще?

... Да вроде. А у тебя?

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

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

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

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

Ждите второй части... Когда-нибудь я её напишу!