Понадобилось мне недавно нарисовать в 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.
Комментарии (5)
gematit
14.03.2022 06:42Очень удобный и простой способ.
Тем не менее, для Python могут быть более удобны пакеты Basemap/Cartopy, которые лишь немного сложнее в рамках описанной задачи. Но зато там не придётся пересчитывать координаты.
ScapeKP Автор
14.03.2022 09:46Спасибо, посмотрю. GeoPandas выглядит привычнее, потому что много работаю с Pandas, но может и правда лучше отказаться от него...
bozman
Все-таки карту надо рисовать в картографической проекции, для России это Albers Siberia. Сейчас у вас нечто, отдаленно напоминающее растянутую Equirectangular, но в ней даже в нерастянутой земные карты рисовать не принято.
Geopandas позволяет задавать проекцию для карты, и это надо делать обязательно. Ну, как минимум, чтобы контуры хоть узнавались. И при задании проекции у вас отпадет необходимость двигать кусок — он сам встанет куда надо.
ScapeKP Автор
Я читала эту статью и пробовала разные проекции, прежде чем начала переносить вручную, но не смогла найти подходящую. Можете подсказать, какой код у Albers Siberia? Не могу ее найти. На сайте https://epsg.io ее нет, просто гугл тоже не помогает. Кажется, что задавать ее надо вручную (а очень не хочется для чернового анализа так заморачиваться). Взяла наиболее близкую - Europe Albers Equal Area Conic. С ней пропорции лучше, но без дополнительных обработок выглядит вот так :(
bozman
Извините, что сразу не ответил — мне можно раз в день, рассердил чем-то разум улья. Вам и нужна проекция Альберса, повернутая на 105 градусов и центрированная по координатам 56,100. Для проходной работы достаточно задать поворот и все будет норм. Артефакт непонятный, тут не знаю что делать. Может помочь задание границ карты (по широтам/долготам), но я не знаю, как это делается конкретно в геопандасе.