Время прочтения: 5 мин.
В рамках одного проекта мне потребовалось выяснить, какие объекты нашей организации находятся в зоне вечной мерзлоты.
Выявить эти объекты удалось с помощью библиотеки geopandas, которая прекрасно, хоть и не всегда интуитивно, работает с географическими объектами, а также используя API Яндекс.
Работал с данным в Google Colab, так как Anaconda для Windows потребует актуализации слишком большого количества зависимостей, и не факт, что это получится сделать корректно. В очередной раз убеждаюсь, что Python рожден для работы в Linux. В другом похожем проекте у меня получилось работать с этой библиотекой, установив ее в Windows Linux Subsystem, и подключившись к ядру Jupyter Lab из браузера, запущенного в Windows.
Для начала подключим свой диск в Google Drive. На нем будем сохранять полученные данные, так как при длительном неиспользовании Colab сбрасывает вашу сессию и сохраненные внутри сессии данные теряются.
from google.colab import drive
drive.mount('/content/drive')
Теперь установим необходимые библиотеки:
!pip install geopandas
!pip install pygeos
Импортируем используемые в процессе исследования библиотеки:
import numpy as np
import pandas as pd
#import shapefile as shp
import matplotlib.pyplot as plt
import seaborn as sns
import geopandas as gpd
sns.set(style="whitegrid", palette="pastel", color_codes=True)
sns.mpl.rc("figure", figsize=(20,20))
from shapely import wkt
%matplotlib inline
Теперь загрузим необходимые датасеты.
Нам потребуются shp-файлы, которые geopandas умеет преобразовывать в нужный ему формат. К сожалению, не все shp-файлы корректны, и в некоторых случаях geopandas выбрасывает ошибку импорта.
Датасет geopandas представляет собой обычный датасет pandas, который содержит дополнительные колонки, содержащие информацию о географических точках и полигонах.
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
permice_nasa = gpd.read_file("/content/drive/My Drive/Colab Notebooks/eco/maps/permaice.shp")
Датасет permice_nasa скачан с сайта Nasa.
Он содержит информацию в виде полигонов о типах вечной мерзлоты в северном полушарии Земли. К сожалению, Роскосмос и российские университеты такими данными не делятся, так что приходится использовать иностранные источники.
Давайте нарисуем полученные данные:
# Выделяем отдельным датасетом Европу и Россию
europe = world[world.continent == 'Europe']
russia = europe[europe.name =='Russia']
#приводим наши датасеты к одной проекции
russia = russia.to_crs(epsg=3035)
permice_nasa = permice_nasa.to_crs(epsg = 3035)
#совмещаем датасеты на одном изображении
base = permice_nasa.plot(figsize=(50,50))
russia.plot(ax = base, color = 'red', alpha = 0.1)
Обратите внимание на конструкцию gpd.to_crs(epsg = 3035). Дело в том, что существуют разные способы проецирования круглого земного шара на плоский лист бумаги (или монитор компьютера), и для того, чтобы успешно работать с разными датасетами в рамках одного проекта, нужно привести их к одному виду. В данном случае мы привели их к проекции под номером 3035, с видом на земной шар со стороны полюса.
По большому счету, визуализация в данном проекте нам не нужна, так как мы работаем не с картинками, а с данными. Но если вам нужно красиво нарисовать какие-то данные на географической карте, geopandas сделает это для вас (возможно, не без помощи PhotoShop).
Теперь нам нужно получить координаты объектов, которые содержатся в отдельном датасете. Карта вечной мерзлоты не отличается абсолютной точностью, поэтому достаточно понять, в какой зоне находится населенный пункт.
Для этого мы воспользуемся сервисом геокодирования от Yandex — tech.yandex.ru/maps/geocoder. Он позволяет отрабатывать бесплатно до 25 тысяч запросов в сутки. Придется проигнорировать требование Yandex о размещении полученных данных на картах. Впрочем, мы же используем эти данные не в коммерческих целях, верно? Сначала получите ключ для доступа к API. Данные для поиска я заранее привел в формат типа “р-н Троицкий+с Бобровка”, что позволило снизить количество ошибок. Но все равно, некоторые населенные пункты из российской глубинки могут внезапно очутиться, по мнению Яндекса, в Карибском бассейне или на Африканском континенте.
Подождав некоторое время, мы получаем еще один датасет, теперь уже содержащий координаты интересующих нас объектов:
import requests
#функция получения координат по адресу с Яндекса
def yandex_geocode(adress):
url = 'https://geocode-maps.yandex.ru/1.x'
params = {'geocode': adress,
'apikey': 'здесь вставьте ваш api-ключ', 'format': 'json'}
r = requests.get(url, params=params)
results = r.json()
#координаты
try:
pos = results['response']['GeoObjectCollection']['featureMember'][0]['GeoObject']['Point']['pos']
except:
pos = 'error'
#адрес для проверки
try:
ya_adress = results['response']['GeoObjectCollection']['featureMember'][0]['GeoObject']['metaDataProperty']['GeocoderMetaData']['text']
except:
ya_adress = 'error'
return pd.Series([pos, ya_adress])
def convertAdressToGPD(cities):
#получение координат из Яндекса
#первое - долгота Longitude, второе - широта Latitude
cities[['coorditate', 'ya_adress']] = cities.apply(lambda x: yandex_geocode(x['Адрес']), axis = 1)
cities[['Longitude', 'Latitude']] = cities.coorditate.str.split(" ", expand = True)
cities_correct = cities[cities['coorditate']!='error']
gpd_cities = gpd.GeoDataFrame(
cities_correct, geometry=gpd.points_from_xy(cities_correct.Longitude, cities_correct.Latitude), crs={'init': 'epsg:4326'})
gpd_cities = gpd_cities.to_crs(epsg = 3035)
return gpd_cities
#читаем список населенных пунктов из файла Excel
#Адрес сохраняем с столбце «Адрес»
cities = pd.read_excel('/content/drive/My Drive/Colab Notebooks/eco/cities/7000.xlsx')
gpd_cities = convertAdressToGPD(cities)
Осталось немного — определить, в какой зоне находится нужный нам объект (а вечная мерзлота, как выяснилось, имеет весьма широкую градацию).
Для этого пишем еще одну функцию, которая определяет, для каких из полигонов датасета permise_nasa объект city имеет значение True:
def frost_zone(city):
try:
zone_df = permice_nasa.geometry.contains(city.geometry)
return permice_nasa.loc[zone_df[zone_df == True].index].COMBO.iloc[0]
except:
return None
gdf_cities['frost_zone'] = gdf_cities.apply(lambda x: frost_zone(x), axis = 1)
Теперь наш датафрейм gpd_cities содержит колонку, в которой указан тип полигона из датасета permice_nasa
Выгружаем его в Excel и изучаем полученный результат.
gpd_cities.to_excel('/content/drive/My Drive/Colab Notebooks/eco/cities/cities_frost.xlsx')
Модуль geopandas имеет весьма широкие возможности по анализу и визуализации географических данных, но, к сожалению, в рускоязычном сегменте Интернета хороших материалов по нему не так много. Для изучения других его возможностей рекомендую почитать стандартную документацию к библиотеке и книгу Э.Вейстра “Разработка приложений на языке Python”.