Spatial join: соединение по положению

Урок про spatial join — операцию, которая соединяет два слоя не по общему ключу, а по тому, как их геометрии расположены друг относительно друга.

Spatial join (пространственное соединение) присоединяет атрибуты объектов одного слоя к объектам другого на основании их пространственного отношения (содержит, пересекает, ближайший).

Если вы знаете SQL JOIN, идея знакома: соединить две таблицы. Но в обычном JOIN ключ — совпадающее значение (id = id). В spatial join «ключ» — это пространственное отношение: точка попадает в полигон, линия пересекает зону. Классический пример: есть тысячи точек-инцидентов (ДТП, заявки) и слой районов; spatial join проставит каждой точке, в каком она районе, — после чего можно группировать и считать по районам.

Логика операции

Для каждого объекта левого слоя ищем объект правого слоя, удовлетворяющий предикату (чаще всего within — «точка внутри полигона»), и приклеиваем его атрибуты. Реализуем на основе нашего point_in_polygon:

def point_in_polygon(x, y, poly):
    n = len(poly)
    inside = False
    j = n - 1
    for i in range(n):
        xi, yi = poly[i]
        xj, yj = poly[j]
        if (yi > y) != (yj > y):
            if x < (xj - xi) * (y - yi) / (yj - yi) + xi:
                inside = not inside
        j = i
    return inside

# Слой районов (полигоны)
districts = {
    "Север": [(0, 5), (10, 5), (10, 10), (0, 10)],
    "Юг":    [(0, 0), (10, 0), (10, 5), (0, 5)],
}
# Слой инцидентов (точки)
incidents = [("заявка-1", 2, 7), ("заявка-2", 8, 2),
             ("заявка-3", 5, 9), ("заявка-4", 6, 1)]

from collections import Counter
counts = Counter()
for name, x, y in incidents:
    district = next((d for d, poly in districts.items()
                     if point_in_polygon(x, y, poly)), "вне зоны")
    print(f"{name} -> {district}")
    counts[district] += 1

print("--- по районам ---")
for d, c in counts.items():
    print(f"{d}: {c}")

Вывод:

заявка-1 -> Север
заявка-2 -> Юг
заявка-3 -> Север
заявка-4 -> Юг
--- по районам ---
Север: 2
Юг: 2

Мы только что выполнили полный конвейер: spatial join (точка → район) плюс агрегация (счётчик по районам). В промышленном коде это две строки geopandas, но логика ровно та же.

Виды spatial join

Как работает под капотом

Наивный spatial join — это вложенный цикл: для каждой из $m$ точек проверить каждый из $n$ полигонов, итого $O(m \cdot n)$ проверок «точка в полигоне», каждая ещё и $O(k)$ по вершинам. На реальных данных это неприемлемо, поэтому правый слой индексируют R-деревом: для точки сначала по bounding box находят полигоны-кандидаты (обычно один-два), и только их проверяют точно. Так $O(m \cdot n)$ сжимается почти до $O(m \log n)$. Именно поэтому geopandas с пакетом rtree делает join по миллионам точек за секунды.

import geopandas as gpd

points = gpd.read_file("incidents.geojson")
zones = gpd.read_file("districts.geojson")

# присоединить атрибуты района к каждой точке
joined = gpd.sjoin(points, zones, predicate="within", how="left")
print(joined.groupby("district").size())

Частые ошибки

• Разные проекции у слоёв. sjoin требует одинаковый CRS; иначе все совпадения «потеряются».
• Точка на границе двух полигонов. Может присоединиться к обоим или ни к одному — границы заранее продумывают.
• Забыть про индекс. Без R-дерева join больших слоёв упирается в $O(m \cdot n)$ и «висит».

Spatial join в реальной аналитике

Пространственное соединение — рабочая лошадка любого отчёта «по территориям». Город получает поток заявок 112 как точки с координатами, но мэру нужна сводка по районам: spatial join проставляет каждой заявке район, и дальше обычная группировка даёт таблицу «район — число заявок». Ритейл соединяет чеки (привязанные к магазинам) с переписными участками, чтобы понять, из каких кварталов приходят покупатели. Телеком джойнит координаты вышек с зонами покрытия и считает, сколько домохозяйств остаётся без связи. Во всех случаях соединение идёт не по совпадению ключа, а по геометрии, и именно это делает spatial join незаменимым там, где общего идентификатора между наборами попросту нет.

Важно понимать варианты соединения, потому что от них зависит результат. Для точек в полигонах берут предикат within, и тут возникает вопрос дубликатов: если точка попала на стык двух районов или в перекрывающиеся зоны, она присоединится к нескольким записям и «размножит» строки — это надо осознанно обрабатывать. Для линий, пересекающих область, берут intersects. А когда прямого попадания нет (остановку надо привязать к ближайшей станции метро), используют nearest — соединение по минимальному расстоянию. Выбор предиката — это и есть постановка вопроса, и ошибка здесь тихо искажает всю последующую статистику.

Итог

• Spatial join соединяет слои по пространственному отношению, а не по ключу.
• Самый частый случай — точки внутри полигонов (within) с последующей агрегацией.
• Наивно это $O(m \cdot n)$; R-дерево ускоряет почти до $O(m \log n)$.
• Оба слоя должны быть в одной системе координат.