Разметка и анализ геоданных с использованием QGIS + Python

Доброго времени суток!

На связи Мурачёв Виктор и Князев Платон.

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

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

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

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

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

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

Разметка и анализ геоданных с использованием QGIS + Python

Территория России находится в четырех климатических поясах: арктический, субарктический, умеренный и тропический. Можно увидеть это разделение на карте. Но для детального анализа этих данных мало изображения в формате .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. Спасибо за внимание.