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

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

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

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

Содержание 1 Общие сведения 2 Конвертация 3 Стэки, композиты, мозаики 4 Работа с NODATA 5 Работа с системами координат 5.1 Переназначение системы координат 5.2 Конвертация с перепроецированием 6 Обрезка 7 Работа с рельефом 8 Работа с цветовой таблицей 8.1 С использованием VRT 9 Построение изолиний 10 netCDF 11 HDF4 12 Расчеты 13 Создание растров уменьшенного разрешения (т.н. quicklook, preview) 14 Ссылки по теме

[править] Общие сведения

Утилиты 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

[править] Ссылки по теме

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