Geo data in Python

Понадобилось мне недавно нарисовать в Python данные на карте, благо в данных есть координаты. Казалось бы, что может быть сложного... Но обо всем по порядку.

Мои точки - это города. Города имеют разные признаки, как качественные (используем разные цвета точек на карте), так и количественные (используем разный размер кружочка).

Первым делом нужно установить и импортнуть библиотеку geopandas - она показалась мне самой простой и понятной среди всех вариантов работы с геоданными.

import geopandas as gpd

Нужно нарисовать контуры России. Для этого нужен файл с границами - geojson. Его можно очень быстро нагуглить (я взяла отсюда). Мне нужны лишь контуры страны, поэтому я взяла файл 2 уровня - admin_level_2.geojson

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

map_df=gpd.read_file('my_env/admin_level_2.geojson') 
map_df=gpd.GeoDataFrame(map_df) 
map_df.plot(figsize=(14,10), color='black', alpha=0.1)

3 простых шага, которых достаточно достаточно почти в любой стране мире. Но есть одно "но" - Россия реально большая страна.

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

Поэтому сначала нужно сделать преобразование координат, прибавив к оторвавшемуся куску 360^о.

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

Но нет. Тут мой код зависал всерьез и надолго. Геоданные - огромный тяжелый массив и перелопачивать их все jh явно не желал. Недолго думая я приняла решение пожертвовать Кемской волостью мелкими полигончиками.

Я отрисовала все полигоны по отдельности и выяснила, что основная часть - это полигон 28. Отрезанный кусок дальнего востока - 26. Сдвинуть один на 360^о и отрисовать только эти два полигона занимает буквально секунду. Итоговый код занимает 5 строчек, а результат вполне меня устраивает.

map_df=gpd.read_file('my_env/admin_level_2.geojson') 
map_df=gpd.GeoDataFrame(map_df)
p = gpd.GeoSeries(shapely.affinity.translate(map_df['geometry'][0][26], xoff=360, yoff=0))
p=p.append(gpd.GeoSeries(map_df['geometry'][0][28]))
p.plot(figsize=(14,10), color='black', alpha=0.1)

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

select city
		,lat
    ,case when lon<0 then lon+360 else lon end as lon
    ,category
    ,qty
from my_table

Преобразуем поля широты и долготы в одно поле "координаты" и преобразуем DataFrame в GeoDataFrame.

df['coordinates'] = df[['lon', 'lat']].values.tolist()
df['coordinates'] = df['coordinates'].apply(Point)
df = gpd.GeoDataFrame(df, geometry='coordinates')

Всё. Осталось только нарисовать оба набора данных вместе. Важно помещать обе строки в одну ячейку jupyterhub'а, иначе не сработает.

ax_all=p.plot(figsize=(14,10), color='black', alpha=0.1)
df.plot(ax=ax_all, column = 'category',markersize=20,marker='o')

Или пример с размером кружочков (на урезанном наборе данных для наглядности).

ax_all=p.plot(figsize=(14,10), color='black', alpha=0.1)
df_size.plot(cmap='autumn',ax=ax_all, markersize='qty',marker='o',alpha=0.6)

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