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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   

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

   

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

   

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

   

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

   

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

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

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

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

  

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

pip install gdal

или

conda install gdal

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

Начинаем писать код! В репозитории вы найдёте его в формате 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. Именно в этой СК мы получили все данные от 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, scikit-learn и другие инструменты Python, мы можем начать многое извлекать из прочитанных данных - кластеризовать, сегментировать, и классифицировать спутниковые данные, поверх данных цифровой модели рельефа считать различные геоморфометрические производные, отправлять это всё в модели машинного обучения и так далее.

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

# данные из архива (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, заточенный на конвертацию с попутными преобразованиями, и gdalwarp, универсальный инструмент поддерживающий огромное число типов преобразования растров.

Эти методы позволяют в одно простое действие выполнять целые цепочки сложных преобразований, например “Перепроецируй этот растр в систему координат Pseudo Mercator (EPSG:3857), затем вырежи определенный кусок, загруби пространственное разрешение в два раза и сохрани результат в файл в формате 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 была выбрана как минимально искажающая площади объектов в регионе данных.

Второе упражнение - загрубим пространственное разрешение топографической карты из файла 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 появляется много нюансов и деталей. Например:

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

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

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

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

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

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

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

GDAL GitHub | GDAL Sponspring FAQ