Выбор кадастрового инженера с помощью Data Science10.06.2022 13:16

Для кого эта статья

Статья может быть полезна для:

Начинающих Data Scientist’ов (как я), чтобы посмотреть/попробовать инструменты интерактивной визуализации (включая геоданные) для РФ
Людей, кто столкнулся с земельными проблемами и пытается понять, что происходит в сфере деятельности кадастровых инженеров
Чиновников, отвечающих за деятельность Росрееста и кадастровых инженеров
Чиновников, отвечающих за развитие бизнеса и пытающихся понять, где есть «узкие места»

Автор: Предприниматель, начинающий Data Scientist на досуге.

Предыстория

Заканчивался 1 квартал 2020 года, ажиотаж вокруг пандемии ковид в РФ был на своем пике. Симптоматика первых переболевших показывала, что даже в случае относительно легко перенесенной болезни вопрос реабилитации и восстановления работоспособности (в том числе и психологическо-когнитивной) — встает на первое место. И мы наконец-то решили «Хватит сидеть, пора делать свое дело. Если не сейчас, то когда?!». В условиях повсеместной удаленки нашли иностранного профильного партнера-инвестора и разработали адаптированный к РФ концепт клиники/пансионата по реабилитации пациентов после перенесенного ковида. Ключевым риском для инвесторов была возможная скорость реализации проекта (после пандемии предполагалась реконцепция клиники в многопрофильный реабилитационный центр -, а это существенно большие инвестиции и сроки окупаемости) — поэтому было важно стартовать как можно быстрее.

Команда проекта была преисполнена энтузиазма, готова соинвестировать и мы договорились с инвесторами, что основной транш инвестиций пойдет не на стройку, а на расширение и оборудование приобретенных командой площадей. Мы достаточно быстро нашли несколько подходящих объектов в Московской области, но самым интересным показался объект, реализуемый Агентством по Страхованию Вкладов в рамках банкротство одного из банков РФ. Взвесив все «за» и «против», мы приняли решение об участии в публичных торгах и выкупили объект. Окрыленные победой на торгах, мы быстро заключили ДКП, произвели оплату и подали документы в Росреестр на регистрацию сделки. Не ожидая никаких подвохов с регистрацией (все-таки продавец — АСВ, торги — публичные, имущество — банковское) мы сразу же начали переговоры с подрядчиками по реновации и строительству. Как же мы тогда ошибались…

Сейчас, спустя уже полтора года после сделки, мы все еще не зарегистрировали право собственности на купленные объекты. Зато:

Пообщались (и получили от этого колоссальное количество негативного опыта) с десятками кадастровых инженеров.
Подали в Росреестр бессчётное количество межевых планов с исправлением их же ошибок, но на каждый получили бессодержательные отписки.
Площадь одного из участков Росреестр в одностороннем порядке уменьшил (!) на 2 гектара. Без уведомлений ни АСВ, ни нас как покупателя.
Ходим на заседания судов как на работу (за эту работу еще и платим юристам).
Познакомились с соседями, которые неожиданно «почуяв» деньги, подумали, что у них есть претензии территориального характера.

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

И мы подумали: ну не могут же быть все плохие (спойлер: по статистике — могут). Ну не может же Росреестр блокировать все решения (спойлер: по статистике — может), ведь из каждого утюга сейчас говорят о том, как же важно поддерживать бизнес и предпринимателей. А о каком бизнесе может идти речь, если ты даже не можешь зарегистрировать землю?
Наверное, мы просто неправильных кадастровых инженеров выбирали (мы поработали уже с 4я) — давайте найдем объективные данные и по ним выберем хорошего кадастрового инженера.

Данные Росреестра

Если зайти на сайт Росрееста и покопаться в его в глубинах можно найти реестр аттестованных кадастровых инженеров. Далее, по каждому кадастровому инженеру можно посмотреть основную информацию, членство в СРО, информацию о дисциплинарных взысканиях и, главное, — статистику его деятельности:

На момент написания статьи (Апрель-Май 2022) были доступны данные по 4 кв. 2021 года включительно. К сожалению, удобных методов выгрузить данные сайт не предоставляет. Поэтому на фриланс-бирже был найден профессиональный исполнитель, который всего за несколько дней (сайт Росреестра работает очень медленно, поэтому этот результат считаю очень хорошим смог собрать данные. Наверняка он без сна и отдыха ручками прокликал без малого 40 тыс. кадастровых инженеров в интуитивно понятном и дружественном интерфейсе сайта.

Приступаем к анализу

Импортируем необходимые библиотеки

import bokeh.io
import geopandas as gpd
import numpy as np
import pandas as pd
import pandas_bokeh
from bokeh.io import output_notebook, reset_output, show
from bokeh.models import ColumnDataSource, HoverTool, NumeralTickFormatter
from bokeh.palettes import Viridis3, Viridis256, viridis
from bokeh.plotting import ColumnDataSource, figure, output_file, save, show
from bokeh.resources import INLINE
from bokeh.transform import linear_cmap

bokeh.io.output_notebook(INLINE)

# Замьютим предупреждения от shapely и определим вывод графиков в ноутбук
import warnings
from shapely.errors import ShapelyDeprecationWarning
warnings.filterwarnings("ignore", category=ShapelyDeprecationWarning) 
output_notebook()

dt_dict = {
    "general_info" : {"path" :"./PARSED DATA/general.xlsx"},
    "statistics_1" : {"path" :"./PARSED DATA/statistics_1.xlsx"},
    "statistics_2" : {"path" :"./PARSED DATA/statistics_2.xlsx"},
    "sro_membership": {"path" :"./PARSED DATA/sro.xlsx"},
    "penalties": {"path" :"./PARSED DATA/discipline.xlsx"},
}

# Читаем данные, смотрим базовую информацию
for data_name, data_name_dict in dt_dict.items():
    data_path = data_name_dict.get("path")
    data_raw = pd.read_excel(data_path)
    data_name_dict["data_raw"] = data_raw
    display(data_raw.head(3))
    display(data_raw.info())

Предобработка данных

Необходимые шаги предобработки данных:

Таблица: Общая информация:

Сплит данных аттестата: att_number, att_date
Сплит ФИО: first_name, last_name , middle_name
reg_number: float → int
Коррекция типа данных для колонок с датами
Переименование колонок по словарю

Таблица: Членство в СРО

Таблица: Диспицлинарные взыскания

Коррекция типа данных для колонок с датами;
Переименование колонок по словарю;
Перевод в нижний регистр тип взыскания.

Таблица: Статистика:

Объединение 2х файлов статистики деятельности;
Создание колонки statistics_period из year + period (квартал) в формате дат pandas;
Переименование колонок по словарю.

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

# Определим словарь переименования
rename_columns_dict = {
    # Все таблицы
    "ID":"id",
    # Таблица дисциплинарных взысканий
    "Мера ДВ": "penalty_type",
    "Дата решения о применении меры ДВ": "penalty_decision_date",
    "Основание применения меры ДВ":	"penalty_decision_reason",
    "Дата начала ДВ": "penalty_start_date",
    "Дата окончания ДВ": "penalty_end_date",
    # Таблица членства в СРО
    "date_sro_incl" : "sro_inclusion_date",
    "date_sro_excl" : "sro_exclusion_date",
    "sro_excl_reason": "sro_exclusion_reason",
    # Таблица общей информации
    "date_added" : "added_date",
    "date_sro": "sro_date",
    "name": "full_name",
    # Таблица статистики
    "total_decisions": "decisions_total",
    "rejections_27fz": "decisions_27fz",
}

def clean_general_df(df_to_clean:pd.DataFrame) -> pd.DataFrame:
    """Function to clean raw general info. Returns cleaned df"""
    df_clean = (
        df_to_clean.copy()
        # Разбираем attestat на необходимые поля
        .assign(att_number = lambda x: x.attestat.str.split("_").str[0])#.str.split(" ").str[1])
        .assign(att_number = lambda x: x.att_number.str.split(" ").str[1])
        .assign(att_date = lambda x: x.attestat.str.split("дата выдачи: ").str[1])
        .drop("attestat", axis=1)
        # Разбираем ФИО. При такой реализации могут быть ошибки в нестандартных именах
        .assign(first_name = lambda x: x.name.str.split(" ").str[1])
        .assign(last_name = lambda x: x.name.str.split(" ").str[0])
        .assign(middle_name = lambda x: x.name.str.split(" ").str[-1])
        # Переименовываем колонки по словарю
        .rename(columns=rename_columns_dict)
        # Меняем формат данных 
        .assign(reg_number = lambda x: x.reg_number.astype("Int64"))
        # Меняем формат дат
        .assign(added_date = lambda x: pd.to_datetime(x.added_date, format="%d.%m.%Y", errors="ignore").dt.date)
        .assign(sro_date = lambda x: pd.to_datetime(x.sro_date, format="%d.%m.%Y", errors="ignore").dt.date)
        .assign(att_date = lambda x: pd.to_datetime(x.att_date, format="%d.%m.%Y", errors="ignore").dt.date)
    )
    return df_clean



def clean_sro_membership_df(df_to_clean:pd.DataFrame) -> pd.DataFrame:
    """Function to clean raw SRO membership info. Returns cleaned df"""
    df_clean = (
        df_to_clean.copy()
        # Переименовываем колонки по словарю
        .rename(columns=rename_columns_dict)
        # Меняем формат дат
        .assign(sro_inclusion_date = lambda x: pd.to_datetime(x.sro_inclusion_date, format="%d.%m.%Y", errors="coerce").dt.date)
        .assign(sro_exclusion_date = lambda x: pd.to_datetime(x.sro_exclusion_date, format="%d.%m.%Y", errors="coerce").dt.date)
    )
    return df_clean


def clean_penalties_df(df_to_clean:pd.DataFrame) -> pd.DataFrame:
    """Function to clean raw penalties info. Returns cleaned df"""
    df_clean = (
        df_to_clean.copy()
        # Переименовываем колонки по словарю
        .rename(columns=rename_columns_dict)
        # Меняем формат дат
        .assign(penalty_decision_date = lambda x: pd.to_datetime(x.penalty_decision_date, format="%d.%m.%Y", errors="coerce").dt.date)
        .assign(penalty_start_date = lambda x: pd.to_datetime(x.penalty_start_date, format="%d.%m.%Y", errors="coerce").dt.date)
        .assign(penalty_end_date = lambda x: pd.to_datetime(x.penalty_end_date, format="%d.%m.%Y", errors="coerce").dt.date)
        # Переводим в нижний регистр тип взыскания
        .assign(penalty_type = lambda x: x.penalty_type.str.lower())
    )
    return df_clean

def clean_statistics_df(df_to_clean:pd.DataFrame) -> pd.DataFrame:
    """Function to clean raw statistics info. Returns cleaned df"""
    df_clean = (
        df_to_clean.copy()
        # Переименовываем колонки по словарю
        .rename(columns=rename_columns_dict)
        # Создадим колонку с периодами деятельности
        .assign(statistics_period = lambda x: x.period.astype(str)+ "-"+ x.year.astype(str)) 
        .assign(statistics_period = lambda x: (pd.to_datetime(x.statistics_period, format="%m-%Y",errors="coerce") + pd.offsets.MonthEnd(0)).dt.date)
        .assign(quarter = lambda x: (x.period/3).astype("int64"))
    )
    return df_clean

# Сохраним очищенные данные в общий словарь
dt_dict["general_info"]["data_clean"] = clean_general_df(dt_dict["general_info"]["data_raw"])
dt_dict["statistics_1"]["data_clean"] = clean_statistics_df(dt_dict["statistics_1"]["data_raw"])
dt_dict["statistics_2"]["data_clean"] = clean_statistics_df(dt_dict["statistics_2"]["data_raw"])
dt_dict["sro_membership"]["data_clean"] = clean_sro_membership_df(dt_dict["sro_membership"]["data_raw"])
dt_dict["penalties"]["data_clean"] = clean_penalties_df(dt_dict["penalties"]["data_raw"])
dt_dict["statistics"] = {"data_clean": pd.concat([dt_dict["statistics_1"]["data_clean"], dt_dict["statistics_2"]["data_clean"], ])}

for k, v in dt_dict.items():
    display(v.get("data_clean").head(3))

Анализ

Анализ выданных аттестатов кад. инженеров

Посмотрим внимательнее на признак «Номер аттестата» att_number и попробуем понять, значат ли что-то цифры его составляющие. Больше всего мы бы хотели вытащить информацию о регионах выдачи аттестатов и, может быть, одна из цифр кодирует регион. Если это так, то регионов должно быть около 85, а максимальное количество инженеров ожидается в 50-м, 77-м, 78-м регионах.

_att_number_df = dt_dict["general_info"].get("data_clean")["att_number"]
_att_number_df = _att_number_df.str.split("-", expand=True)
_att_number_df = _att_number_df.dropna()
_att_number_df = _att_number_df.astype("int64", errors="raise")
_att_number_df = _att_number_df.rename(columns={0: "smt_0", 1: "smt_1", 2: "smt_2"})
display(_att_number_df["smt_0"].nunique())
display(_att_number_df["smt_1"].nunique())
display(_att_number_df["smt_2"].nunique())

83

7

1577

Отлично! Гипотеза пока не опровергнута.

Количество уникальных объектов первой части атт.номера близко кол-ву субъектов рф;
Вторая часть номера тоже представляет интерес: уникальных значений всего 7. Но что это может быть — пока не понятно.

Проверим количество аттестатов по предполагаемому признаку региона

# Посчитаем количество аттестатов по регионам
regions_att_count = (
    _att_number_df.groupby(by="smt_0")
    .agg(attestat_count=("smt_0", "count"))
    .sort_values(by="attestat_count", ascending=False)
).reset_index().rename(columns={"smt_0":"region"})
regions_att_count["rank_by_count"] = regions_att_count["attestat_count"].rank(ascending=False)
regions_att_count.head(5)

Визуализируем данные. Для образовательных целей данного проекта графики будут строиться с помощью библиотеки bokeh. Да, есть более удобные высокоуровневые библиотеки. Например, pandas-bokeh или hvplot. Последний даже мождет выступать в качестве бэкенда для графиков pandas вместо дефолтного matplotlib (pandas plotting backend docs). Однако хочется лучше понять bokeh и, в первую очередь, его возможности по более низкоуровневой кастомизации графиков.

source = ColumnDataSource(regions_att_count)

# Определяем цветовой mapper для раскраски в зависимости от количества аттестатов
mapper = linear_cmap(
    field_name="attestat_count",
    palette=Viridis256,
    low=min(regions_att_count.attestat_count),
    high=max(
        regions_att_count.attestat_count,
    ),
)

# Определим данные, показываемые при наведении на график
# Наведем красоту в форматах представления данных
TOOLTIPS = [
    ("Код региона", "@region"),
    ("Количество аттестатов", "@attestat_count{0.0a}"),
    ("Ранг по кол-ву аттестатов", "@rank_by_count{0o}"), #1 -> 1st, 2 -> 2nd
]

p = figure(
    plot_height=400,
    plot_width=1000,
    tooltips = TOOLTIPS,
    title = "Количество аттестатов по предполагаемому атрибуту региона"
)

p.vbar(
    x="region",
    top="attestat_count",
    color=mapper,
    source=source,
)

# Кастомизируем оси
p.xaxis.axis_label = "Код региона"
p.yaxis.formatter = NumeralTickFormatter(format='0a')
p.yaxis.axis_label = "Количество аттестатов"

show(p)

Интерактивный график доступен по ссылке

Гипотеза о том, что в номере аттестата закодирован регион выдачи, кажется, подтвердилась.

Наибольшее количество аттестатов выдано в Москве (77);
Удивительно, но на 2-м и 3-м месте, обгоняя Московскую Область (50), находятся Краснодарский край (23) и Республика Башкортостан (02).

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

# Скачаем данные о регионах из репозитория HFLabs. Спасибо ребятам за инфо в удобном формате
region_naming = pd.read_csv(
    "https://raw.githubusercontent.com/hflabs/region/master/region.csv",
    dtype=object,
)

# Достаем необходимые поля из таблицы регионов
geoname_df = region_naming.loc[:, ["kladr_id", "geoname_name", "iso_code"]]
geoname_df["code"] = geoname_df["kladr_id"].str[0:2]
#geoname_df["iso_code"] = geoname_df["iso_code"].str.replace("-", ".")
display(geoname_df.head(3))

# Смерджим данные в датафрейм с основной информацией
general_info_clean = dt_dict["general_info"].get("data_clean").copy()
general_info_clean["att_region"] = (
    general_info_clean["att_number"].str.split("-", expand=False).str[0]
)
general_info_clean = general_info_clean.merge(
    geoname_df[["geoname_name", "code", "iso_code"]],
    how='left',
    left_on="att_region",
    right_on="code",
).drop("code", axis=1)

# Посмотрим на результат и сохраним в словаре с данными
display(general_info_clean.head(3))
dt_dict["general_info"]["data_clean"]  = general_info_clean

Анализ статистики деятельности кад.инженеров

Проанализируем статистику деятельности кад.инженеров во времени:

Посчитаем суммарное количество отказов по всем причинам;
Посчитаем долю отказов.

(!) Так как обработка документов Росреестром растянуто во времени, могут быть кварталы, когда количество полученных в периоде отказов (по сути, по поданным ранее документам) превышает количество поданных в периоде документов

statistics_df = dt_dict["statistics"]["data_clean"]
cols_to_sum = ["decisions_27fz", "decisions_mistakes", "decisions_suspensions"]
statistics_df["rejections_total"] = statistics_df[cols_to_sum].sum(axis=1)
statistics_df["acceptions_total"] = statistics_df["decisions_total"] - statistics_df["rejections_total"]
statistics_df["rejections_share"] = (
    statistics_df["rejections_total"] / statistics_df["decisions_total"]
)
statistics_df["acceptions_share"] = (
    statistics_df["acceptions_total"] / statistics_df["decisions_total"]
)
display(statistics_df.head(3))

Посмотрим на агрегированную статистику отказов во времени. Так как отказы Росрееста получаются с временным лагом, возможна ситуация, когда доля отказов > 1. Данные по кварталам представлены накопленным итогом.

# Подготовим данные для Bokeh
_df = statistics_df.replace([np.inf, -np.inf], np.nan, inplace=False).dropna()
# Посчитаем суммарное количество принятых и отвергнутых документов по кварталам
period_stat_df = (
    _df[["acceptions_total", "rejections_total", "statistics_period"]]
    .groupby("statistics_period")
    .sum()
    .reset_index()
)
period_stat_df["rejections_share"] = period_stat_df["rejections_total"] / (
    period_stat_df["rejections_total"] + period_stat_df["acceptions_total"]
)
source = ColumnDataSource(period_stat_df)

# Определим данные, показываемые при наведении на график
# И наведем красоту в форматах отображения данных
hover_1 = HoverTool(
    tooltips=[
    ("Период", "@statistics_period{%m-%Y}"),
    ("Количество аксептов", "@acceptions_total{0.0a}"),
    ("Количество отказов", "@rejections_total{0.0a}"),
    ("Доля отказов", "@rejections_share{0.0%}"),
],
    formatters={
        "@statistics_period"        : 'datetime', # use 'datetime' formatter for '@date' field
    },
)

# Возьмем 2 цвета из палитры Viridis
colors = viridis(2)

# Определяем график
p = figure(
    width=1000,
    height=400,
    title="Поквартальное (накопленное за год) количество одобренных и отвегнутых документов",
    toolbar_location=None,
    x_axis_type="datetime",
    tools=[hover_1],
)

# Добавляем на график данные
p.vbar_stack(
    ["acceptions_total", "rejections_total"],
    # Ось Х - это ось времени, где базовая единица миллисекунда.
    # Поэтому ширину столбцов необходимо указывать достаточно большую
    width=5e9,
    x="statistics_period",
    color=colors,
    source=source,
)

# Кастомизируем названия осей
p.xaxis.axis_label = "Период"
p.yaxis.formatter = NumeralTickFormatter(format='0.0a')
p.yaxis.axis_label = "Количество документов"

show(p)

Интерактивный график доступен по ссылке

Мы видим, что с начала 2019 года явно поменялась структура и/или подход к проверке документов: при сохранении общей динамики и сезонности, количество отказов возрасло многократно. В вики ведомства ничего примечательного относительно 2018–2019 годов не написано, да и на TAdvisor тоже ничего примечательного в данные периоды.
Также достаточно странным выглядит околонулевой выброс отказов в 1 кв. 2020 года.

Посмотрим на динамику доли отказов

# Подготовим данные дополнив расчетом доли принятых документов
period_stat_df["acceptions_share"] = period_stat_df["acceptions_total"] / (
    period_stat_df["rejections_total"] + period_stat_df["acceptions_total"]
)

source_2 = ColumnDataSource(period_stat_df)

# Определим данные, показываемые при наведении на график
# И наведем красоту в форматах отображения данных
hover_2 = HoverTool(
    tooltips=[
    ("Период", "@statistics_period{%m-%Y}"),
    #("Количество аксептов", "@acceptions_total{0.0a}"),
    #("Количество отказов", "@rejections_total{0.0a}"),
    ("Доля отказов", "@rejections_share{0.0%}"),
    ("Доля акцептов", "@acceptions_share{0.0%}"),
],
    formatters={
        "@statistics_period"        : 'datetime', # use 'datetime' formatter for '@date' field
    },
)

# Возьмем 2 цвета из палитры Viridis
colors = viridis(2)

# Определяем график
p2 = figure(
    width=1000,
    height=400,
    title="Поквартальная (накопленная за год) доля одобренных и отвегнутых документов",
    toolbar_location=None,
    x_axis_type="datetime",
    tools=[hover_2],
)

# Добавляем на график данные
p2.vbar_stack(
    ["acceptions_share", "rejections_share"],
    # Ось Х - это ось времени, где базовая единица миллисекунда.
    # Поэтому ширину столбцов необходимо указывать достаточно большую
    width=5e9,
    x="statistics_period",
    color=colors,
    source=source_2,
)

# Кастомизируем названия осей
p2.xaxis.axis_label = "Период"
p.yaxis.formatter = NumeralTickFormatter(format='0%')
p2.yaxis.axis_label = "Доля документов"

show(p2)

Интерактивный график доступен по ссылке

Геовизуализация

Визуализируем данные на карте России. Для этого получим границы регионов с помощью OpenStreetMaps и его Overpass API для запросов. Большой аналитической ценности в данном этапе анализа нет — все можно увидеть в таблице, но так как цель данного анализа образовательная, то очень хотелось научиться визуализировать данные именно по РФ и реализовать и этот функционал.

Для тех, кто хотел бы чуть больше узнать о геоданных/геовизуализации крайне рекомендую архив курса Университета Хельсинки Henrikki Tenkanen and Vuokko Heikinheimo, Digital Geography Lab, University of Helsinki. Наверное, это один из лучших и комплексных мануалов, расказывающий (кратко, но по делу) весь процесс end-to-end

import pandas as pd
import requests
import geopandas as gpd
from osm2geojson import json2geojson

# Создадим запрос административных границ регионов
overpass_url = "http://overpass-api.de/api/interpreter"
overpass_query = """
    [out:json];
    rel[admin_level=4]
    [type=boundary]
    [boundary=administrative]
    ["ISO3166-2"~"^RU"];
    out geom;
    """
# Запрашиваем данные и формируем GeoDataFrame
response = requests.get(overpass_url, 
                        params={'data': overpass_query})
response.raise_for_status()
data = response.json()
geojson_data = json2geojson(data)
gdf_osm = gpd.GeoDataFrame.from_features(geojson_data) 

# Конвертируем словари тэгов ответа на запрос в колонки
df_tags = gdf_osm["tags"].apply(pd.Series)

# Определим, какие колонки оставить для дальнейшего анализа
# По сути - удалим переводы наименовая регионов на разные языки, оставив ru и en
cols_keep = []
for col in list(df_tags.columns):
    if "name:" not in col:
        cols_keep.append(col)
cols_keep.extend(["name:en", "name:ru"])

# Получим финальный геодатафрейм с нужными колонками
gdf_full = pd.concat([gdf_osm, df_tags.loc[:,cols_keep]], axis=1)
display(gdf_full.head())

Наиболее быстрый и простой метод построить интерактивную карту с Bokeh — воспользоваться более высокоуровневой библиотекой pandas_bokeh, которая поддерживает данные геоформата. Приведем пример и выведем интерактивную карту регионов РФ.

# Установим координатную систему Coordinate Reference System (CRS)
gdf_full_mercator = gdf_full.set_crs('epsg:4326')

gdf_full_mercator.plot_bokeh(
    figsize = (1000, 600),
    simplify_shapes=20000,
    hovertool_columns=["name:ru"],
    title="Пустая карта РФ",
    xlim=[20, 180],
    ylim=[40, 80],

)

Анализ общего количества документов

Подготовим статистику для визуализации и объединим с данными о границах регионов:

_df = statistics_df.replace([np.inf, -np.inf], np.nan, inplace=False).dropna()
_df_general = dt_dict["general_info"]["data_clean"]
_df = _df.merge(_df_general, how="left", on="id")
_df.head(3)

_df3 = (_df[["decisions_total","acceptions_total","rejections_total", "year", "att_region","iso_code"]]
.groupby(["year", "iso_code"])
.sum())

# Пересчитаем доли на агрегатах
_df3["rejections_share"] = (
    _df3["rejections_total"] / _df3["decisions_total"]
)
_df3["acceptions_share"] = (
    _df3["acceptions_total"] / _df3["decisions_total"]
)

annual_reg_stat = _df3.reset_index()
display(annual_reg_stat.head(3))

reg_stat_2021 = annual_reg_stat.loc[annual_reg_stat["year"]==2021,]
reg_stat_2021 = reg_stat_2021.replace("",np.nan).dropna()
points_to_map = gdf_full_mercator.merge(reg_stat_2021, how="left", left_on="ISO3166-2", right_on="iso_code")


#Replace NaN values to string 'No data'.
points_to_map.loc[:,["year","decisions_total","acceptions_total","rejections_total", "rejections_share"]].fillna('No data', inplace = True)
points_to_map.head()

plot_total_counts = points_to_map.plot_bokeh(
    figsize = (1000, 600),
    simplify_shapes=20000,
    hovertool_columns=["name:ru", "decisions_total","acceptions_total","rejections_total", "rejections_share"],
    dropdown = ["decisions_total","acceptions_total","rejections_total"],
    title="2021",
    colormap="Viridis",
    colorbar_tick_format="0.0a",
    xlim=[20, 180],
    ylim=[40, 80],
    return_html=True,
    show_figure=True,
)

# Export the HTML string to an external HTML file and show it:
with open("plot_total_counts.html", "w") as f:
    f.write((r"""""" + plot_total_counts))

Интерактивный график доступен по ссылке

Анализ отказов

plot_rejections = points_to_map.plot_bokeh(
    figsize=(1000, 600),
    simplify_shapes=20000,
    hovertool_columns=[
        "name:ru",
        "decisions_total",
        "acceptions_total",
        "rejections_total",
        "rejections_share",
    ],
    dropdown=["rejections_share"],
    title="2021",
    colormap="Viridis",
    colorbar_tick_format="0%",
    xlim=[20, 180],
    ylim=[40, 80],
    return_html=True,
    show_figure=True,
)

# Export the HTML string to an external HTML file and show it:
with open("plot_rejections.html", "w") as f:
    f.write(r"""""" + plot_rejections)

Интерактивный график доступен по ссылке

Визуализируем изменение доли отказов по регионам во времени. Ранее мы определили, что отказы «поперли» только с 2019 года. Соответственно статистику отобразим с этого момента. Для этого агрегируем данные по годам/регионам и подготовим dataframe в wide формате для возможности отображения на географике.

display(annual_reg_stat.head())
statistics_df_wide = annual_reg_stat.pivot(index="iso_code", columns=["year",])
# Убираем мультииндекс и объединяем название колонок с годами
statistics_df_wide.columns = ['_'.join((col[0],str(col[1]))) for col in statistics_df_wide.columns]
statistics_df_wide.reset_index(inplace=True)
statistics_df_wide.head()

# Replace NaN values to string 'No data'.
statistics_df_wide.fillna('No data', inplace = True)

# Combine statistics with geodataframe
history_to_map = gdf_full_mercator.merge(statistics_df_wide, how="left", left_on="ISO3166-2", right_on="iso_code")

#Specify slider columns:
slider_columns = ["rejections_share_%d"%i for i in range(2019, 2022)]

#Specify slider columns:
slider_range = range(2019, 2022)

plot_rejections_slider = history_to_map.plot_bokeh(
    figsize = (1000, 600),
    simplify_shapes=20000,
    hovertool_columns=["name:ru"]+slider_columns,  
    slider=slider_columns,
    slider_range=slider_range,
    slider_name="Year",
    title="Изменение доли отказов по регионам/годам",
    colormap="Viridis",
    colorbar_tick_format="0%",
    xlim=[20, 180],
    ylim=[40, 80],
    return_html=True,
    show_figure=True,
)

# Export the HTML string to an external HTML file and show it:
with open("plot_rejections_slider.html", "w") as f:
    f.write(r"""""" + plot_rejections_slider)

Интерактивный график доступен по ссылке

В 2020 году лидером по доле отказов была Астраханская область. «Зарезано» 90% поданных документов. В 2021 году в лидеры вырывается Московская область с 73% отказов.

Цифры колоссальные, если учесть сколько труда стоит за каждым из документов:

как минимум, несколько часов работы кадастрового инженера;
сходить в МФЦ (в лучшем случае) и подать их;
время сотрудников Росреестра на формирование отказа.

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

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

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

Анализ работы индивидуальных инженеров

# Объединим датасеты статистики работы и общей информации по кад.инженерам
kadeng_stat = statistics_df.copy()
_df_general = dt_dict["general_info"]["data_clean"]
kadeng_stat = kadeng_stat.merge(_df_general, how="left", on="id")

# Сгруппируем данные по кадастровым инженерам и годам
kadeng_stat_agg = (kadeng_stat[["decisions_total","acceptions_total","rejections_total", "year","id", "att_region"]]
.groupby(["att_region","id", "year"])
.sum())

# Пересчитаем доли на агрегатах
kadeng_stat_agg["rejections_share"] = (
    kadeng_stat_agg["rejections_total"] / kadeng_stat_agg["decisions_total"]
)
kadeng_stat_agg["acceptions_share"] = (
    kadeng_stat_agg["acceptions_total"] / kadeng_stat_agg["decisions_total"]
)

kadeng_stat_agg.replace([np.inf, -np.inf], np.nan, inplace=True)
kadeng_stat_agg = kadeng_stat_agg.reset_index(drop=False)

kadeng_stat_agg.loc[kadeng_stat_agg["year"] == 2021].describe()

decisions_total_hist = kadeng_stat_agg.loc[kadeng_stat_agg["year"] == 2021].plot_bokeh(
    kind="hist",
    bins = 100,
    y=["decisions_total"],
    xlim=(0, 3000),
    vertical_xlabel=True,
    show_average = True,
    title = "РФ_2021: Количество поданных документов",
    show_figure=False,

)

rejections_share_hist = kadeng_stat_agg.loc[kadeng_stat_agg["year"] == 2021].dropna().plot_bokeh(
    kind="hist",
    bins=np.arange(0, 3.5, 0.1),
    y="rejections_share",
    xlim=(0, 2),
    vertical_xlabel=True,
    show_average = True,
    title = "РФ_2021: Доля отказов",
    show_figure=False,
)

plot_kad_eng_stat = pandas_bokeh.plot_grid([[decisions_total_hist, rejections_share_hist]], width=400, height=300, return_html=True,)


# Export the HTML string to an external HTML file and show it:
with open("plot_kad_eng_stat.html", "w") as f:
    f.write(r"""""" + plot_kad_eng_stat)

Интерактивный график доступен по ссылке

В 2021 году средняя доля отказов в группировке по кадастровым инженерам составляет почти 25% — вдвое больше, чем доля отказов по суммарному количеству документов. Гипотеза: есть небольшое количество «супер-успешных» кадастровых инженеров, с большим количеством поданных документов, которые «проходят» на отлично и которые вытягивают среднюю статистику

Посмотрим аналогичную статистику по Московской области

def plot_hist_by_region(year, region_num, kadeng_stat_agg):

    if isinstance(region_num, int):
        region_num=str(region_num)

    _decisions_total_hist = kadeng_stat_agg.loc[((kadeng_stat_agg["year"] == year) & (kadeng_stat_agg["att_region"] == region_num))].plot_bokeh(
        kind="hist",
        bins = 30,
        y=["decisions_total"],
        #xlim=(0, 3000),
        vertical_xlabel=True,
        show_average = True,
        title = f"{region_num}_{year}: Количество поданных документов", #"РФ 2021: Количество поданных документов",
        show_figure=False,
    )

    _rejections_share_hist = kadeng_stat_agg.loc[((kadeng_stat_agg["year"] == year) & (kadeng_stat_agg["att_region"] == region_num))].dropna().plot_bokeh(
        kind="hist",
        #bins=np.arange(0, 2.5, 0.1),
        bins = 30,
        y=["rejections_share"],
        #xlim=(0, 2.5),
        vertical_xlabel=True,
        show_average = True,
        title = f"{region_num}_{year}: Доля отказов",
        show_figure=False,
    )
    return [_decisions_total_hist, _rejections_share_hist]

plots_list = plot_hist_by_region(2021, 50, kadeng_stat_agg)
plot_kad_eng_stat_50 = pandas_bokeh.plot_grid([plots_list], width=400, height=300, return_html=True,)


# Export the HTML string to an external HTML file and show it:
with open("plot_kad_eng_stat_50.html", "w") as f:
    f.write(r"""""" + plot_kad_eng_stat_50)

Интерактивный график доступен по ссылке

Для жителей Московской области или москвичей, кто хотел бы решить земельные вопросы, статистика неутешительная. «Средний» кадастровый инженер получил в 2021 году 98% отказов.

Определим лидеров и аутсайдеров среди кадастровых инженеров. Дальнейший анализ сделаем для данных по Московской области за последние 3 года (2019–2021). Нас интересует рэнкинг по количеству документов (больше — лучше) и доле отказов (больше — хуже).

years = [2019, 2020, 2021]
region_num = "50"

kadeng_stat_50 = kadeng_stat_agg.loc[((kadeng_stat_agg["year"].isin(years) ) & (kadeng_stat_agg["att_region"] == region_num))]

# Переведем в wide форму
kadeng_stat_50_wide = kadeng_stat_50.pivot(index=["att_region", "id"], columns=["year",])
# Убираем мультииндекс и объединяем название колонок с годами
kadeng_stat_50_wide.columns = ['_'.join((col[0],str(col[1]))) for col in kadeng_stat_50_wide.columns]
kadeng_stat_50_wide.reset_index(inplace=True)

# Дополним данными ФИО и аттестата кад.инженера
kadeng_stat_50_wide = kadeng_stat_50_wide.merge(general_info_clean, how="left", left_on="id", right_on="id")
kadeng_stat_50 = kadeng_stat_50.merge(general_info_clean, how="left", left_on="id", right_on="id")
display(kadeng_stat_50_wide.head())
display(kadeng_stat_50.head())

cols_to_plot = [
    "decisions_total",
    "acceptions_total",
    "rejections_total",
    "rejections_share",
    "full_name",
    "att_number",
    "att_date",
    "year",
]

#Use the field name of the column source
mapper = linear_cmap(field_name='year', palette=Viridis3, low=2019 ,high=2021)

source = ColumnDataSource(data=kadeng_stat_50.loc[:, cols_to_plot])

TOOLTIPS = [
    ("Год", "@year"),
    ("ФИО", "@full_name"),
    ("Всего решений", "@decisions_total"),
    ("Доля отказов", "@rejections_share{(0%)}"),
    ("Аттестат", "@att_number"),
]

p = figure(width=800, height=400, tooltips=TOOLTIPS, title="Количество решений Росреестра: всего и отказов",)

p.circle(
    "rejections_total",
    "decisions_total",
    source=source,
    color = mapper,
    legend_group = "year"
)

p.legend.location = "top_left"
p.xaxis.axis_label = "Количество отказов"
p.xaxis.formatter = NumeralTickFormatter(format='0')
p.yaxis.axis_label = "Количество решений"

# Output to file / notebook
reset_output()
#output_file("plot_kad_eng_50_2019_2021.html")
#save(p)
output_notebook()
show(p)

Интерактивный график доступен по ссылке

cols_to_plot = ["decisions_total", "acceptions_total", "rejections_total",  "rejections_share", "full_name", "att_number", "att_date", "year", ]
kadeng_stat_50.loc[:,cols_to_plot]

Выбор лучших кад.инженеров по Московской области

Определим лучших кадастровых инженеров по Московской области

kadeng_stat_50_wide.head(3)

# Определим правила ренкинга: для колонок где "меньше-лучше" используем параметр ascending = True
cols_to_rank = {
    "decisions_total": False,
    "acceptions_total": False,
    "rejections_share": True,
}
years_to_rank = [2019, 2020, 2021]

# Функция создания колонок с рангом по разным показателям/годам
def rank_kad_eng(df, cols_to_rank, years_to_rank):
    for year in years_to_rank:
        for col, asc_param in cols_to_rank.items():
            df[f"rank_{col}_{year}"] = df[f"{col}_{year}"].rank(
                na_option="bottom",
                ascending=asc_param,
                method="min",
            )
    return df


kadeng_stat_50_ranked = rank_kad_eng(kadeng_stat_50_wide, cols_to_rank, years_to_rank)

# Определим интегральный взвешенный показатель, учитывающий активность и доли отказов
# Веса определил исходя из своего видения важности критериев
# Для своих целей я решил 2019 не использовать, его вес будет 0

ranking_weights = {
    # Sum to 1
    "years_weights": {
        2019: 0,
        2020: 0.25,
        2021: 0.75,
    },
    # Sum to 1
    "indicator_weights": {
        "decisions_total": 0.25,
        "rejections_share": 0.75,
        "acceptions_total": 0,
    },
}


def integral_rank(df, ranking_weights):
    df["integral_rank"] = 0
    for year, year_weight in ranking_weights.get("years_weights").items():
        for indicator, ind_weight in ranking_weights.get("indicator_weights").items():
            df["integral_rank"] = (
                df["integral_rank"]
                + df[f"rank_{indicator}_{year}"] * year_weight * ind_weight
            )
    return df

kadeng_stat_50_integral = integral_rank(kadeng_stat_50_ranked, ranking_weights)


display(
    kadeng_stat_50_integral.sort_values(by="integral_rank").head(10).T
)

Предварительные выводы

Статистика по Московской Области крайне печальная:

Росреестр отказывает по ¾ поданных документов;
Среднестатистический кад. инженер в 2021 году в Московской области получил 98% отказов.

Я не буду размышлять и строить гипотезы о корупциогенной природе такого уровня отказов. В любом случае, такие цифры выглядят, как минимум, как саботаж нормальной работы и перекладывание ответственности с Росреестра на Суды (куда вынуждены идти люди, которым надо решить земельные вопросы). Возникают сомнения в целесообразности существования структуры (включая немаленькие бюджеты, в том числе на ИТ и сотрудников), чья суть деятельности сводится к практически 100% отказу.

Допущения и ограничения анализа

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

Следующие шаги

ML/Статистика:

Анализ второй цифры в аттестате и ее значимости влияния на на эффективность деятельности;
Анализ домена почты кад.инженера и его влияния на эффективность деятельности;
Анализ даты выдачи аттестата и ее влияния на эффективность деятельности;
Анализ СРО по размеру и эффективности деятельности;