Примеры использования инструментов GDAL

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/gdal-examples.html

Перечь примеров для справки

GDAL/OGR - библиотека для работы с географическими форматами данных. GDAL представляет собой набор утилит для обработки растровых данных, в то время, как OGR предназначена для работы с векторными форматами. В статье рассматриваются некоторые практические примеры применения утилит этой библиотеки для работы с растровыми данными.

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

Если у вас есть свои часто используемые примеры - присылайте автору или дописывайте прямо здесь.

Общие сведения

Утилиты GDAL предназначены для конвертации растровых данных из одного формата в другой и выполнения над ними различных операций. Установить утилиты GDAL для Windows можно с помощью OSGeo4W. Удобный визуальный интерфейс для утилит имеется в GDALTools для QGIS. В комплект GDAL входят следующие утилиты:

• gdal_calc.py - растровый калькулятор, арифметические операции с растрами;
• gdal-config - получить опции необходимые для создания ПО использующего GDAL;
• gdal_contour - получение изолиний по цифровым моделям рельефа (ЦМР);
• gdal_edit.py - редактировать информацию о растре;
• gdal_fillnodata.py - заполнение областей имеющих значение NODATA;
• gdal_merge.py - создание мозаик и композитных изображений;
• gdal_polygonize.py - векторизовать растр с получением полигонального слоя;
• gdal_proximity.py - расчитать растр близости;
• gdal_rasterize - растеризация векторных данных;
• gdal_retile.py - создать новый набор тайлов и/или перестроить пирамидные слои;
• gdal_grid - создание ЦМР из векторных данных;
• gdal_sieve.py - фильтрация осколочных объектов растра;
• gdal_translate - конвертация растров из формата в формат;
• gdal2tiles.py - создание тайловой структуры, KML и простого просмотровщика;
• gdaladdo - добавление пирамидных слоёв (overview);
• gdalbuildvrt - создание виртуального растра (VRT) из набора;
• gdalcompare.py - сравнение двух изображений;
• gdaldem - набор инструментов для анализа и визуализации ЦМР;
• gdalinfo - информация о растре;
• gdallocationinfo - запросы информации к растровыми файлам;
• gdalmanage - управление растровыми файлами (копирование, переименование, удаление и т.д.);
• gdalmove.py - трансформирование системы координат растра без ресэмплирования;
• gdalsrsinfo - показывается информацию по системе координат в разных форматах (WKT, PROJ.4 и др.);
• gdaltindex - построить индекс фрагментов (тайлов) MapServer;
• gdaltransform- трансформация координат;
• gdalwarp - трансформация изображения в новую систему координат;
• nearblack - конвертация черных/белых границ в нужное значение;
• pct2rgb.py - конвертация 8-битных изображений с палитрой в 24-битные RGB изображений;
• rgb2pct.py - конвертация 24-битных RGB изображений в 8-битные с палитрой;

Поддерживаемые форматы и используемые ключи можно узнать просто набрав в командной строке имя одной из утилит.

gdalinfo

В результате будет получена справка по использованию этой программы:

Usage: gdalinfo [--help-general] [-mm] [-stats] [-nogcp] [-nomd]
                [-noct] [-checksum] [-mdd domain]* datasetname

Версию GDAL можно посмотреть командой:

gdalinfo --version

Список форматов поддерживаемых утилитами GDAL можно посмотреть следующим образом:

gdalinfo --formats

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

• GRASS (ro): GRASS Database Rasters (5.7+)
• VRT (rw+): Virtual Raster
• GTiff (rw+): GeoTIFF
• HFA (rw+): Erdas Imagine Images (.img)
• AIG (ro): Arc/Info Binary Grid
• AAIGrid (rw): Arc/Info ASCII Grid
• JPEG (rw): JPEG JFIF
• MEM (rw+): In Memory Raster
• GIF (rw): Graphics Interchange Format (.gif)
• BMP (rw+): MS Windows Device Independent Bitmap
• DIMAP (ro): SPOT DIMAP
• PCIDSK (rw+): PCIDSK Database File
• SRTMHGT (rw): SRTMHGT File Format
• GMT (rw): GMT NetCDF Grid Format
• HDF4 (ro): Hierarchical Data Format Release 4
• HDF4Image (rw+): HDF4 Dataset
• ENVI (rw+): ENVI .hdr Labelled
• EHdr (rw+): ESRI .hdr Labelled

Конвертация

Извлечь три канала с номерами 1, 2, 3 в новый файл из исходного с перекомбинацией, в котором каналов может быть больше.

gdal_translate -b 3 -b 2 -b 1 output.tif input.tif

В результате в текущем каталоге появится результат 3-х канальный файл output.tif. Или в цикле для например 46-канального (win):

for /L %i in (1,1,46) DO gdal_translate -b %i input.tif output_%i.tif

Конвертация с обрезкой по заданным координатам:

gdal_translate -of GTiff -projwin 75.081940 57.250275 89.869980 49.083084 input.tif output.tiff

Конвертация с компрессией и созданием world-файла:

gdal_translate -co "COMPRESS=LZW" -co "worldfile=yes" input.tif output.tiff

Конвертация 16 битного одноканального растра в 8 битный:

gdal_translate -scale -ot Byte input_16bit.tif output.tif

Пакетная конвертация всех JPG в TIF (Windows):

for %i in (*.jpg) do gdal_translate %i %~ni.tif 

Стэки, композиты, мозаики

Создание композитного изображения из серии отдельных растров, каждый из которых в своем файле TIF. Разрешение выходного файла устанавливается по первому их растров. Таким образом, если первый канал 15 м, а остальные 30 м, то последние будут пересчитаны на 15 м. Чтобы указать, что каждый исходный растр должен попасть в отдельный слой (layer-stacking), а не мозаицирование, используется ключ -separate:

gdal_merge.py -o output.tif band1.tif band2.tif band3.tif band4.tif band5.tif -separate

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

gdal_merge -o altay.tif srtm_53_02.tif srtm_53_03.tif srtm_54_02.tif srtm_54_03.tif

А так можно указать разрешение выходного растра и то, что использовать определенное значение - не нужно (чтобы поля не перекрывали значащие части):

gdal_merge -n 0 -ps 0.00416 0.00416 -o output.tif in1.tif in2.tif in3.tif in4.tif

Работа с NODATA

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

gdalwarp -srcnodata -999 -dstnodata 0 input.tif output.tif

Если в исходном растре nodata записано как nan NoData Value=nan, то конвертировать лучше так:

gdalwarp -dstnodata nan input.tif output.tif

Как снять флаг NODATA? Допустим есть растр в котором часть пикселей имеют значение 0 и на них установлен флаг NODATA и нужно его убрать. Делается это так:

gdal_translate -a_nodata none input_with_NoData.tif output_without_NoData.tif

Работа с системами координат

Переназначение системы координат

Если вам нужно просто перепрописать систему координат, без пересчета самого растра:

gdal_edit -a_srs "EPSG:4326" input.tif

Если код неизвестен, можно указать описание системы координат в формате Proj.4:

gdal_edit -a_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" input.tif

Узнать строку описания системы координат можно следующим образом:

gdalsrsinfo template.tif

Конвертация с перепроецированием

gdalwarp позволяет не только конвертировать данные из одного формата в другой, но и одновременно произвести перепроецирование данных из одной системы координат в другую. Для этого используются параметры:

• -a_srs используется для указания системы координат (СК) для данных
• -s_srs используется для перезаписи информации о системе координат
• -t_srs перепроецирования данных в требуемую систему координат

В самом простейшем случае это делается следующим образом (ключ -tr указывает разрешение целевого растра):

gdalwarp.exe -tr 0.0083 0.0083 -t_srs "EPSG:4326" in.tif out.tif

Если EPSG-кода у СК конечного растра нет, то указать целевую СК можно так:

gdalwarp.exe -t_srs "+proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=45 +x_0=8500000 +y_0=0 +ellps=krass +units=m +towgs84=28,-130,-95,0,0,0,0 +no_defs" in.tif out.tif

Обрезка

Обрезка по векторному контуру c уменьшением охвата растра (реальная обрезка, а не просто заполнение ненужных областей значениями NODATA):

gdalwarp -cutline aoi.shp -crop_to_cutline input.tif output.tif

Работа с рельефом

Теневая отмывка рельефа:

gdaldem hillshade altay.tif altayhill.tif -z 5 -s 111120

Ключ -s 111120 используется для пересчета для растров сделанных в EPSG:4326 в метровые СК. Если исходник уже находится в проекции, то он не нужен.

Цветовая отмывка рельефа:

gdaldem color-relief altay.tif ramp.txt altay-color.tif

Пример файла ramp.txt:

5000 255 255 255
1000 168 112 0
650 198 165 48
400 229 218 97
200 218 229 97
0 112 168 0

Работа с цветовой таблицей

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

С использованием VRT

1. Создаем файл VRT:

gdal_translate -of VRT input.tif input.vrt

2. Открываем файл VRT в любом текстовом редакторе, находим:

<ColorInterp>Gray</ColorInterp>

и меняем на:

<ColorInterp>Palette</ColorInterp>
<ColorTable>
        <Entry c1="0" c2="0" c3="128"/>
        <Entry c1="0" c2="128" c3="0"/>
        <Entry c1="0" c2="255" c3="0"/>
        <Entry c1="153" c2="204" c3="0"/>
        <Entry c1="131" c2="66" c3="37"/>
        <Entry c1="51" c2="153" c3="102"/>
</ColorTable>

Каждая Entry - это цвет, кодированный в RGB и соответствующий по порядку значениям от 0 до 255. Если цветов меньше 255, то остальные будут добавлены автоматически с RGB = 0,0,0 (черный).

После того как таблица отредактирована, нужно сохранить ее обратно в растр:

gdal_translate input.vrt result.tif

Построение изолиний

Утилита gdal_contour используется для получения изолиний - линий равных значений по растровым данным. Полученные линии пересекают все пиксели с одинаковым значением, очерчивая при этом некоторую область. Чаще всего применяется для построения горизонталей рельефа из ЦМР.

Построение контуров с интервалом в 5 единиц (Единица указывается в единицах измерения исходного растра):

gdal_contour -i 5 mydem.tif contour.shp

Построение контуров из первого канала растра, с интервалом в 100 единиц начиная с 1200 и записью значения в поле elev:

gdal_contour -b 1 -a elev -i 100 -off 1200 mydem.tif contour.shp

Построение только контуров с фиксированными значениями 1000, 1100 и 1120 и выводом результата в таблицу PostGIS

gdal_contour -a elev -f PostgreSQL -fl 1000 1100 1120 -nln cont mydem.tif "PG:host=localhost user=iampg password=iampgpass dbname=iamgis"

netCDF

Конвертирование в GeoTIFF:

gdal_translate -of GTiff -b 1 NETCDF:precip.mon.mean.nc:precip b1.tif

Конвертирование в GeoTIFF с обрезкой по исходным координатам и созданием TFW (world-)файла:

gdal_translate -of GTiff -srcwin 0 0 72 72 -co TFW=YES -b 1 NETCDF:precip.mon.mean.nc:precip b1.tif

Разбиение поднаборов данных ("SUBDATASET") на отдельные файлы netCDF (на выходе — файлы типа "example1", "example2" и т.д.):

gdal_translate -sds example.nc example

HDF4

В HDF4 распространяется множество данных дистанционного зондирования, например MODIS и ASTER.

Импорт данных ASTER L1A:

gdalwarp -overwrite -of GTiff HDF4_EOS:EOS_SWATH:"110601_081441.hdf":VNIR_Band1:ImageData b1.tif 

Импорт данных MODIS Level 3,4:

gdal_translate HDF4_EOS:EOS_GRID:"MCD12Q1.A2001001.h00v08.051.2014287161513.hdf":MOD12Q1:Land_Cover_Type_1 output.tif

Если одна из частей названия SDS (массива данных внутри HDF) имеет пробелы, ее нужно взять в кавычки. Использовать -geoloc для перепроецирования не нужно.

Расчеты

Для разнообразных пересчетов значений пикселей можно использовать gdal_calc. Например, такая команда сбросит в NODATA все пиксели чьё значение больше 16000:

gdal_calc.bat -A input.tif --outfile=output.tif --calc="A*(A<16000)" --NoDataValue=0

Читается такое выражение следующим образом: для каждого пикселя растра input.tif, если его значение меньше 16000 - оставить его таким же, иначе - сбросить в 0.

В выражениях могут использоваться логические "и" и "или", т.е. например в примере ниже: если значение пикселя больше или равно 249, но меньше 255 - сделать его единицей, остальные значения (включая 255) установить в 0.

gdal_calc.bat -A input.tif --outfile=output.tif --calc="1*(logical_and(A>=249,A<255)) " --NoDataValue=0'

Следует обратить внимание, что часть ELSE выражения по умолчанию всегда будет сбрасывать остальные значения в 0. При этом, часть --NoDataValue=0 только говорит к какому значению нужно "приклеить" ярлык NODATA. Таким образом, если в выражении выше сказать, например, --NoDataValue=255, то, несмотря на то, что NODATA установится на 255, значения не удовлетворяющие условию в calc вовсе не станут равны 255, а останутся равны 0. Что бы установить остальные значения в число отличное от 0, нужно сделать это явным образом, например:

gdal_calc --overwrite -A input.tif --outfile=input.tif --calc="1*(A==0) + 255*(A>0)" --NoDataValue=255

Присвоить пикселям со значением 0 значение 1, а остальным присвоить 255, установив на это значение тег NODATA.

Создание растров уменьшенного разрешения (т.н. quicklook, preview)

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

gdal_translate -of "JPEG" -outsize 20% 20%  ALOS_example.tif ALOS_example_preview.jpg -co "WORLDFILE=YES"

Или с помощью gdalwarp, ему можно передать нужный размер пикселя в единицах системы координат, так же можно указать метод, которым будет происходить объединение значений (например усреднение):

gdalwarp -tr 8000 8000 -r average input.tif output.tif

или, если на выходе нужен растр определенного размера:

gdalwarp -ts 8000 8000 -r average input.tif output.tif

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

• Библиотеки GDAL/OGR
• Примеры использования ogr2ogr