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

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

Итак поехали…

Загрузка датасета

На сайте https://data.gov.ru/ находим и скачиваем датасет Объекты культурного наследия на территории Санкт-Петербурга. Последнее обновление датировано 2016 годом, но для объектов возрастом в пару веков это не будет проблемой.

Давайте оценим количество данных, с которыми будем работать:

import pandas
data = pandas.read_csv('spb_memo.csv')
print('Количество строк:', len(data))

Количество строк: 9275

Ого, более девяти тысяч объектов культурного наследия! Недаром наш город называют культурной столицей России. Что есть внутри?

data.head()

Выборка по районам

У нас есть названия объектов, даты постройки, адреса и даже имена авторов, отлично!

Как правило, городские достопримечательности, по большей части, располагаются ближе к центру города. Но любая гипотеза требует доказательств. Посмотрим статистку по районам:

districts = list(data['district'])
districts_unique = list(set(districts))
total_per_district = []

for district in set(districts):
    district_counter = 0
    for index in range(len(districts)):
        if districts[index] == district:
            district_counter += 1
    total_per_district.append(district_counter)

seaborn.barplot(x=total_per_district, y=districts_unique)

Никаких сюрпризов. Большинство памятников архитектуры располагаются в Центральном, Петродворцовом и Адмиралтейском районах города. Центральный и Адмиралтейский районы — это самый что ни на есть центр. На периферии объектов значительно меньше, за исключением Петродворцового и Пушкинского района. Но вряд ли за свои короткие выходные друг поедет так далеко, хотя и много потеряет. Как бы там ни было, вычёркиваем Пушкинский район, а также другие районы, находящиеся за КАДом.

outside_districts = ['Пушкинский', 'Кронштадтский', 'Кронштадт', 
    'Колпинский', 'Курортный', 'Петродворцовый', 'Санкт-Петербург']
districts_unique = [item for item in districts_unique 
    if item not in outside_districts ]

for district in districts_unique:
    district_counter = 0
    for index in range(len(districts)):
        if districts[index] == district:
            district_counter += 1
    total_per_district.append(district_counter)

seaborn.barplot(x=total_per_district, y=districts_unique)

Проверка других столбцов

Что ж, круг сужается, двигаемся дальше. В датасете есть такая характеристика объектов, как протекционная категория. Будет ли нам полезна эта колонка?

protection_categories = list(data['protection_category'])
protection_categories_unique = list(set(protection_categories))
total_per_category = []

for category in set(protection_categories):
    category_counter = 0
    for index in range(len(districts)):
        if protection_categories[index] == category:
            category_counter += 1
    total_per_category.append(category_counter)
    
seaborn.barplot(x=total_per_category, y=protection_categories_unique)

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

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

Выборка по авторам

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

authors_all = list(data['author'])
authors = []
total_per_author = []
    
for author_line in authors_all:
    if author_line == author_line:
        if ',' in author_line:
            for author in author_line.replace(';', ',').replace('арх. ', '').replace('худ. ', '').replace('гражд.инж. ', '').replace('архитекторы ', '').replace('фонтанный мастер ', '').replace('арх-ры ', '').split(','):
                author = author.strip()
                if author not in authors:
                    authors.append(author.strip())
                    total_per_author.append(1)
                else:
                    index = authors.index(author.strip())
                    total_per_author[index] += 1
        else:
            if author not in authors:
                authors.append(author.strip())
                
authors_df = pandas.DataFrame(authors, columns=['name'])
authors_df['count'] = total_per_author

seaborn.barplot(x=authors_df['count'], 
                y=authors_df['name'])

Этот график совсем не выглядит информативным. Сколько же там авторов:

print('Авторов в выбоке:', len(authors_df))

Авторов в выбоке: 2011

Конечно, для такого количества авторов столбчатую диаграмму не построишь. Чисто для целей визуализации предлагаю посмотреть тех авторов, которые построили более 20 объектов. Меняем последние строчки кода на:

most_frequent_authors_df = authors_df.loc[(authors_df['count'] > 20)]
    
seaborn.barplot(x=most_frequent_authors_df['count'], 
                y=most_frequent_authors_df['name'])

Чем дальше, тем интересней! И кто же этот Г.А. Симонов, который судя по статистике отстроил чуть ли не половину города, а в честь него даже ни одной улицы не назвали? Заглянем в Википедию:

Григорий Александрович Симонов (23 января 1893, Ташкент — 31 января 1974, Москва) — советский архитектор, инженер, педагог.
Г. А. Симонов родился в Ташкенте. Детство провел в Троицке, там окончил гимназию.
Окончил Петроградский Институт гражданских инженеров в 1920 году, где преподавал с 1929 года.
В 1919—1922 гг. обучался в Академии Художеств.
С 1924 года руководил Проектным Бюро Стройкома.
С 1943 года — заместитель председателя Государственного Комитета по делам архитектуры при Совете Народных Комиссаров СССР.
С 1947 по 1949 год — председатель Комитета по делам архитектуры при Совмине СССР.
С 1955 года — преподаватель Московского архитектурно-строительного института.

Теперь понятно, это послереволюционный, советский архитектор, руководивший сооружением множества объектов в тогдашнем Ленинграде. Если отставить личность товарища Симонова в покое, что еще мы здесь видим? Имеет место выброс статистических данных — Симонов упоминается так часто, что затмевает других. К тому же друг попросил старинные памятники архитектуры. Так что, ни в коем случае не умаляя творцов последнего века, давайте исключим из выборки все сооружения старше 1900 года. Правда тут есть небольшая проблема...

В нашем датасете столбец с датами — текстовый, произвольной длинны и содержания. Где-то стоит одна дата, где-то период, где-то просто двухзначная цифра века постройки. В данном случае мы не можем сделать числовую выборку, а можем работать только со строками. Если бы точность была важна, можно было бы преобразовать этот столбец, как мы сделали с авторами, можно было бы отбирать строки регулярными выражениями. Но это всего лишь поездка друга, так что решаем не усложнять, и волевым решением отсекаем все даты, где встречается число 19. Мы упустим несколько построек, где в поле дата указано "конец 19 в.", "1819" и тому подобных. Примем это как допустимые потери. Снова, мы не сможем визуализировать всех авторов выборки, но отобразим хотя бы наиболее часто встречающихся.

districts_authors_df = pandas.DataFrame(districts_unique, columns=['district'])

for col in most_frequent_authors_df['name']:
    districts_authors_df[col] = 0
    
districts_authors_df.set_index('district', inplace=True)

top_poi = []

for district in districts_unique:
    district_df = data.loc[(data['district'] == district) & 
        (data['protection_category'] == 'Объект культурного наследия федерального значения') & 
        (~data["date"].str.contains('19', na=True))]
    for index, row in district_df.iterrows():
        for author in most_frequent_authors_df['name']:
            if row['author'] == row['author'] and author in row['author']:
                districts_authors_df[author][district] += 1
        for author in authors_df['name']:
            if row['author'] == row['author'] and author in row['author']:
                top_poi.append(row['number'])

seaborn.heatmap(districts_authors_df, xticklabels=True, yticklabels=True) #, annot=True

Ситуация проясняется. У нас теперь есть список наиболее интересных для нашей задачи районов и авторов. Переменная top_poi (Top Points of Interests) содержит номера наших призеров. Но где именно находятся объекты? Пришло время для геокодирования...

Геокодирование

У нас есть колонка с адресами объектов, но если мы хотим определить, как далеко или близко что-то находится, нам нужны координаты. Процесс конвертации адресов в координаты называется геокодирование или геокодинг. У Яндекса есть отличный сервис, который все сделает за нас.

Для начала сделаем копию нашего датасета и добавим туда два новых поля: долгота и широта:

df = pandas.read_csv('spb_memo.csv')
df['lat'] = float('nan')
df['lon'] = float('nan')
df.to_csv("spb_memo_geo.csv", index=False)

У геокодера Яндекса бесплатно только 1000 запросов в сутки. В нашей переменной top_poi содержится чуть меньше 500 объектов, так что на это исследование у нас сегодня хватает:

df = pandas.read_csv('spb_memo_geo.csv')
from decimal import Decimal
import os
from dotenv import load_dotenv
from yandex_geocoder import Client

load_dotenv('.env')
yandex_geo_api_key = os.environ.get("yaGeoApi")
client = Client(yandex_geo_api_key)

coordinates = 0
api_limit_per_day = 1000

for poi in top_poi:
    if api_limit_per_day > 0:
        poi_row = df.loc[(df['number'] == poi)]
        if poi_row.empty:
            lat = float('nan')
        else:
            lat = list(poi_row['lat'])[0]
        addr = list(poi_row['address'])
        if lat != lat and len(addr) > 0 and addr[0] == addr[0] and len(addr[0]) > 5:
            coords = client.coordinates("Санкт-Петербург, " + addr[0])
            df.loc[poi_row.index, 'lon'], df.loc[poi_row.index, 'lat'] = coords
            api_limit_per_day -= 1

df.to_csv("spb_memo_geo.csv", index=False)

Последний отсев

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

df = pandas.read_csv('spb_memo_geo.csv')
lat_lon = df[df['number'].isin(top_poi)]
x = list(lat_lon['lat'])
y = list(lat_lon['lon'])

seaborn.scatterplot(x=x, y=y, c=['red'])

Вы будете правы, если скажите, что такое расположение точек не соответствует сторонам света — по X должна быть долгота, а по Y — широта. Но на данном этапе нам это и не нужно. Мы ищем зависимости, отклонения, совпадения и т.д.

И что это за одинокая точка в левом верхнем углу, из-за которой у нас верхняя половина графика пустая?

print(df.loc[(df['lat'] < 59.8) & (df['lon'] > 31)])

7922    Летний сад
Name: name, dtype: object

Ах, это Летний сад, который друг упоминал в своей просьбе. Туда он пойдет независимо от расстояния. Вычеркиваем его из выборки, к тому же для нас это небольшой выброс данных, делающий картину менее понятной.

lat_lon = df.loc[df['number'].isin(top_poi) & (df['lat'] > 59.6) & 
                 (df['lat'] < 62) & (df['lon'] > 29) & (df['lon'] < 30.8)]
x = list(lat_lon['lat'])
y = list(lat_lon['lon'])

seaborn.scatterplot(x=x, y=y, c=['red'])

Итог

Итак, пришло время ответить на самый главный вопрос нашего исследования. Раз друг хочет поселиться где-то в центре, чтобы независимо от того какой объект он выберет, это было относительно недалеко, просто возьмем среднее значение всех координат.

xy_center = (sum(x) / len(x), sum(y) / len(x)) 
seaborn.scatterplot(x=y, y=x, c=['red'])
seaborn.scatterplot(x=[xy_center[1]], y=[xy_center[0]], c=['green'])

Зеленая точка и есть наше заветное место. Интересно, а где это вообще?

print(xy_center[1], ',', xy_center[0])

59.926768, 30.294057

Снова прибегнем к геокодеру Яндекса, но на этот раз в обратном направлении — для преобразования координат обратно в адрес:

address = client.address(Decimal("30.294057"), Decimal("59.926768"))
print(address)

Россия, Санкт-Петербург, набережная Крюкова канала

Вот оно! Когда все только началось, лично мне казалось, что центральная точка окажется где-то на Дворцовой площади или в Петропавловской крепости, но нет, другие архитектурные объекты оттянули ее на себя, и мы получили набережную Крюкова канала.

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

Финальная визуализация

Это было чисто статистическое исследование, — максимум, что удалось выжать из этого датасета. Исследование не имело отношение к истории. Мы не узнали, какие архитектурные стили преобладают в Санкт-Петербурге, и как они менялись со временем; не увидели в выборках архитекторов-основателей. Надеюсь обо всем этом друг узнал из посещенных им экскурсий.

Напоследок, давайте сделаем все красиво, и действительно сопоставим наши данные с реальной картой, используя библиотеку OSMnx, основанную на данных OpenStreetMap (спасибо Carlos Lannister за подсказку):

import osmnx as ox

# Center of map
latitude = 59.939099
longitude = 30.315877

point = (latitude, longitude)
G = ox.graph_from_point(point, dist=10000, retain_all=True, simplify = True, 
                        network_type='all')

u = []
v = []
key = []
data = []
for uu, vv, kkey, ddata in G.edges(keys=True, data=True):
    u.append(uu)
    v.append(vv)
    key.append(kkey)
    data.append(ddata)    

# List to store colors
roadColors = []
roadWidths = []

for item in data:
    if "length" in item.keys():
        if item["length"] <= 100:
            linewidth = 0.10
            color = "#a6a6a6" 
            
        elif item["length"] > 100 and item["length"] <= 200:
            linewidth = 0.15
            color = "#676767"
            
        elif item["length"] > 200 and item["length"] <= 400:
            linewidth = 0.25
            color = "#454545"
            
        elif item["length"] > 400 and item["length"] <= 800:
            color = "#d5d5d5"
            linewidth = 0.35
        else:
            color = "#ededed"
            linewidth = 0.45
    else:
        color = "#a6a6a6"
        linewidth = 0.10
            
    roadColors.append(color)
    roadWidths.append(linewidth)
            
bgcolor = "#061529"

fig, ax = ox.plot_graph(G, node_size=0,figsize=(27, 40), 
                        dpi = 300,bgcolor = bgcolor,
                        save = False, edge_color=roadColors,
                        edge_linewidth=roadWidths, edge_alpha=1,
                        show=False, close=False)

for i in range(len(x)): #
    ax.scatter(y[i], x[i], s = 100, c='red')

ax.scatter(xy_center[1], xy_center[0], s = 100, c='green')

Как видно из кода, красные точки — это объекты культурного наследия Санк-Петербурга (в данном случае порядка 1 300 объектов, без ограничения по годам), а зеленая точка — рекомендованное место старта.

Где же мистика?

Собрав все правки, пожелания и рекомендации из комментариев я еще раз переделал расчёты:

• Убрал из выборки Петродворцовый район

• Добавил Приморский район

• Ограничил выборку только объектами имеющими федеральное значение

• Убедился, что значимые архитекторы города, такие как Бенуа, Растрелли, Трезини и др. включены в выборку

В результате всех этих коррекций точка сместилась на адрес Россия, Санкт-Петербург, Инженерная улица, 2-4А, тот самый мистический Инженерный замок, с его призраком, легендами и загадочной надписью на фасаде. Ну чем не мистика?

Ниже обновленная финальная визуализация со всеми коррекциями и новой центральной точкой:

Если захотите проверить мои расчёты, можете найти все файлы в моём репозитории Github.