**Распределение стационарных торговых объектов по округам Москвы**

Хочу поделиться примером‑инструкцией как получить инсайты из геоданных без регистрации, смс и ML (только open‑source и бесплатные инструменты: OSM, python, Портал открытых данных Правительства Москвы, DataLens). Будь то точки продаж или очаги заболеваемости covid, прежде всего нам нужно понять общую картину и увидеть аномалии.

В рамках примера мы:

добудем очертания округов Москвы и данные по точкам продаж
создадим из точек сетку полигонов для отображения и анализа концентрации
посмотрим как соотнести полигоны с точками по координатам
отобразим данные в BI (DataLens) (базовый дашборд)

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

1. Получаем данные об очертаниях Москвы

Используем библиотеку osmnx и получаем данные о границах на уровне административных округов в формате geopandas.geodataframe.GeoDataFrame

df = ox.geometries_from_place('Moscow, Russia',tags={'admin_level': '5'})
df = df[df['element_type']=='relation'] #оставляем только полигоны

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

# с сохранением неразрывных границ между собой
df['wkt_simplified'] = df['geometry'].simplify(tolerance=0.0001,
                                                preserve_topology=True)

Во сколько раз упрощается геометрия для загрузки в BI

Округ	Изначальная длина WKT	Упрощенный WKT	Экономия
Северо-Западный АО	60284	10144	5.9
Зеленоградский АО	34815	7295	4.8
Центральный АО	19488	4316	4.5
Юго-Западный АО	42534	10021	4.2
Северо-Восточный АО	17031	5156	3.3
Новомосковский АО	79806	26605	3.0
Южный АО	20588	7343	2.8
Юго-Восточный АО	17816	6625	2.7
Северный АО	19695	7392	2.7
Восточный АО	18636	7971	2.3
Западный АО	71960	31793	2.3
Троицкий АО	35132	22929	1.5

Вторым действием мы переводим геометрию из WKT (well-known text) в формат, понятный для BI, где каждая пара координат в полигоне выделена квадратными скобками и разделена запятой (DataLens,Qlik). Также необходимо поменять местами широту и долготу.

Из WKT в формат BI

def bi_format(x):
    '''функция для перевода из wkt в формат для BI'''
    if 'MULTIPOLYGON' in x:
        return x.replace('MULTIPOLYGON','')\
        .strip()\
        .replace(', ','],[')\
        .replace('(','[')\
        .replace(')',']')\
        .replace(']]],[[[',']],[[')\
        .replace(' ',',')
    elif 'POLYGON' in x:
        return '[' + x.replace('POLYGON','')\
        .strip()\
        .replace(', ','],[')\
        .replace('(','[')\
        .replace(')',']')\
        .replace(']]],[[[',']],[[')\
        .replace(' ',',') + ']'
    return x
  
#добавялем колонку с форматом полигона для BI
df['wkt_simplified_bi'] = df['wkt_simplified'].apply(bi_format) 

#меняем местами широту и долготу для формата BI
df['wkt_simplified_bi'] = df['wkt_simplified_bi']\
  .apply(lambda x: str([[lvl2[::-1] for lvl2 in lvl1] for lvl1 in eval(x)]))

В результате формируем в БД (PostgreSQL) таблицу округов Москвы "msc_polygons" с исходной геометрией WKT, osmid, упрощенными полигонами в формате для BI и их названиями.

2. Забираем данные о торговых точках

Для примера я выбрал Стационарные торговые объекты на портале открытых данных Москвы, но там есть еще много интересных подборок для анализа.

В результате выгрузки мы получаем 65882 записей с координатами (сейчас мы не используем другие данные, в т.ч. и существующее справочное деление по округам).

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

Поэтому, после создания таблицы "msc_points" (global_id,широта, долгота) мы добавляем колонки с округленными значениями широты и долготы для агрегации их в сетку.
Долготу округляем до сотых, а широту до тысячных с шагом в 0.005, чтобы получить не вытянутый прямоугольник с севера на юг, а более квадратную форму.
Данный метод требует меньше вычислительных мощностей для составления сетки, чем, например, сопоставление точек с сотами H3: Uber’s Hexagonal Hierarchical Spatial Index.

Добавление столбцов широты, долготы для сетки

ALTER TABLE geo_project.msc_points ADD COLUMN grid_lat float;
ALTER TABLE geo_project.msc_points ADD COLUMN grid_lon float;
UPDATE geo_project.msc_points SET grid_lat = (floor(lat*100) +
case 	when lat*100 - floor(lat*100) <0.25 then 0
	when lat*100 - floor(lat*100) <=0.75 then 0.5
	when lat*100 - floor(lat*100) >0.75 then 1 end)/100;        
UPDATE geo_project.msc_points SET grid_lon = round(lon::numeric,2);

3. Сопоставление исходных точек с округами Москвы

В примере я пользуюсь PostgreSQL так как с ростом количества данных снижается удобство пользования python, появляются проблемы с нехваткой оперативной памяти.
Решение с использованием sql позволит выполнить задачу на кластере с бОльшим количеством ресурсов в PostgreSQL (или практически без изменения синтаксиса можно воспользоваться Spark SQL) .

Для создания справочника соотнесения округов Москвы (osmid) с каждой точкой (global_id) можно воспользоваться стандартными функциями:

st_geomfromtext(geometry) polygon_geom --гео формат для полигона округа
st_point(lon,lat) as point --гео формат для координаты точки
находим произведение таблицы с округами и таблицы с точками, где есть пересечения where st_intersects(point,polygon_geom))

Запрос целиком

with poly as (
	select
	osmid,
	st_geomfromtext(geometry) polygon_geom --преобразуем в гео формат
	from
	geo_project.msc_polygons)
--
,points as (
	select
	st_point(lon,lat) as point, --преобразуем в гео формат
	global_id
	from
	geo_project.msc_points)
--
,duplicates as (
    select 
    osmid,
    global_id,
    row_number() over(partition by global_id order by global_id) rn
    from points,poly 
    where st_intersects(point,polygon_geom)) --произведение множеств при условии пересечения точек с полигоном
select 
osmid,
global_id
from duplicates
where rn = 1

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

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

row_number() over(partition by global_id order by polygon_area asc) rn

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

4. Еще немного обработки данных для BI

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

bi_geom - элемент сетки, до которой агрегируются точки продаж в виде полигона
wkt_simplified_bi - геометрия полигона округа Москвы
grid_lat - широта центра ячейки сетки
grid_lon - долгота центра ячейки сетки
name - название округа
sp_count - количество точек (с группировкой до ячейки сетки)

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

SQL для формирования полигона сетки

'[['||
'['||text(round(   (grid_lat::numeric-0.0025)::numeric    ,3)) ||', ' || text(round(   (grid_lon::numeric-0.005)::numeric    ,3))||'], '
'['||text(round(   (grid_lat::numeric+0.0025)::numeric    ,3)) ||', ' || text(round(   (grid_lon::numeric-0.005)::numeric    ,3))||'], '
'['||text(round(   (grid_lat::numeric+0.0025)::numeric    ,3)) ||', ' || text(round(   (grid_lon::numeric+0.005)::numeric    ,3))||'], '
'['||text(round(   (grid_lat::numeric-0.0025)::numeric    ,3)) ||', ' || text(round(   (grid_lon::numeric+0.005)::numeric    ,3))||']'
||']]'

Т.е. создаем четыре точки, отталкиваясь от центра.
Из grid_lat=55.215, grid_lon=36.97
получаем '[[[55.213, 36.965], [55.218, 36.965], [55.218, 36.975], [55.213, 36.975]]]'

5. Очень короткая инструкция по созданию гео виджета

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

Скучная инструкция по созданию виджета

Инструкция работает, если уже создан аккаунт для работы.

На странице вверху справа через синюю кнопку создаем "воркбук"
Добавляем подключение (PostgreSQL или csv).
Нажимаем синюю кнопку вверху справа "создать подключение" и даем ему имя.
Нажимаем кнопку "создать датасет".
В разделе поля меняем тип данных для bi_geom, wkt_simplified_bi на "Геополигон"
Нажимаем "сохранить" и "создать чарт".
На странице создания чарта воспроизводим такую картинку

, где "Кол-во точек" = sum([sp_count]) добавляется через значок "+", а в цветах настраиваем нужный нам градиент.

6*. Определяем аномалии для каждого округа

Для Москвы, разбитой по округам, определить аномалии не составляет труда, но если мы пытаемся проанализировать всю Россию, то мы сталкиваемся с большими объемами данных и разными порядками цифр.

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

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

Для каждого полигоны определяется (значение 75-ого перцентиля + IQR * coef).
Если это значение ниже, чем количество точек внутри ячейки сетки, то считаем, что это аномальное скопление.
coef определяем следующим образом: начинаем с 5, но если у нас выбросов (ячеек выше границы) меньше чем [количество ячеек внутри округа / 100 + 1], то уменьшаем значение coef на 0.1 (но не более, чем до 1.5).

Реализация на python

dft = pd.read_sql("select * from geo_project.bi_dataset",con=mcfg.conn2())
dft = dft[~dft['name'].isna()] # убираем точки, которые не определились внутри полигонов по координатам

conc_points = [] #список датафреймов с выбросами по округам
i=1

for district in dft['name'].unique():

    coef = 5
    df = dft[dft['name']==district]
    print(f"{i}/{len(dft['name'].unique())}  {district} {' '*20}",end='\r')
    i+=1

    while coef>=1.5: #выбросом мы считаем iqr*1.5, не меньше
        lst = list(sorted(df['sp_count'])) 
        data = np.array(df['sp_count'])
        q3,q1 = np.percentile (data, [75 ,25]) #находим значения 75 и 25 персентилей внутри округа
        iqr = q3-q1 #находим интерквартильный размах (IQR) внутри округа
        if len(df[df['sp_count']>=iqr*coef]) < round(len(df)/100+1): #нам нужно, чтобы в выборку попало не менее 1% от общего числа квдратов + 1 (таким образом мы можем применять этот алгоритм к разным населенным пунктам)
            coef-=0.1
        else:
            df2 = df[df['sp_count']>=iqr*coef] #делаем выборку из квадратов-выбросов
            df2['coef'] = df2['sp_count'].apply(lambda x: round(coef,1)) 
            conc_points.append(df2)
            break
conc_df = pd.concat(conc_points) #объединяем все выбросы в одну таблицу
conc_df['conc_flg'] = conc_df['sp_count'].apply(lambda x: 1)

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

Особенности подхода

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

Спасибо за внимание!

Ссылки к статье:
Дашборд с результатом
Репозиторий с кодом подготовки данных

HowTo: базовая геоаналитика