Python, Анализ данных

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

Время прочтения: 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”.

Советуем почитать