Просмотр кода страницы Работа с растрами при помощи GDAL и Python

{{Статья|Опубликована|gdal-python}}
{{Аннотация|Руководство по использованию GDAL Python API, частично основано на «GDAL API tutorial»}}

GDAL — свободная библиотека для работы с растровыми данными. Утилиты командной строки, входящие в состав библиотеки широко используются для выполнения разнообразных задач.

Благодаря наличию развитого API можно работать с функциями GDAL из многих языков программирования. та статья описывает работу с GDAL API через Python. При написании статьи использовались материалы [http://gdal.org/gdal_tutorial.html GDAL API tutorial].

=== Подготовка ===

Будем исходить из того, что все программное обеспечение установлено при помощи установщика OSGeo4W. Для использования GDAL совместно с Python необходимо наличие соответствующий байндингов (bindings). В стеке программ OSGeo4W нужный пакет называется gdal-python (для GDAL 1.5.x) или gdal16-python (для GDAL 1.6.x). Рекомендую использовать версию 1.6.x. Кроме того, для выполнения различных операций с растром, нам потребуется пакет numpy.

GDAL Python API состоит из пяти основных модулей и пяти дополнительных (существуют для совместимости со старыми версиями):

* gdal — Python интерфейс к библиотеке GDAL
* ogr — Python интерфейс к библиотеке OGR
* osr — работа с системами координат
* gdal_array — вспомогательные функции
* gdalconst — константы

Подключить их можно командами

<syntaxhighlight lang="python">
# основные
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo import gdal_array
from osgeo import gdalconst

# версии для совместимости. Будут удалены в версии 2.0
import gdal
import ogr
import osr
import gdalnumeric
import gdalconst
</syntaxhighlight>

Если используется GDAL версии 1.5 и выше, рекомендуется использовать «основные» модули. А для случаев, когда необходимо использовать код написанный ранее можно сделать проверку

<syntaxhighlight lang="python">
try:
  from osgeo import gdal
except ImportError:
  import gdal
</syntaxhighlight>

В большинстве случаев достаточно подключить только модуль gdal.

=== Открытие файла ===

Для открытия растра используем функцию <tt>gdal.Open()</tt>, в качестве аргумента принимающую полный путь к растру и необзятельную константу, описывающую режим открытия. Если константа опущена, то подразумевается режим только для чтения. В случае успеха функция вернет объект <tt>GDALDataset</tt>, в противном случае — <tt>None</tt>.

<syntaxhighlight lang="python">
import osgeo.gdal as gdal
gdalData = gdal.Open( "/home/alex/test/input.tiff", GA_ReadOnly )
# или так
# gdalData = gdal.Open( "/home/alex/test/input.tiff" )

# проверяем все ли в порядке
if gdalData is None:
  print "ERROR: can't open raster"
  sys.exit( 1 )
</syntaxhighlight>

=== Получение информации ===

Объект <tt>GDALDataset</tt> позволяет получить информацию о количестве каналов растра, извлечь данные и метаданные. Все используемые функции имеют понятные названия, поэтому подробно их описывать не вижу смысла, вместо этого приведу небольшой пример:

<syntaxhighlight lang="python">
print "Driver short name", gdalData.GetDriver().ShortName
print "Driver long name", gdalData.GetDriver().LongName
print "Raster size", gdalData.RasterXSize, "x", gdalData.RasterYSize
print "Number of bands", gdalData.RasterCount
print "Projection", gdalData.GetProjection()
print "Geo transform", gdalData.GetGeoTransform()
</syntaxhighlight>

Последняя функция возвращает очередь из 6 чисел:

* X координата левого верхнего пикселя
* ширина пикселя
* вращение
* Y координата левого верхнего пикселя
* вращение
* высота пикселя

=== Извлечение растра/канала ===

Любой растр можно представить в виде массива, одноканальный в виде двумерного размерностью X на Y; многоканальный — в виде многомерного. GDAL позволяет получить доступ как ко всему растру сразу, так и к любому каналу. Результатом такой операции будет массив numpy соответсвующей размерности.

<syntaxhighlight lang="python">
# получаем весь растр целиком
raster = gdalData.ReadAsArray()
# получаем отдельный канал
gdalBand = gdalData.GetRasterBand( 1 )
band_1 = gdalBand.ReadAsArray()
</syntaxhighlight>

Кроме того, можно обрабатывать растр построчно или блоками

<syntaxhighlight lang="python">
line = gdalBand.ReadAsArray( xoffeset, yoffset, xsize, ysize )
</syntaxhighlight>

Здесь

* <tt>xoffset</tt> — смещение по X
* <tt>yoffset</tt> — смещение по Y
* <tt>xsize</tt> — размер блока по X
* <tt>ysize</tt> — размер блока по Y

Т.е. если оба смещения буду равны нулю, а <tt>xsize</tt> и <tt>ysize</tt> — высоте и ширине растра, то в переменной <tt>line</tt> мы получим первую строку массива. А изменяя значение <tt>yoffset</tt> можно последовательно пройтись по всем строкам массива.

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

<syntaxhighlight lang="python">
band_1 = gdalBand.ReadAsArray().astype( numpy.float32 )
</syntaxhighlight>

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

=== Операции с растром ===

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

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

<syntaxhighlight lang="python">
result = numpy.add( array1, array2 )
</syntaxhighlight>

Аналогично выполняется вычитание, умножение, деление и другие действия. Унарные операции тоже поэлементные, например после выполнения сроки

<syntaxhighlight lang="python">
r = sin( array1 )
</syntaxhighlight>

в массив <tt>r</tt> будут записаны значения синуса каждого элемента массива <tt>array1</tt>. Вот более сложный пример — поэлементная (попиксельная) обработка растра

<syntaxhighlight lang="python">
gdalData = gdal.Open( "/home/alex/test/input.tiff" )
# размер растра
xsize = gdalData.RasterXSize
ysize = gdalData.RasterYSize
# получаем растр в виде массива
raster = gdalData.ReadAsArray()
# перебираем все пиксели растра
for col in range( xsize ):
  for row in range( ysize ):
    # если значение пикселя равно 5, то меняем его на 10
    # иначе значение остается без изменений
    if raster[ row, col ] == 5:
      raster[ row, col ] = 10
</syntaxhighlight>

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

<syntaxhighlight lang="python">
temp1_bool = numpy.equal(raster,5)
numpy.putmask(raster,temp1_bool,10)
</syntaxhighlight>

=== Сохранение файла ===

Существует два способа сохранить растр: используя <tt>CreateCopy()</tt> или <tt>Create()</tt>. При использовании <tt>CreateCopy()</tt> необходимо указать растр-источник, параметры которого будут использованы при создании нового растра. В случае использования <tt>Create()</tt> необходимо вручную сформировать метаданные и выполнить запись растра.

Следует помнить, что не все драйверы поддерживают метод <tt>Create()</tt>, поэтому вначале необходимо проверить, есть ли поддержка нужного метода в драйвере:

<syntaxhighlight lang="python">
format = "GTiff"
driver = gdal.GetDriverByName( format )
metadata = driver.GetMetadata()
if metadata.has_key( gdal.DCAP_CREATE ) and metadata[ gdal.DCAP_CREATE ] == "YES":
  pass
else:
  print "Driver %s does not support Create() method." % format
  sys.exit( 1 )
# аналогично выполняется проверка для CreateCopy
if metadata.has_key( gdal.DCAP_CREATECOPY ) and metadata[ gdal.DCAP_CREATECOPY ] == "YES":
  pass
else:
  print "Driver %s does not support CreateCopy() method." % format
  sys.exit( 1 )
</syntaxhighlight>

Некоторые драйверы могут работать в режиме только чтения и не поддерживают ни один из этих методов.

При использовании метода <tt>CreateCopy()</tt> вся необходимая информация берется из «эталонного» растра, но есть возможность задать специфичные для формата параметры

<syntaxhighlight lang="python">
# эталонный растр
inputData = gdal.Open( "/home/alex/test/input.tiff" )
# создаем свой растр "по образу и подобию" inputData
outputData = driver.CreateCopy( "/home/alex/test/output.tiff", inputData, 0 )
# то же самое, но заданы дополительные параметры
# outputData = driver.CreateCopy( "/home/alex/test/output.tiff", inputData, 0, [ 'TILED=YES', 'COMPRESS=PACKBITS' ] )
# закрываем датасеты и освобождаем память
inputData = None
outputData = None
</syntaxhighlight>

При использование <tt>Create()</tt> необходимо вручную задать количество каналов, размеры растра и желаемый тип данных (байт, длинное целое, число с плавающей запятой...)

<syntaxhighlight lang="python">
# размеры растра
cols = 512
rows = 512
# количество каналов
bands = 1
# тип данных
dt = gdal.GDT_Byte
# создаем растр
outData = driver.Create( "/home/alex/test/out.tiff", cols, rows, bands, dt )
</syntaxhighlight>

После того, как растр создан можно добавить информацию о проекции.

<syntaxhighlight lang="python">
outData.SetProjection( proj )
outData.SetGeoTransform( transform )
</syntaxhighlight>

И записать данные в канал

<syntaxhighlight lang="python">
outData.GetRasterBand( 1 ).WriteArray( raster )
# после записи данных закрываем датасет
outData = None
</syntaxhighlight>

Если каналов несколько, записывать удобнее в цикле

<syntaxhighlight lang="python">
bands = 3
for i in range( bands ):
  outData.GetRasterBand( i + 1 ).WriteArray( raster[ i ] )
outData = None
</syntaxhighlight>

Само собой, массив raster должен иметь столько измерений, сколько каналов в растре (или больше).

=== Пример. Простой растровый калькулятор ===

Вооружившись этими сведениями попробуем написать простенький скрипт на Python для сложения двух каналов растра. Скрипт принимает два параметра: исходный растр и путь для сохранения результата.

<syntaxhighlight lang="python">
# -*- coding: utf-8 -*-

#!/usr/bin/env python

import sys
import numpy
import osgeo.gdal as gdal

# функция проверяющая растр на многоканальность
def isMultiband( path ):
  gdalData = gdal.Open( path )
  if gdalData.RasterCount < 2:
    print "ERROR: raster must contain at least 2 bands"
    return False
  return True

# извлекает из заданного растра канал с заданным номером
def bandAsArray( path, bandNum ):
  gdalData = gdal.Open( path )
  gdalBand = gdalData.GetRasterBand( bandNum )
  array = gdalBand.ReadAsArray().astype( numpy.float32 )
  gdalBand = None
  gdalData = None
  return array

# сохраняет массив array как растр с именем outPath в формате
# GeoTiff. Данные о проекции берутся из растра etalonPath
def saveRaster( outPath, etalonPath, array ):
  gdalData = gdal.Open( etalonPath )
  projection = gdalData.GetProjection()
  transform = gdalData.GetGeoTransform()
  xsize = gdalData.RasterXSize
  ysize = gdalData.RasterYSize
  gdalData = None

  format = "GTiff"
  driver = gdal.GetDriverByName( format )
  metadata = driver.GetMetadata()
  if metadata.has_key( gdal.DCAP_CREATE ) and metadata[ gdal.DCAP_CREATE ] == "YES":
    outRaster = driver.Create( outPath, xsize, ysize, 1, gdal.GDT_Float32 )
    outRaster.SetProjection( projection )
    outRaster.SetGeoTransform( transform )
    outRaster.GetRasterBand( 1 ).WriteArray( array )
    outRaster = None
  else:
    print "Driver %s does not support Create() method." % format
    return False

if __name__ == '__main__':
  args = sys.argv[ 1: ]
  inPath = args[ 0 ]
  outPath = args[ 1 ]
  if isMultiband( inPath ):
    band1 = bandAsArray( inPath, 1 )
    band2 = bandAsArray( inPath, 2 )
    res = numpy.add( band1, band2 )
    if not saveRaster( outPath, inPath, res ):
      print "ERROR: saving failed"
      sys.exit( 1 )
</syntaxhighlight>

=== Ссылки по теме ===

# [http://gis-lab.info/qa/ogr-python.html Работа с векторными данными при помощи OGR и Python]
# [http://gdal.org/gdal_tutorial.html GDAL API tutorial]
# [http://gdal.org/gdal_8h.html GDAL C API]
# [http://trac.osgeo.org/gdal/wiki/GdalOgrInPython GDAL/OGR in Python]
# [http://www.scipy.org/Tentative_NumPy_Tutorial Tentative NumPy Tutorial]
# [http://docs.scipy.org/doc/ NumPy and SciPy documentation]