Shapely: геометрические объекты и операции

Введение

Shapely — это специализированная библиотека с открытым исходным кодом для создания, анализа и модификации планарных геометрических объектов на языке Python. Она построена вокруг проверенного промышленного C-движка вычислительной геометрии GEOS, на который также опираются настольная QGIS, расширение PostGIS и библиотека GeoPandas.

Каждая геометрическая запись в таблицах GeoPandas на базовом уровне является объектом Shapely. Поэтому владение Shapely дает прямое понимание того, как на уровне алгоритмов вычисляются буферные зоны, пространственные пересечения или топологические проверки вхождения точек в полигоны. Библиотека устанавливается командой pip install shapely, а при установке GeoPandas подключается автоматически.

Типы геометрий

Shapely строго следует спецификации OGC Simple Features, предоставляя базовые классы: Point (точка), LineString (ломаная линия) и Polygon (полигон). Их составные коллекции — MultiPoint, MultiLineString и MultiPolygon — описывают составные пространственные объекты, состоящие из нескольких разрозненных частей (например, островной архипелаг или разделенный район).

from shapely import Point, LineString, Polygon

quduq = Point(420350, 4575120)
quvur = LineString([(420000, 4575000), (420600, 4575300)])
uchastka = Polygon([(420200, 4575000), (420500, 4575000),
                    (420500, 4575250), (420200, 4575250)])

print(uchastka.area, quvur.length)   # 75000.0  670.82...

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

Для текстовой сериализации и переносимости геометрий используется открытый стандарт WKT (Well-Known Text): инструкция print(quduq.wkt) выдаст строку POINT (420350 4575120). Эта запись нативно интерпретируется SQL-запросами PostGIS и калькулятором полей QGIS, выступая универсальным мостом обмена геометрией между программами.

Базовые типы геометрий Shapely (Point, LineString, Polygon) и их составные мультипримитивы

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

Пространственный предикат (spatial predicate) — это булева функция, проверяющая топологическое отношение между двумя геометриями по правилу «да/нет». Наиболее часто используются:

Предикат	Смысл вопроса
`intersects`	Имеют ли объекты хотя бы одну общую точку?
`contains` / `within`	Находится ли один объект целиком внутри другого?
`touches`	Соприкасаются ли объекты исключительно по границе?
`disjoint`	Полностью ли разобщены объекты в пространстве?

print(uchastka.contains(quduq))     # True
print(quvur.intersects(uchastka))   # True

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

Пространственные отношения между двумя геометриями: пересекаются, содержат, касаются и не пересекаются

Геометрические операции

В отличие от предикатов, геометрические операции конструируют новые векторные формы:

intersection — вычисление общей области двух объектов (пересечение);
union — слияние контуров в единый объект (объединение);
difference — геометрическое вычитание второго объекта из первого (разность);
buffer(d) — построение эквидистантной полосы на удалении d (буферная зона).

from shapely import box

zona = quduq.buffer(100)                  # санитарная зона 100 м
tegishli = uchastka.intersection(zona)    # часть участка внутри зоны
print(round(tegishli.area))               # 31365

qoshni = box(420500, 4575000, 420700, 4575250)   # смежный участок
print(uchastka.touches(qoshni))                  # True
print(uchastka.union(qoshni).area)               # 125000.0
print(round(uchastka.difference(zona).area))     # 43635

Буферный контур представляет собой не идеальную аналитическую окружность, а аппроксимирующий многоугольник: по умолчанию четверть окружности формируется 16 линейными сегментами. Поэтому площадь составляет 31 365 м² вместо теоретических π·100² = 31 416 м²; повышенная гладкость задается параметром quduq.buffer(100, quad_segs=64).

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

Геометрические операции над полигонами: пересечение, объединение, разность и построение буферной зоны вокруг точки

Работа с массивами объектов

В производственных задачах обрабатываются тысячи объектов. Начиная с релиза Shapely 2.0 (декабрь 2022 года), функции векторизованы на уровне C-массивов: вызов shapely.contains(uchastkalar, quduq) сканирует массив объектов в скомпилированном цикле в десятки раз быстрее стандартного цикла Python.

Для масштабных массивов прямое попарное сравнение всё равно неэффективно. В таких сценариях применяют пространственный индекс (spatial index) — структуру STRtree. Она индексирует ограничивающие прямоугольники (BBOX) объектов в иерархическое дерево, быстро возвращая на этапе грубого фильтра только потенциальных кандидатов:

import shapely
from shapely import STRtree

uchastkalar = [uchastka, qoshni]
print(shapely.contains(uchastkalar, quduq))    # [ True False]

daraxt = STRtree(uchastkalar)
idx = daraxt.query(zona, predicate="intersects")
print(idx)                                     # [0]
tanlangan = [uchastkalar[i] for i in idx]

Метод query возвращает не копии объектов, а их целочисленные индексы в исходном списке, по которым извлекаются строки таблицы атрибутов.

Именно этот движок задействован внутри метода sjoin библиотеки GeoPandas. Понимание индексов Shapely позволяет разработчику прогнозировать вычислительную сложность пространственных алгоритмов.

Системы координат и невалидная геометрия

Shapely не содержит информации о системах координат: библиотека оперирует абстрактными планарными числами. Если передать координаты в градусах, buffer(100) сгенерирует зону радиусом 100 градусов, охватывающую континент, а площадь рассчитается в квадратных градусах. Поэтому перед вычислением метрических буферов данные трансформируют в прямоугольную проекцию — например, UTM 42N (EPSG:32642) для центрального Узбекистана. Для этого используют библиотеку pyproj:

from pyproj import Transformer
from shapely.ops import transform

t = Transformer.from_crs("EPSG:4326", "EPSG:32642", always_xy=True)

nuqta_grad = Point(69.2401, 41.2995)            # Ташкент в градусах
xato = nuqta_grad.buffer(500)                   # радиус 500 градусов — ошибка
print(xato.bounds[0])                           # -430.7599 (долгота меньше -180)

nuqta_m = transform(t.transform, nuqta_grad)    # спроецировано в метры
print(round(nuqta_m.buffer(500).area))          # 784137 м²

Второй критический фактор — невалидная геометрия: самопересекающиеся контуры или дублирующиеся вершины приводят к фатальным ошибкам топологических алгоритмов. Проверка выполняется свойством is_valid, а функция make_valid из документации Shapely 2 исправляет сломанную топологию:

from shapely import is_valid_reason, make_valid

kapalak = Polygon([(0, 0), (2, 2), (2, 0), (0, 2)])   # самопересекающийся полигон («бабочка»)
print(kapalak.is_valid)              # False
print(is_valid_reason(kapalak))      # Self-intersection[1 1]

tuzatilgan = make_valid(kapalak)
print(tuzatilgan.geom_type, tuzatilgan.area)   # MultiPolygon 2.0

У исходного самопересекающегося контура площадь area равна 0, поскольку треугольники имеют противоположную ориентацию обхода вершин и математически уничтожают друг друга. После применения make_valid объект разделяется на корректный MultiPolygon с правильной суммарной площадью.

Заключение

Shapely образует фундаментальный каркас вычислительной геометрии в Python: библиотека создает векторные примитивы, оценивает пространственные отношения предикатами и генерирует новые производные контуры.

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