Обзор готовых скриптов Python для работы с геоданными

Привет всем! Я Татьяна, работаю ГИС-инженером в ИнноГеоТехе, хочу поделиться одним инструментом, который помогает оптимизировать процессы в ГИС-проектах — и это Python.

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

В ГИС-проектах мы часто используем Python, рутинные задачи решаются гораздо быстрее — с помощью скриптов мы проводим постобработку космических снимков.

Перейдем к примерам общих скриптов, картографы могут использовать их для оптимизации своей работы.

1 скрипт — обработка и анализ геоданных, при помощи функции head() мы можем посмотреть первые несколько строк данных.

Обработка и анализ геоданных import geopandas as gpd # Загрузка данных data = gpd.read_file('data.geojson') # Анализ данных print(data.head()) # Выполнение пространственного анализа data['area'] = data.geometry.area data.to_file('output_data.geojson', driver='GeoJSON')

2 скрипт — преобразование систем координат, при помощи функции to_crs() мы можем перевести данные из исходной системы координат in_proj в целевую систему координат out_proj.

Преобразование систем координат import geopandas as gpd from pyproj import Proj, transform # Загрузка данных data = gpd.read_file('data.geojson') # Преобразование координат in_proj = Proj(init='epsg:4326') out_proj = Proj(init='epsg:3857') data['geometry'] = data['geometry'].to_crs(out_proj) # Сохранение данных data.to_file('output_data.geojson', driver='GeoJSON')

3 скрипт — при помощи данной функции Python мы можем визуализировать данные на карте — folium.Map — создаем карту, folium.Geojson — добавляем данные на карту, save — сохраняем карту.

Создание карт import folium import geopandas as gpd # Загрузка данных data = gpd.read_file('data.geojson') # Создание карты m = folium.Map(location=[data['geometry'].centroid.y.mean(), data['geometry'].centroid.x.mean()], zoom_start=10) # Добавление геоданных на карту for idx, row in data.iterrows(): folium.GeoJson(row['geometry']).add_to(m) # Сохранение карты m.save('map.html')

База Python, которая понадобится для автоматизации ГИС-задач:

— понимание концепций Python, таких как переменные, типы данных, циклы, условия;

— умение работать с файлами и данными в Python, включая чтение и запись геоданных в различных форматах например, shapefile, GeoJSON;

— навыки работы с библиотеками таких как GDAL, Fiona, Shapely, GeoPandas, Pyproj.

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

Несколько рекомендаций по используемым версиям библиотек Python, которые помогут ГИС-специалистам написать скрипт и избежать ошибок.

1. Используйте Python 3, многие библиотеки уже перешли на эту версию и большинство новых проектов также используют Python 3.

2. Библиотека GDAL — рекомендуется использовать GDAL версии 2.x или выше, так как она поддерживает множество форматов данных и функций.

3. Последние версии понадобятся для библиотек: Shapely (для выполнения геометрических операций), GeoPandas, Rasterio (библиотека для работы с растровыми данными).

Переходим к готовым скриптам — для автоматизации выгрузки растровых слоев лучше проще всего использовать QGIS.

Преимущество — единое рабочее пространство, в котором ГИС-специалисты смогут быстрее проводить пакетную обработку растровых и векторных данных с сохранением в целевые директории.

Для выполнения небольших скриптов в QGIS применяем встроенную консоль Python. Консоль открывается как немодальное окно. Используем специальный интерфейс, можно обращаться к карте, меню, панелям инструментов и другим модулям QGIS.

В качестве исходных данных у нас был заранее открытый проект QGIS с набором растровых слоев. С чего начать ?

— Создаем папку, название папки = названию растрового слоя;

— Создаем новый растровый слой в другой проекции;

— Выполняем сжатие lzw;

— Задаем «Тип данных» — 8 бит, задаем режим сохранения «Изображение»;

— Сохраняем растры в папки, соответствующие названиям растров. Для взаимодействия с QGIS предназначена переменная qgis.utils.iface, которая является экземпляром класса QgisInterface.
На основе данного интерфейса, создаем команды:

x_form = QgsCoordinateTransform(source_crs, dest_crs, QgsProject.instance()) - перепроецирование слоя opts = ["COMPRESS=LZW"] file_writer.set Create Options(opts) - задать сжатие pipe = QgsRasterPipe() pipe.set(provider.clone()) pipe.set(renderer.clone()) - задать режим сохранения "Изображение" file_writer.writeRaster( pipe, raster.width(), raster.height(), x_form.transform(raster.extent()), dest_crs, ) - сохранить растр

Готовый скрипт для работы с выгрузкой растровых файлов в QGIS

import os layers = iface.mapCanvas().layers() path_folder = 'C:/Users/User/Documents/Нурское_СП/Нурское_СП/' opts = ["COMPRESS=LZW"] dest_crs = QgsCoordinateReferenceSystem("EPSG:4326") for raster in layers: source_crs = raster.crs() x_form = QgsCoordinateTransform(source_crs, dest_crs, QgsProject.instance()) path = path_folder + raster.name() if not os.path.exists(path): os.makedirs(path) print(f'Путь {path} создан') else: print(f'Путь {path} уже существует') file_name = path_folder + raster.name() + '/' + raster.name() +'.tif' print(path_folder) file_writer = QgsRasterFileWriter(file_name) file_writer.setCreateOptions(opts) provider = raster.dataProvider() renderer = raster.renderer() pipe = QgsRasterPipe() pipe.set(provider.clone()) pipe.set(renderer.clone()) projector = QgsRasterProjector() projector.setCrs(source_crs, dest_crs, QgsProject.instance().transformContext()) pipe.insert(2, projector) file_writer.writeRaster( pipe, raster.width(), raster.height(), x_form.transform(raster.extent()), dest_crs, ) del file_writer

Еще один вариант готового скрипта — обрезка растровых файлов с помощью библиотек Python — Rasterio и Geopandas.

Используя библиотеку Geopandas, обрабатываем векторный слой AOI, фильтруем геометрии по типам, оставив только полигональные объекты.

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

out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)

и сохраняем в соответствии со значением из атрибутивной таблицы объекта.

Разделили скрипт на этапы:

Шаг 1. Импорт библиотек

# импорт библиотек import sys import glob import geopandas as gpd import rasterio import rasterio.mask

Шаг 2. Указываем путь до исходных файлов

# указываем путь до исходных файлов boundary = glob.glob('boundary/*.shp') raster = glob.glob('input_raster/*.tif')

Шаг 3. Проверка наличия пути с векторным файлом

# проверка наличия пути с векторным файлом if not boundary: print('Отсутствует векторный файл') sys.exit()

Шаг 4. Проверка наличия пути с растровым файлом

# проверка наличия пути с растровым файлом if not raster: print('Отсутствует растровый файл') sys.exit()

Шаг 5. Читаем векторный файл

# читаем векторный файл boundary = gpd.read_file(boundary[0]) if len(boundary) == 0: print('В векторном файле отсутствуют данные') sys.exit()

Шаг 6. Производим фильтрацию по типам геометрии

# производим фильтрацию по типам геометрии boundary = boundary[((boundary.geometry.geom_type == 'Polygon') | (boundary.geometry.geom_type == 'MultiPolygon')) & (~boundary.geometry.astype(str).str.contains("EMPTY"))] boundary.reset_index(drop=True, inplace=True)

Шаг 7. Фильтрация векторного слоя по атрибуту «title» для исключения пустых ячеек

# фильтрация векторного слоя по атрибуту 'title' для исключения пустых значений boundary = boundary[~boundary.title.isna()]

Шаг 8. Проверка наличия столбца «title» в векторном слое

# проверка наличия столбца 'title' в векторном слое if not 'title' in boundary.columns: print('Отсутствуют непустые значения "title"') sys.exit()

Шаг 9.Читаем растровый слой

# читаем растровый файл with rasterio.open(raster[0]) as src: # еслив векторном слое есть объекты с непустым значением атрибута 'title' if len(boundary)>0: # пробегаемся по геометриям векторном слое for row in range(len(boundary)): shapes = gpd.GeoSeries(boundary.loc[row, 'geometry']) boundary = boundary[~boundary.title.isna()] title = boundary.loc[row, 'title'] try: # обрезаем растровый слой по границам геометрий из векторного слоя out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True) out_meta = src.meta out_meta.update({"driver": "GTiff", "height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform}) # сохраняем "разрезанные" растры в файл с наименованием соответствущим значению атрибуту title в векторном слое with rasterio.open(f"mask/{title}.tif", "w", **out_meta) as dest: dest.write(out_image) except ValueError: print('Проверьте соответсвие систем координат растра и вектора') break

Готовый скрипт по обрезке растровых файлов с помощью библиотек Python выглядит так:

# импорт библиотек import sys import glob import geopandas as gpd import rasterio import rasterio.mask # указываем путь до исходных файлов boundary = glob.glob('boundary/*.shp') raster = glob.glob('input_raster/*.tif') # проверка наличия пути с векторным файлом if not boundary: print('Отсутствует векторный файл') sys.exit() # проверка наличия пути с растровым файлом if not raster: print('Отсутствует растровый файл') sys.exit() # читаем векторный файл boundary = gpd.read_file(boundary[0]) if len(boundary) == 0: print('В векторном файле отсутствуют данные') sys.exit() # производим фильтрацию по типам геометрии boundary = boundary[((boundary.geometry.geom_type == 'Polygon') | (boundary.geometry.geom_type == 'MultiPolygon')) & (~boundary.geometry.astype(str).str.contains("EMPTY"))] boundary.reset_index(drop=True, inplace=True) # фильтрация векторного слоя по атрибуту 'title' для исключения пустых значений boundary = boundary[~boundary.title.isna()] # проверка наличия столбца 'title' в векторном слое if not 'title' in boundary.columns: print('Отсутствуют непустые значения "title"') sys.exit() # читаем растровый файл with rasterio.open(raster[0]) as src: # еслив векторном слое есть объекты с непустым значением атрибута 'title' if len(boundary)>0: # пробегаемся по геометриям векторном слое for row in range(len(boundary)): shapes = gpd.GeoSeries(boundary.loc[row, 'geometry']) boundary = boundary[~boundary.title.isna()] title = boundary.loc[row, 'title'] try: # обрезаем растровый слой по границам геометрий из векторного слоя out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True) out_meta = src.meta out_meta.update({"driver": "GTiff", "height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform}) # сохраняем "разрезанные" растры в файл с наименованием соответствущим значению атрибуту title в векторном слое with rasterio.open(f"mask/{title}.tif", "w", **out_meta) as dest: dest.write(out_image) except ValueError: print('Проверьте соответсвие систем координат растра и вектора') break

Теперь вы можете использовать готовые скрипты в собственных ГИС-проектах:

— Автоматизировать анализ пространственных данных, изменения координат и создания карт, внедряя в работу стандартные скрипты;

— Использовать встроенные библиотеки в QGIS и на основе пакетного метода обрабатывать растровые и векторные данные с сохранением в целевые директории;

— С помощью библиотек Python — Rasterio и Geopandas сможете автоматически «нарезать» растровые файлы.

Обзор готовых скриптов Python для работы с геоданными

Python и типовые скрипты

Что нужно знать ГИС-специалисту, чтобы писать скрипты

Готовый Скрипт — выгрузка растровых слоев из QGIS

Обрезка растров с помощью библиотеки Python — Rasterio

Вывод