Привет, хабр! Сегодня мы хотели бы продолжить тему обработки пространственных данных средствами Python библиотеки estaty. Мы уже рассказывали о том как можно Объединять открытые данные Open Street Map и Landsat для уточнения площадей зеленых зон вокруг объектов недвижимости.

Теперь поговорим о другом типе анализа: proximity analysis или обеспеченность объектов (или доступность) какими-либо другими полезными объектами, например теми же парками, больницами, детскими садами, и др. (Рисунок 1). 

Немного про анализ доступности

Предлагаем начать рассказ с небольшого литературного обзора. Так, например, анализ доступности, или зоны доступности, широко используются в геомаркетинге. Тема довольно сильно переплетена с анализом доступности каких-либо сервисов для выбранной локации, когда например сети магазинов хочется открыть новую точку. Не могу не сослаться на хорошую (на мой скромный взгляд) статью Геоаналитика с помощью Python и открытых данных: пошаговое руководство на хабре в контексте этой тематики. Для более искушенных читателей предпочту сослаться на Как мы используем модель Хаффа для открытия новых магазинов.

Также подобный тип анализа довольно широко распространен в крупных агрегаторах вроде Яндекс.Недвижимости (см. вкладку “Транспортная доступность” или “Инфраструктура”) при покупке / аренде жилой недвижимости. Часто речь идет про транспортную доступность, однако спектр инфраструктурных объектов для анализа огромный, вплоть до площадок для выгула собак.

Библиотека (которую мы используем)

Теперь перейдем к более предметному разговору: про инструменты которые можно для этого использовать (помимо геоинформационных систем вроде QGIS или ArcGIS и проприетарных сервисов в качестве конечного потребителя продукта). 

Мы, как разработчики, хотели бы анализировать не агрегированные показатели, а более сырые данные при помощи удобных библиотек или сервисов (конечно желательно с открытым исходным кодом). Подходящих для нас инструментов мы не нашли, поэтому разрабатываем и используем в анализе нашу собственную open-source библиотеку написанную на языке Python. Называется она estaty и если хотите, то можете посмотреть страничку с документацией после прочтения статьи:

• https://estaty.readthedocs.io/en/latest/?badge=latest 

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

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

Анализ (который мы можем сделать)

Сначала установим библиотеку с помощью команды

pip install estaty==0.1.0

или, если вы используете poetry:

poetry add estaty==0.1.0

Теперь мы готовы использовать версию 0.1.0

И приступаем к реализации. Для примера будем производить оценку обеспеченности по следующему адресу: “Биржевая линия, 16, Санкт-Петербург, Россия, 199034” - координаты {широта: 59.944843895537566, долгота: 30.294778398601856} (привет бывшим коллегам из Университета ИТМО). Будем анализировать окрестности этого объекта в радиусе 2 км.

Для того, чтобы узнать сколько времени нам потребуется для того чтобы дойти до ближайшего бара (по данным OSM - тут стоит оговориться) напишем и запустим следующий код:

from estaty.analysis.action import Analyzer
from estaty.data_source.action import DataSource
from estaty.main import EstateModel
from estaty.preprocessing.action import Preprocessor

import warnings
warnings.filterwarnings('ignore')


def launch_proximity_analysis_for_bars_objects():
   point_for_analysis = {'lat': 59.944843895537566, 'lon': 30.294778398601856}
   osm_source = DataSource('osm', params={'category': 'bar'})

   osm_reprojected = Preprocessor('reproject', params={'to': 'auto'}, from_actions=[osm_source])

   analysis = Analyzer('distance', params={'network_type': 'walk', 'visualize': True, 'color': 'orange',
                                           'edgecolor': 'black', 'title': 'Bars'},
                       from_actions=[osm_reprojected])

   model = EstateModel().for_property(point_for_analysis, radius=2000)
   founded_routes = model.compose(analysis)

   print(founded_routes.lines)
   print(f'Min length: {founded_routes.lines["Length"].min():.2f}, m')
   print(f'Mean length: {founded_routes.lines["Length"].mean():.2f}, m')
   print(f'Max length: {founded_routes.lines["Length"].max():.2f}, m')


if __name__ == '__main__':
   launch_proximity_analysis_for_bars_objects()

Здесь происходит следующее - мы конструируем пайплайн для анализа объекта недвижимости. В данном случае пайплайн обработки будет выглядеть как граф состоящий из последовательных преобразований над данными: 

Рассмотрим подробнее что происходит “под капотом”. Начнем с того, что архитектура estaty построена на 5 абстрактных слоях (в случае с пайплайном выше используются узлы только трех типов):

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

• Preprocessor - предобработка векторных или растровых данных, например назначение новой проекции;

• Merger - объединение данных требуемым образом. Возможные комбинации: векторные данные с векторными, растровые с растровыми и векторные с растровыми;

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

• Report - опциональный завершающий анализ модуль, который позволяет формировать отчет в удобном для пользователя формате.

Комбинируя эти узлы определенным образом можно составлять различные пайплайны для анализа данных без модификации способов обработки исходных данных - подменяя исключительно исходные источники или способы анализа (см. Анимация 1).

Итак, на рисунке 2 показано, что при анализе используется только 3 узла - это довольно простой пайплайн. В результате его работы мы получим geopandas GeoDataFrame (таблицу с геометриями объектов) с линейными объектами - маршрутами до объектов нашего интереса. В данном случае, - до баров. Визуализация результата, к слову, будет выглядеть следующим образом: 

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

• Маршрут с минимальной длиной: 450.11, м

• Средняя длина маршрутов: 2790.71, м

Как истинный эксперт в расположении баров в Санкт-Петербурге, могу ответственно заявить, что баров в области анализа, в действительности, конечно больше. Но мы сразу договорились в начале, что анализ будем делать на основе данных OSM, так что…

Сам пайплайн получился при этом независимым от источника данных, на которых производится анализ. Это обеспечивается за счет слабосвязности модулей библиотеки. Начнем с DataSource узла - в нем встроен механизм приведения данных к одному типу (например, векторные данные тут будут в формате geopandas GeoDataFrame и всегда будут только трех видов неважно какой источник мы для этого использовали, все они будут приведены к точкам, линиям, или полигонам). Затем следует предобработчик, который автоматически определит подходящую метрическую проекцию для наших данных. Для предобработчика важно лишь одно - данные должны прийти из DataSource узла, об остальном должен позаботится предшествующий ему узел. И завершает наш небольшой пайплайн блок анализа, который будет работать на любых геометрических объектах неважно площадные они, линейные или точечные. 

Таким образом, точно такие же расчеты мы можем делать для других категорий данных, которые мы будет доставать из OSM (естественно, лишь только OSM список не ограничен), например анализируя доступность школ, парков, водных объектов, туалетов, мусорных баков, отделений полиции и так далее - да чего угодно (какие типы объектов можно определить по тегам смотрите в вики страничке Map features или в документации estaty про источники данных): 

За счет чего подход можно расширять 

Мы упомянули, что анализ можем проводит на пространственных данных любой природы: социальные это данные, природные объекты или искусственные постройки. Заглянем ещё глубже в исходный код. Из рисунков видно, что рассчет расстояний производится на маршрутах построенных с учетом дорог, что логично, ведь обычно мы до баров не летаем вертолетами, а ходим пешком. Поэтому ‘distance’ анализ включает в себя поиск маршрутов на графе дорог, используя таковой из OSM. То есть алгоритм осуществляет поиск маршрута на следующих абстракциях (мы опирались на следующую статью, очень рекомендуем - OSMnx: Python for Street Networks): 

аф представляет собой географически привязанные (имеют атрибуты широты и долготы) узлы и ребра. Чтобы найти кратчайший маршрут от объекта A до объекта B, нужно сделать три шага:

1. Найти ближайший к объекту A узел на графе (назовем его узел 1) 

2. Найти ближайший к объекту B узел на графе (назовем его узел 2)

3. При помощи алгоритма обхода графа найти оптимальный маршрут от узла 1 к узлу 2

Если хочется узнать расстояние, то придется сложить длину маршрута от узла 1 к узлу 2 с расстояниями от объекта A к узлу 1 и от объекта B к узлу 2. Последние две компоненты можно рассчитывать по прямой так как сеть с узлами довольно густая. Но и тут есть трудность - если объект A и B точки, то найти ближайшие к ним узлы относительно просто (рассчитать расстояния от точки до точки). Если же геометрия объекта A например площадная, то придется искать “ближайший” узел к полигону. А что делать в случае линейного объекта. 

Далее ответим на вопрос “А как считать расстояния от узлов графа до векторных объектов разного типа: линий, точек и полигонов унифицированным образом?”. Для этого внутри библиотеки перед началом расчета маршрутов осуществляется преобразование объектов к точечному типу (мы ведь помним что находить расстояние между точками очень просто). Мы называем это репрезентациями (Рисунок 6).

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

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

Более продвинутый анализ 

Одними маршрутами до баров никого не удивить. Поэтому усложним задачу. Предположим, что мы хотим сделать более сложный анализ с объединением данных из различных источников. Например, мы хотим рассчитать сколько в среднем нам придется идти до парков, в которых растут дубы (ну любим мы после баров погулять и полюбоваться на раскидистые деревья). Для этого данные объединять будем из двух источников:

• Выгрузим векторные данные из GBIF | Global Biodiversity Information Facility (см. https://www.gbif.org/, на 19 октября 2023)

• OpenStreetMap - нативная интеграция выгрузки данных в библиотеке

Данные GBIF представляют собой векторный слой с точечными геометриями, которые показывают места, где были обнаружены те или иные виды живых организмов. Нас в данном случае будет интересовать Дуб черешчатый (или Quercus robur на латыни). 

Ссылка, цитирующая эти данные: GBIF.org (19 October 2023) GBIF Occurrence Download  https://doi.org/10.15468/dl.fzzxam. Так будут выглядеть наши данные, когда мы наложим точечный слой с дубами поверх полигонов парков:

Код для анализа будет выглядеть следующим образом: 

from estaty.analysis.action import Analyzer
from estaty.data_source.action import DataSource
from estaty.main import EstateModel
from estaty.merge.action import Merger
from estaty.preprocessing.action import Preprocessor
from estaty.report.action import Report

import warnings
warnings.filterwarnings('ignore')


def launch_parks_with_quercus_proximity_analysis():
   """
   Calculate proximity to parks using OpenStreetMap data and
   GBIF.org (19 October 2023) GBIF Occurrence Download  https://doi.org/10.15468/dl.fzzxam
   """
   point_for_analysis = {'lat': 59.944843895537566, 'lon': 30.294778398601856}
   # 1 Stage - define data sources and get data from them
   osm_source = DataSource('osm', params={'category': 'park'})
   bio_source = DataSource('csv', params={'path': './data/quercus_spb.csv',
                                          'lat': 'decimalLatitude', 'lon': 'decimalLongitude',
                                          'crs': 4326, 'sep': '\t'})

   # 2 Stage - re project layers into metric projection
   osm_reprojected = Preprocessor('reproject', params={'to': 'auto'}, from_actions=[osm_source])
   bio_reprojected = Preprocessor('reproject', params={'to': 'auto'}, from_actions=[bio_source])

   # 3 Stage - merge data into OSM polygons (method match) - order is important
   merged_vector = Merger('vector', params={'method': 'match', 'buffer': 10},
                          from_actions=[osm_reprojected, bio_reprojected])

   # 4 Stage - calculate distances to objects
   analysis = Analyzer('distance', params={'network_type': 'walk', 'visualize': True, 'color': 'green',
                                           'edgecolor': 'black', 'title': 'Parks with Quercus robur'},
                       from_actions=[merged_vector])

   # 5 Stage - display all prepared output calculations into console
   report = Report('stdout', from_actions=[analysis])

   # Launch model
   model = EstateModel().for_property(point_for_analysis, radius=2000)
   model.compose(report)

    
if __name__ == '__main__':
   launch_parks_with_quercus_proximity_analysis()

Из важного - отметим, что точки не обязательно должны точно попадать в полигон парка. Все таки мы объединяем данные из разных источников и они могут немного отличаться. Поэтому мы задаем буфер в 10 метров, и смотрим на результат (Рисунок 8) - будем считать, что в парке растет дуб, если полигон парка имеет пересечение с полигоном буфера по данным GBIF:

Ещё немного про актуальность

Тема анализа, как говорится, “на виду”. Следовательно, например функциональность построения маршрутов есть во многих решениях с исходным кодом (во время подготовки статьи наткнулся на любопытную тетрадку “Exploring routing tools for calculating paths for a set of origin-destination pairs” - как раз в тему), так и в популярных пользовательских сервисах вроде Google Maps или карт Яндекса, так и сервисах, связанных с картами опосредованно. Если же говорить про анализ обеспеченности, то тут тоже первопроходцем не получится стать. Можно например упомянуть такой инструмент как pandana. 

Все эти инструменты безусловно полезны, однако их использование в прикладных примерах требует написания некоторого количества вспомогательного кода. Мы же пытаемся сделать удобный инструмент для нашей собственной команды и хотим поделиться находками с сообществом. В приведенном в данной статье подходе мы же можем отметить независимость от источника данных и гибкий (на наш взгляд) способ эти векторные данные объединять. Будем рады комментариям и предложениям, оставляем несколько полезных ссылок: 

• https://github.com/red5ai/estaty - Python open-source библиотека для обработки пространственных данных и подготовки прототипов алгоритмов для анализа объектов недвижимости. Библиотека совсем молодая, но мы планируем её развивать и улучшать;

• https://estaty.readthedocs.io/en/latest/ - страничка с документацией;

• https://github.com/red5ai/estaty_examples - open-source репозиторий с примерами запуска для данной статьи

Данные, используемые в этом посте (& лицензии к ним): 

1. Map data copyrighted OpenStreetMap contributors and available from https://www.openstreetmap.org
   Лицензия — Open Data Commons Open Database License
   Ссылки актуальны на 25 октября 2023

2. GBIF.org (20 October 2023) GBIF Occurrence Download https://doi.org/10.15468/dl.f487j5
   Лицензия — Attribution-NonCommercial 4.0 International
   Ссылки актуальны на 25 октября 2023

С рассказом про анализ доступности с использованием данных OpenStreetMap и не только был Михаил Сарафанов и Wiredhut team

P.S. Если любопытно почитать эту же статью на английском, но на примере Берлина, то можете изучить Proximity Analysis to Find the Nearest Bar Using Python на Towards Data Science