Время прочтения: 10 мин.

Геоданные нужны в основном для создания карт. Зачастую разработчики используют их для создания различных сервисов, отслеживания маршрутов, прогнозирования поведения потребителей и решения прочих задач на рассматриваемых территориях.

В том числе, использование такого типа данных помогает в финансовом планировании для предприятий, например, сельскохозяйственного или страхового сектора. Такие предприятия могут оценить риски, сопутствующие с владением землей на той или иной территории, подготовиться к возможным катаклизмам. Географические данные могут помочь в создании инструментов для предсказания будущих тенденций, таких как рост или снижение количества жителей в конкретном регионе, оценка будущих экологических трендов и других важных факторов.

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

В данной публикации нами будет рассмотрено использование бесплатного приложения QGIS для получения координат объектов из карт различного формата. Мы также проведем исследование и полученных данных с использованием модуля GeoPandas для Python, помимо этого, данные будут представлены графически с использование модуля Plotly.

В качестве примера попробуем рассмотреть вхождение районов и областей России в различные климатические зоны, опоясывающие нашу страну:

Территория России находится в четырех климатических поясах: арктический, субарктический, умеренный и тропический. Можно увидеть это разделение на карте. Но для детального анализа этих данных мало изображения в формате .jpg. С помощью приложения QGIS можно наложить изображение на любой участок и, разметив его, получить координаты интересующих объектов на карте.

Первым шагом после несложной установки приложения было создание и выбор формата шаблона карты. В рамках данной задачи будет исследована территорию России. В нашей стране, согласно государственным стандартам, в качестве системы координат допускается использование референсной СК Pulkovo-1942.

Для начала, в приложении QGIS был создан новый проект, где в качестве источника был использован слой карты OpenStreetMap Standard из модуля QuickMapServices:

В свойствах слоя была изменена стандартная система координат на СК, в которой отрисована исследуемая карта (Pulkovo 1942):

Как можно видеть на рисунке снизу, выбор такой системы координат оказался верным. Формат карты очень похож на формат источника:

После выбора СК проекта, необходимо было совместить изображение-источник со слоем карты. Это можно сделать, выбрав пункт «Слои» — «Привязка растров» на панели управления.  Откроется меню привязки. После этого можно загрузить растровое изображение исследуемого участка:

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

Чтобы правильно связать объекты, необходимо несколько мест привязки (имеющаяся карта была привязана по 21 позиции). Советую сохранить координаты точек на случай ошибок в работе приложения. Можно запускать процесс привязки:

После успешной привязки, появился новый слой с исходным изображением. Сделав его видимым, мы удостоверились, в корректности совмещения источника и карты:

Если все привязалось правильно, можно переходить к разметке. В ходе решения этой задачи, климатические зоны из источника были размечены вручную. Для этого был использован слой формата ShapeFile, с типом геометрии «Полигон»:

После создания слоя, был активирован режим редактирования и добавлен полигональный объект:

Далее на слой карты были нанесены полигоны, повторяющие климатические зоны из источника. Например, эта часть острова «Северный» находится в арктическом климатическом поясе:

Таким способом были размечены все климатические зоны из источника:

Чтобы продолжить работу с размеченными полигонами, они были сохранены в файл формата .shp – для этого нужно нажать правой кнопкой на слой, выбрать «экспорт» — «сохранить объекты как…». И в открывшемся меню задать путь к файлу и используемую систему координат:

Теперь можно «прочитать» экспортированные полигоны. Эта задача решалась в интерфейсе Jupyter Notebook. Код проекта и все вложения есть в репозитории на GitHub: https://github.com/slowpokalipsis/nta_geo

Для работы с файлами были использованы библиотеки Pandas и GeoPandas для Python:

import geopandas as gpd
import pandas as pd

# чтение файлов с полигонами зон
arctic = gpd.read_file('qgis_data/arctic_zone.shp') 
subarctic = gpd.read_file('qgis_data/subarctic_zone.shp') 
umerenny = gpd.read_file('qgis_data/umerenny_zone.shp')
subtropic = gpd.read_file('qgis_data/subtropic_zone.shp')

# объединение разных фреймов в один, дополнение цветом и названиями
zones_list = [arctic, subarctic, umerenny, subtropic]
zone_names = ['Арктический', 'Субарктический', 'Умеренный', 'Субтропический']
zone_colors = ['mediumpurple', 'powderblue', 'yellowgreen', 'red']

climate_zones = gpd.GeoDataFrame()
for zone_ind in range(len(zones_list)):
    zone = zones_list[zone_ind][['geometry']]
    zone['level'] = zone_names[zone_ind]
    zone['zone_color'] = zone_colors[zone_ind]
    climate_zones = climate_zones.append(zone)
climate_zones.reset_index(inplace=True)
climate_zones.rename(columns={'index': 'id'}, inplace=True)
climate_zones = climate_zones.to_crs('ESRI:102027') # удобный формат отображения

Теперь полученные полигоны можно отобразить. Для этого были применены возможности библиотеки Plotly:

import plotly.graph_objects as go
from plotly.graph_objs import Layout

# отрисовка Plotly
def plot_climate_zones(climate):
    
    layout = Layout(
        paper_bgcolor='rgba(0,0,0,0)',
        plot_bgcolor='rgba(0,0,0,0)'
    )
    
    fig = go.Figure(layout=layout)
    
    for i, vals in tqdm(climate.iterrows()):
        x, y = np.array(vals.geometry.exterior.coords.xy)
        dist = vals['level']
        zone = vals['zone_color']
        zone_color = zone.translate(remove_digits)
        fig.add_trace(go.Scatter(x=x, y=y, text=dist, line_color=zone_color,
                                fill='toself',fillcolor=zone_color, line_width=0.3, 
                                 mode='lines', opacity=0.3))
        
    fig.update_layout(showlegend=False, dragmode='pan', width=1000, height=1000)
    fig.update_xaxes(visible=False)
    fig.update_yaxes(visible=False, scaleanchor="x", scaleratio=1)
    
    return fig

# создание объекта карты
mymap = plot_climate_zones(climate_zones)

Для отображения полученной карты используем следующую команду:

config = dict({'scrollZoom': True})
mymap.show(config=config)

Имеем картинку с полигонами климатических зон в границах РФ:

Но самих границ пока нет. В свободном доступе был найден файл с координатами границ и областей РФ. Добавим его:

# карта России geojson
ru_map = gpd.read_file('map_data/russia_regions_stats.geojson')
ru_map = ru_map.to_crs('ESRI:102027')

# объекты в файле хранятся в формате «мультиполигонов»
# для отрисовки необходимо преобразовать их в наборы координат X, Y
def geoseries2shapedict(gdf, name_col, gdf_type):
    regions = {}
    for und, reg in tqdm(gdf.iterrows()):
        if type(reg.geometry) == MultiPolygon:
            x, y = np.array([[], []])
            for poly in reg.geometry.geoms:
                x_, y_ = poly.exterior.coords.xy
                x, y = (np.append(x, x_), np.append(y, y_))
                x, y = (np.append(x, None), np.append(y, None))
            x = x[:len(x)-1]
            y = y[:len(y)-1]
        else:
            x, y = np.array(reg.geometry.exterior.coords.xy)
        
        regions[reg[name_col]] = x, y, reg.federal_district
    return regions

В функцию для отрисовки, использующей модуль Plotly, добавим код для отображения регионов на карте:

for reg, vals in tqdm(ru_regions.items()):
        x, y, dist = vals
        fig.add_trace(go.Scatter(x=x, y=y, text=dist, 
                                    line_color='grey',
fill='toself',
fillcolor='rgba(255,255,255,0)', line_width=1))

Применим новые функции:

# преобразование MultiPolygon в координаты X, Y 
ru_map['area'] = ru_map.geometry.apply(lambda x: x.area)
ru_map_regions = geoseries2shapedict(ru_map.sort_values(by=['area'], 
                                                        axis=0, 
                                                        ascending=False), 
                                     'name_ru', 'rumap')
# создание объекта карты
mymap = plot_ru_zones(ru_map_regions, climate_zones)

Получилась карта, на которой обозначены регионы и климатические зоны:

В рамках ранее упомянутой задачи по поиску сельскохозяйственных участков в регионах с повышенной пожароопасностью, из разных источников был агрегирован набор данных с координатами и кодами всех районов в регионах России. Данные для его создания и скрипты разработки находятся в репозитории. В рамках этой публикации не будет рассмотрен процесс создания и сопоставления данных из разных источников. Используем готовый набор данных:

# чтение и преобразование данных
all_ru_regions = pd.read_csv('map_data/ru_regions_kadastr.csv')
all_ru_regions = gpd.GeoDataFrame(all_ru_regions, geometry=wkt.loads(all_ru_regions.geometry),
     crs='ESRI:102027')

Полученный набор имеет поля, содержащие название, геометрию (координаты вершин полигона), кадастровый код региона, в котором находится район, код самого района:

Добавим в функцию отрисовки Plotly:

for reg, vals in tqdm(regions.items()):
        x, y, dist = vals
        fig.add_trace(go.Scatter(x=x, y=y, text=reg, line_color='black', line_width=0.3,
                                 mode='lines'))

mymap = plot_reg_zone(climate_zones, all_ru_regions)

При наведении на карту будет отображаться название выбранного района. Функцию отображения при наведении можно настроить множеством различных способов:

Теперь можно определить, какие районы входят в ту или иную климатическую зону, используя функцию geopandas.geometry.contains():

all_ru_regions[['zone', 'zone_id']] = [None, None] # новые столбцы для указания зоны

# уменьшение площади района для более точного определения его вхождения в зону
all_ru_regions.geometry = all_ru_regions.geometry.scale(xfact=0.1, yfact=0.1, origin='center')

# функция определения вхождения района России в климатическую зону
for i, all_val in all_ru_regions.iterrows():
    for j, f_val in climate_zones.iterrows():
        try:
            if f_val.geometry.contains(all_val.geometry):
                all_ru_regions.iloc[i,-2:] = [f_val.level, f_val.id]
                list_j.append(j)
                cnter += 1
            else:
                pass
        except:
            pass

# после обработки увеличиваем площадь района до исходных значений
all_ru_regions.geometry = all_ru_regions.geometry.scale(xfact=10, yfact=10, origin='center')

После анализа данных была выявлена следующая статистика нахождения регионов в климатических зонах:

print(all_ru_regions.zone.value_counts(normalize=True))

Абсолютное большинство районов РФ находятся в умеренной климатической зоне. Это 97.95 % районов, что составляет 72.55% от общей площади страны. Перейдем к другому, более показательному примеру.

Возьмем данные об очагах пожаров. Датасет был найден в открытом доступе на сайте DSWorks:

# чтение и преобразование данных
fire_df = pd.read_csv('map_data/train_raw.csv')
fire_df = gpd.GeoDataFrame(fire_df, geometry=gpd.points_from_xy(fire_df.lon, fire_df.lat),
                           crs='EPSG:4326')
fire_df = fire_df.to_crs('ESRI:102027')

# фильтрация данных, удаление данных о контролируемых пожарах
fire_df = fire_df[fire_df.type_name != 'Контролируемый пал']

В наборе данных содержатся записи об очагах пожаров с 2012 по 2021 год с координатами и типом пожара:

Совместим эти данные с картой России. Для этого изменим функцию отображения:

# heat_map preparation
# преобразование координат пожара к единому формату X, Y    
    fires['x'] = fires.geometry.x
    fires['y'] = fires.geometry.y

    ny = 100 # количество пикселей (шагов) по оси X (100 чтобы не перегружать карту)
    nx = int(ny*(fires.x.max() - fires.x.min()) / (fires.y.max() - fires.y.min()))
    
# отрисовка тепловой карты
fig = go.Figure(go.Histogram2d(x=fires.x, y=fires.y,
                                nbinsx=nx,
                                nbinsy=ny,
                                colorscale='hot',
                                showlegend=False, showscale=False, reversescale=True
                            ), layout=layout)

Получилось следующее распределение. На карте цветными точками обозначены скопления очагов пожаров – чем больше пожаров, тем темнее:

Попробуем определить самые пожароопасные районы России по соотношению количество пожаров и площади:

# добавление столбца для подсчета пожаров
all_ru_regions['fire_count'] = None

# функция подсчета пожаров 
for i, val in tqdm(all_ru_regions.iterrows(), total = all_ru_regions.shape[0]):
    all_ru_regions.iloc[i, -1] = np.sum(all_ru_regions.iloc[i,3].contains(fire_df.geometry))

# вычисление отношения количества очагов пожаров к территории района
all_ru_regions['area_frac'] = all_ru_regions['area'] / 1e5 
all_ru_regions['fires_area_rate'] = all_ru_regions['fire_count']/all_ru_regions['area_frac']
# создание переменной знаменателя для нормализации отношения 
rate_multiply_value = 1 / all_ru_regions['fires_area_rate'].max()
# нормализация отношения
all_ru_regions['fires_area_rate'] = all_ru_regions['fires_area_rate'] * rate_multiply_value
all_ru_regions['fires_area_rate'] = all_ru_regions['fires_area_rate'].apply(lambda x: round(x, 3))
all_ru_regions.drop(columns=['area_frac'], inplace=True)

В результате анализа данных были выявлены наиболее пожароопасные регионы. В их числе Находкинский, Дальнереченский, Спасск-Дальний, Большой Камень и прочие районы на Дальнем Востоке России:

В результате были решены задачи разметки и анализа данных с помощью QGIS и модуля Geopandas.  У обоих утилит, конечно же, есть и другие применения. Например, используя Plotly и GeoPandas, можно находить пересечения объектов, создавать более интерактивные карты с применением, рисовать 3D проекции и прочее.

В ходе решения задач были использованы ресурсы:

https://gis-lab.info/

https://geopandas.org/en/stable/

Код проекта размещен в репозитории github. Спасибо за внимание.