Определяем принадлежность объектов к зоне вечной мерзлоты с помощью GeoPandas

Как GeoPandas поможет определить принадлежность объектов к зоне вечной мерзлоты.

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

Выявить эти объекты удалось с помощью библиотеки 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

По большому счету, визуализация в данном проекте нам не нужна, так как мы работаем не с картинками, а с данными. Но если вам нужно красиво нарисовать какие-то данные на географической карте, 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».