Начало работы с растровыми геоданными средствами GDAL-Python

Введение в растровую модель геоданных и работу с ней средствами GDAL в Python. Содержание статьи:

• Концепция растровой модели геоданных

• Примеры растровых геоданных

• Свойства растровых геоданных

• Хранение растровых геоданных

• Знакомство с GDAL

• Чтение данных и их свойств

• Базовые манипуляции и преобразования

• Запись данных

• Комментарий к статье

Концепция растровой модели геоданных

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

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

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

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

Примеры растровых геоданных

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

1. Численные поля географических параметров: высота, средняя температура воздуха за определенный период, влажность почвы, концентрация фитопланктона на участке акватории и т.п.

   

2. Более абстрактные численные поля: интерполированное распределение цен на недвижимость, расстояние до ближайшего населенного пункта и т.д.

   

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

   

4. Кодовые значения, опосредованно описывающие категорию территории: типы землепользования, породы леса, бинарные классификации (например, вода/не вода).

   

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

   

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

Основные свойства

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

• Непрерывность внутри пространственного охвата. К какой бы точке внутри растрового набора геоданных мы не обратились, мы обязательно попадём в одну из его ячеек, и, соответственно, получим значения. Поэтому растровая модель незаменима при описании непрерывных процессов (например, метеорологических - температура воздуха ведь определена везде, непрерывно), и это сильно отличает её от объектной векторной модели.

  

• Гомогенность каждой отдельной ячейки. Каждая ячейка описывает не точку местности, а область на ней. Понятно, что внутри этой области, какой бы небольшой она ни была, в реальном мире идеальной однородности не будет, а описывать мы её собрались единым набором значений. Мы вынуждены принимать эту условность и всегда останавливаться на том или ином уровне обобщения.

В практическом смысле для нас важнее технические свойства георастров. То, каким набором параметров они описываются:

1. Размер матрицы, то есть количество её столбцов и строк.

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

3. Тип данных - какие значения хранятся в ячейках, например целочисленные или вещественные.

4. Географическое положение - в каких географических координатах лежит эта матрица.

5. Система координат - в какой системе координат описано это географическое положение.

6. Пространственное разрешение - размер каждой ячейки в заданной системе координат.

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

• Значение, используемое для того, чтобы помечать ячейки с отсутствием информации (NoData).

• Тип сжатия, если оно поддерживается.

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

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

И так далее. Мы не будем подробно рассматривать всё, сконцентрируемся на основном.

Хранение растровых геоданных

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

• GeoTIFF. Расширение популярного формата для хранения изображений TIFF, позволяющее дополнительно к его возможностям хранить все важные метаданные растровых геоданных (системы координат и прочее). Можно считать основным стандартом и базовой рекомендацией, формат очень прозрачный и универсальный. С некоторыми оговорками позволяет использовать данные одновременно в геоинформационном ПО и в графических редакторах. Имеет очень важную модификацию - Cloud Optimized GeoTIFF ^[1], которая помогает эффективно использовать растровые геоданные напрямую в веб-приложениях.

• JPEG, PNG, GIF - очень распространенные, всем знакомые форматы для работы с изображениями. Для хранения метаданных о географической привязке матриц используется специальный формат вспомогательных файлов, worldfile ^[2]. На мой взгляд, причин использовать такую конструкцию сегодня нет, но встретить такие данные можно довольно часто.

Конечно, есть сотни других форматов для хранения георастров, они рождаются и умирают вместе с различным программным обеспечением, организациями и стандартами. В частных обстоятельствах вы можете столкнуться много с чем ещё, например учёные любят хранить растры в NetCDF или HDF, а некоторые предпочитают работать с человекочитаемыми текстовыми файлами ASCII Grid. Благодаря GDAL, о котором мы будем говорить дальше, зоопарк возможных форматов не является большой проблемой.

Знакомство с GDAL

GDAL ^[3] (Geospatial Data Abstraction Library) это одна из ключевых технологий современной индустрии геоинформационных технологий. Важнейшая работа, которую она выполняет в контексте растров, это создание для них универсальной абстрактной модели, в которую она помогает прочитать десятки форматов, и во многие из которых умеет её записать.

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

Помимо возможностей чтения и записи в GDAL также входит мощный набор методов и утилит ^[4] преобразования и обработки растровых геоданных. Они эффективно решают все базовые задачи, связанные с системами координат, трансформациями, ресемплингом и так далее. Их широко используют на самых разных уровнях (в бэкенде приложений работающих с геоданными, в пайплайнах моделей машинного обучения, в научном анализе), и мы рассмотрим главные из них.

GDAL в основном написан на C/C++ и имеет хорошо задокументированное API ^[5], а также оболочки для разных языков программирования. Мы будем использовать GDAL в Python ^[6]. Утилиты GDAL также доступны в виде консольных утилит, которые можно запускать прямо из командной строки вашей операционной системы. GDAL лицензирован ^[7]по свободной лицензии MIT.

Установить GDAL для использования в Python можно через pip, conda, или использовать один из многочисленных docker-контейнеров, посвященных этой теме.

pip install gdal

или

conda install gdal

Иногда стандартная установка приводит к проблемам из-за компилируемых компонентов, поэтому я всегда советую устанавливать GDAL в отдельный venv или среду conda, или с использованием другого способа разделения пакетов Python.

Начинаем писать код! В репозитории ^[8]вы найдёте его в формате jupyter notebook вместе с используемыми для примера данными.

Чтение данных и их свойств

В папке data у нас есть несколько компактных примеров разных растров:

• multiband_imagery.tif - четырех-канальный растр, фрагмент спутникового снимка Sentinel-2, в каналах 3 видимых диапазона (Blue, Green, Red) и ближний инфракрасный (IR).

• land_cover_types.zip - архив с одноканальным растром с категориями типов земель, на основе данных MODIS MCD12Q1. В нём разными числами закодированы разные типы местности (лес, застройка и т.п.)

• digital_elevation_model.tif - одноканальный растр, цифровая модель рельефа, в каждой ячейке - высота над уровнем моря в метрах.

• topomap.tif - фрагмент топографической карты, RGB изображение с географической привязкой.

Прочитаем основные метаданные digital_elevation_model.tif. Для начала импортируем gdal, он находится за osgeo. Также я сразу импортирую numpy и matplotlib чтобы использовать их для обработки и визуализации прочитанных данных.

from osgeo import gdal
import matplotlib.pyplot as plt
import numpy as np

Теперь читаем файл. После чтения в переменной dem_ds будет объект класса osgeo.gdal.Dataset, то самое универсальное абстрактное представление растров GDAL-ом, для которого доступно множество методов. Сразу прочитаем все основные метаданные, о которых мы говорили выше:

dem_ds = gdal.Open('data/digital_elevation_model.tif', gdal.GA_ReadOnly)
print ('Basic metadata')
print ('Bands count: %s' % dem_ds.RasterCount)
print ('Data type: %s' % dem_ds.RasterCount)
print ('Rows: %s | Cols: %s' % (dem_ds.RasterYSize, dem_ds.RasterXSize))
print ('GeoTransform: %s' % str(dem_ds.GetGeoTransform()))
print ('Coordinate reference system: %s' % dem_ds.GetProjection())

> Basic metadata
> Bands count: 1
> Rows: 632 | Cols: 616
> GeoTransform: (18.70875, 0.0002777777775974056, 0.0, 43.559861111, 0.0, -0.00027777777689872693)
> Coordinate reference system: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

Итак, перед нами одноканальный растр с размером 616х632, это ясно. А вот с тем, что вернули GetGeoTransform() и GetProjection() нужно разобраться отдельно.

GetGeoTransform() возвращает кортеж из 6 элементов, который описывает географическое положение растра. Вот что значат его значения в порядке следования:

1. Координата X верхнего левого угла растра (18.70875)

2. Горизонтальный размер ячейки в системе координат растра (0.00027)

3. Первый параметр поворота растра, в реальности растры почти всегда ориентированы на “север” (в терминах своей системы координат) и это значение равно 0

4. Координата Y верхнего левого угла растра (43.559861111)

5. Второй параметр поворота растра, аналогично на практике всегда 0

6. Вертикальный размер ячейки в системе координат растра (-0.00027)

Таким образом GeoTransoform определяет, откуда растр начинается и каков пространственный размер каждой его ячейки. Фактически из этих данных можно вычислить общий охват растра (зная размер матрицы), а также пространственные границы каждой конкретной ячейки.

GetProjection() возвращает описание системы координат растра в стандартном формате OGC WKT CRS ^[9]. Именно в этой СК мы получили все данные от GeoTransform. Благодаря этим двум методам мы получаем полную информацию о географической принадлежности каждой ячейки растра.

Следующий шаг - достать сами данные. Мы можем достать матрицу конкретного канала, либо достать матрицу растра целиком. В случае одноканального растра результат будет одинаков.

# Растр целиком
dem_ds_data = dem_ds.ReadAsArray()

# Данные конкретного канала
dem_ds_band1 = dem_ds.GetRasterBand(1)
dem_ds_band1_data = dem_ds_band1.ReadAsArray()

print (dem_ds_data.shape)
print (dem_ds_band1_data.shape)

> (632, 616)
> (632, 616)

На выходе - numpy матрица, готовая к работе. На самом деле во многих сценариях на этом работа с GDAL и заканчивается - достали из растра матрицу значений и всё, понесли её куда нужно. Из неё можно достать значение в конкретной ячейке или получить срез:

# значение высоты в ячейке 50, 50
print (dem_ds_data[50,50])

# срез значений высот
print (dem_ds_data[40:48,40:48])

> 754
> [[752 762 773 783 792 799 806 813]
  [747 757 767 777 785 792 799 809]
  [746 754 763 772 780 789 796 803]
  [741 750 759 769 779 787 794 799]
  [736 745 753 763 773 780 786 790]
  [727 738 747 755 762 770 777 781]
  [720 731 741 747 753 761 769 775]
  [713 725 735 741 746 753 761 768]]

Или визуализировать целиком:

# визуализация целиком
plt.imshow(dem_ds_data)
plt.colorbar()
plt.show()

Давайте сразу поупражняемся в вычислениях на основе метаданных. Например, вычислим координаты ячейки 50;50.

# вычисление координат конкретной ячейки (50;50) на основе GeoTransform
gt = dem_ds.GetGeoTransform()
x_coord = gt[0] + 50*gt[1]
y_coord = gt[3] + 50*gt[5]
print (x_coord, y_coord)

> 18.722638888879867 43.54597222215506

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

# средняя высота всего участка территории
print (dem_ds_data.mean())

# маска для ячеек со значением > 1500
high_elevations = dem_ds_data[dem_ds_data > 1500]
print (high_elevations.shape[0])

> 867.7228289906296
> 8111

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

Теперь прочитаем многоканальный растр, multiband_imagery.tif.

imagery_ds = gdal.Open('data/multiband_imagery.tif', gdal.GA_ReadOnly)
print ('Bands count: %s' % imagery_ds.RasterCount)

imagery_ds_data = imagery_ds.ReadAsArray()
print (imagery_ds_data.shape)

Bands count: 4
(4, 177, 191)

Склеим и покажем цветное трехканальное изображение в естественных цветах (RGB) Чтобы изображение получилось красивым, исходные значения можно нормализовать в диапазон [0;1]. Поскольку мы работаем с numpy, это не составит труда.

# описываем функцию, которая будет нормализовать значения канала в диапазон от 0 до 1
def normalize(input_band):
    min_value, max_value = input_band.min()*1.0, input_band.max()*1.0
    return ((input_band*1.0 - min_value*1.0)/(max_value*1.0 - min_value))

# собираем матрицу из нормализованных каналов
rgb_normalized = np.dstack([normalize(imagery_ds_data[2]),
                            normalize(imagery_ds_data[1]),
                            normalize(imagery_ds_data[0])])

plt.imshow(rgb_normalized)
plt.show()

И ещё произведем вычисления внутри многоканального растра - посчитаем популярный вегетационный индекс NDVI. Считается он простой формулой нормализованного разностного индекса, и, опять же, благодаря numpy всё очень просто.

# Вычисляем NDVI
ndvi = (imagery_ds_data[3] - imagery_ds_data[2]) / (imagery_ds_data[3] + imagery_ds_data[2])

# применяем жёлто-зелёную шкалу от 0 до 1. Чем зеленее, тем больше NDVI 
plt.imshow(ndvi, cmap='YlGn', vmin=0, vmax=1)
plt.colorbar()
plt.show()

Подтянув scipy ^[10], scikit-learn ^[11] и другие инструменты Python, мы можем начать многое извлекать из прочитанных данных - кластеризовать, сегментировать, и классифицировать спутниковые данные, поверх данных цифровой модели рельефа считать различные геоморфометрические производные, отправлять это всё в модели машинного обучения и так далее.

В завершение раздела посмотрим ещё на кое-что - как без лишних манипуляций открывать растровые геоданные напрямую из архивов и из сетевых источников подходящего типа. В этом нам помогут т.н. GDAL Virtual File Systems ^[12].

# данные из архива (vsizip)
ds_from_zip = gdal.Open('/vsizip/data/land_cover_types.zip/land_cover_types.tif')
ds_from_zip_data = ds_from_zip.ReadAsArray()
print ('Срез данных в растре из архива:')
print (ds_from_zip_data[100:105,100:105])

# данные из сетевого хранилища (vsicurl)
ds_from_url = gdal.Open('/vsicurl/https://demo.nextgis.com/api/resource/7220/cog')
print ('nРазмер сетевого растра:')
print(ds_from_url.RasterXSize)
print(ds_from_url.RasterYSize)

# Читаем конкретные значения, не скачивая весь растр
print('4 пикселя начиная с адреса 10000, 10000. Все каналы:')
print(ds_from_url.ReadAsArray(10000,10000,2,2))

> Срез данных в растре из архива:
> [[20 20 20 25 25]
  [20 20 20 20 25]
  [20 25 25 25 20]
  [20 20 20 20 20]
  [10 10 10 10 20]]

> Размер сетевого растра:
> 28416
> 24320
> 4 пикселя начиная с адреса 10000, 10000. Все каналы:
> [[[101  99]
   [ 99  97]]

  [[ 98  96]
   [ 98  96]]

  [[ 93  91]
   [ 94  92]]

  [[255 255]
   [255 255]]]

Второй пример особенно интересен. Растровый набор геоданных, хранящийся на сервере NextGIS Web, довольно большой - четырёхканальный с размером 28416х24320, и скачивать его целиком может быть неприятно. GDAL позволяет быстро получить метаданные и загрузить только ту часть данных, которая нам нужна. Обратите внимание также на то, как мы использовали метод ReadAsArray() - на этот раз мы указали конкретный срез пикселей в аргументах.

Теперь посмотрим, что предлагает сам GDAL в части манипуляций с данными.

Базовые манипуляции и преобразования

Мы упоминали утилиты GDAL - мощные разносторонние инструменты для преобразований растров, который можно запускать из консоли. Многие из них доступны также в качестве методов в API GDAL. Самые мощные из них это gdal_translate ^[13], заточенный на конвертацию с попутными преобразованиями, и gdalwarp ^[14], универсальный инструмент поддерживающий огромное число типов преобразования растров.

Эти методы позволяют в одно простое действие выполнять целые цепочки сложных преобразований, например “Перепроецируй этот растр в систему координат Pseudo Mercator (EPSG:3857) ^[15], затем вырежи определенный кусок, загруби пространственное разрешение в два раза и сохрани результат в файл в формате GeoTIFF со сжатием LZW.” Посмотрим, как выглядят подобные операции при работе из Python.

Чуть выше мы хотели посчитать площадь территории с высотой > 1500 м., но не стали этого делать из-за особенностей системы координат данных. Давайте вернёмся к этому примеру, трансформируем данные в подходящую систему координат и вычислим совокупную площадь высокогорных пикселей.

# возвращаемся к набору dem_ds и перепроецируем его
dem_ds_transformed = gdal.Warp('',dem_ds, format='MEM', dstSRS='EPSG:32634')

# И всё готово, в переменной dem_ds_transformed хранящийся в памяти gdal.Dataset готовый к работе
# читаем и смотрим на GeoTransform
gt = dem_ds_transformed.GetGeoTransform()
print ('Новый GeoTransform:')
print (gt)
# читаем перепроецированную матрицу значений и применяем маску
dem_ds_transformed_data = dem_ds_transformed.ReadAsArray()
high_elevations = dem_ds_transformed_data[dem_ds_transformed_data > 1500]

# Вычисляем площадь, умножая количество подходящих пикселей на площадь отдельного пикселя. Её вычисляем на основе данных Geotransform и приводим к км2
total_area = high_elevations.shape[0] * (abs(gt[1])*abs(gt[5])) / 1000000
print ('Общая площадь районов с высотой > 1500 м: %s км2' % round(total_area,3))

del dem_ds_transformed 

> Новый GeoTransform:
> (314409.25304288446, 27.10218162365495, 0.0, 4825539.987535062, 0.0, -27.10218162365495)
> Общая площадь районов с высотой > 1500 м: 5.626 км2

В этом коде можно обсудить сразу несколько вещей. Во-первых, сам запуск gdal.Warp. В качестве первого аргумента этот метод принимает путь до результирующего файла. Мы оставили его пустым, и затем указали специальный формат MEM - такая комбинация очень полезна и позволяет работать с производными растрами в памяти, не прибегая к использованию временных промежуточных файлов. Запуск мог бы выглядеть и так:

dem_ds_transformed = gdal.Warp('data/transformed_dem.tif',dem_ds, format='GTiff', dstSRS='EPSG:32634')

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

dem_ds_transformed = gdal.Warp('data/transformed_dem.tif',dem_ds, format='GTiff', dstSRS='EPSG:32634')
del dem_ds_transformed 

Ещё обратите внимание на новые значения в GeoTransform. Они сильно изменились после перепроецирования - теперь координаты верхнего левого угла растра это большие числа в метрах, и размер пикселей также стал более понятным: ~27 метров вместо изначальных долей градуса. Система координат EPSG:32634 ^[16] была выбрана как минимально искажающая площади объектов в регионе данных.

Второе упражнение - загрубим пространственное разрешение топографической карты из файла topomap.tif и сохраним её в PNG с worldfile рядом. На этот раз используем gdal.Translate. И для разнообразия мы вовсе не будем читать исходный растр в объект GDAL, а посмотрим, как просто направить в метод gdal.Translate файл.

# Посмотрим на исходные данные как на RGB картинку
import matplotlib.image as mpimg
image = mpimg.imread('data/topomap.tif')
plt.imshow(image)
plt.show()

# отправим их в gdal.Translate с двумя заданиями - конвертация формата и загрубление пространственного разрешения
topomap_converted = gdal.Translate('data/topomap_converted.png','data/topomap.tif',format='PNG',xRes=20,yRes=20, creationOptions=["WORLDFILE=YES"])
del topomap_converted

# Смотрим на новый растр
image = mpimg.imread('data/topomap_converted.png')
plt.imshow(image)
plt.show()

По картинке видно, что мы сильно уменьшили разрешение, вручную указав его с помощью аргументов xRes и yRes. Также стоит обратить внимание на аргумент creationOptions=["WORLDFILE=YES"], в котором мы передали инструкцию по созданию Worldfile. Так в методы gdal передаются отдельные дополнительные параметры, которые в документации к утилитам следуют за ключом -co.

У методов Translate и Warp много возможностей, хорошие способы их изучить - посмотреть документацию к одноименным утилитам (translate, warp), или обратиться к описанию метода через gdal.WarpOptions и gdal.TranslateOptions

gdal.WarpOptions?

> Signature:
  gdal.WarpOptions(
    options=None,
    format=None,
    outputBounds=None,
    outputBoundsSRS=None,
    xRes=None,
    yRes=None,
    targetAlignedPixels=False,
    width=0,
    height=0,
    srcSRS=None,
    dstSRS=None,
    coordinateOperation=None,
    ...

При ресемплинге может быть важен используемый алгоритм интерполяции (например, nearest neighbour для категориальных растров или cubic для непрерывных). Задумавшись об этом, несложно найти способом выше аргумент resampleAlg, в который можно записать одно из значений: near|bilinear|cubic|cubicspline|lanczos|average|rms|mode|min|max|med|q1|q3|sum. Подобных настроек доступно много.

Из других наиболее ходовых методов обработки данных: gdal.Rasterize для превращение векторных данных в растровые, gdal.Polygonize для превращения растровых данных в векторные, gdal.BuildVRT для создания виртуальных растров (о которых стоит в будущем поговорить отдельно).

Запись данных

Ещё один небольшой, но важный раздел - это создание новых растров вручную. Покажем механизм на примере расчётов NDVI из данных спутниковой съемки выше. Вернемся к этому коду и реализуем запись матрицы значений NDVI в новый файл.

# открываем данные и считаем NDVI
imagery_ds = gdal.Open('data/multiband_imagery.tif', gdal.GA_ReadOnly)
imagery_ds_data = imagery_ds.ReadAsArray()
ndvi = (imagery_ds_data[3] - imagery_ds_data[2]) / (imagery_ds_data[3] + imagery_ds_data[2])

# записываем матрицу в новый файл
# Выбираем драйвер (формат выходного файла)
driver = gdal.GetDriverByName("GTiff")

# У нас будет 1 канал, только NDVI
band_count = 1

# Выбираем тип данных
data_type = gdal.GDT_Float32

# Создаём файл на основе выбранного драйвера, наследуя размер у исходного набора
dataset = driver.Create('data/ndvi.tif', imagery_ds.RasterXSize, imagery_ds.RasterYSize, band_count, data_type)

# Заимствуем GeoTransform и систему координат у исходного растра (новый такой же, т.к. произведен в его пространственном домене)
dataset.SetProjection(imagery_ds.GetProjection())
dataset.SetGeoTransform(imagery_ds.GetGeoTransform())

# Записываем сами данные
dataset.GetRasterBand(1).WriteArray(ndvi)

del dataset

Процесс достаточно прозрачен. Теперь у нас есть новый файл с результатами расчётов!

Комментарий к статье

При углублении работы с GDAL появляется много нюансов и деталей. Например:

• Многими параметрами процессов обработки управляют специальные переменные окружения ^[17].

• Концепция виртуальных растров ^[18] (когда мы храним описание цепочки манипуляций, а не их результат, а результат вычисляется на лету при обращении) в ряде сценариев может быть очень полезна и достойна отдельного материала.

• Альтернативой GeoTransform для определения пространственного положение растра могут быть GCP (Ground Control Points), когда координаты заданы для конкретных пикселей, а положение остальных вычисляется интерполяцией на основе этих точек. Применение к такому растру Warp без параметров автоматически переведёт GCP в GeoTransform.

• Над GDAL в части работы с растрами есть много упрощающих разные процессы разработок. Из точно стоящих внимания - rasterio ^[19].

Можно продолжать и продолжать. Буду рад обсудить в комментариях наиболее интересные вам темы для дальнейших публикаций.

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

Это открытый проект, и если он приносит вам пользу, поддержите его.

GDAL GitHub ^[20] | GDAL Sponspring FAQ ^[21]