Примеры использования инструментов GDAL
по адресу http://gis-lab.info/qa/gdal-examples.html
Перечь примеров для справки
GDAL/OGR - библиотека для работы с географическими форматами данных. GDAL представляет собой набор утилит для обработки растровых данных, в то время, как OGR предназначена для работы с векторными форматами. В статье рассматриваются некоторые практические примеры применения утилит этой библиотеки для работы с растровыми данными.
С библиотекой GDAL так же поставляется утилита ogr2ogr, предназначенная для работы с векторными данными. Примеры использования этой утилиты приводятся в другой статье.
Если у вас есть свои часто используемые примеры - присылайте автору или дописывайте прямо здесь.
Общие сведения
Утилиты GDAL предназначены для конвертации растровых данных из одного формата в другой и выполнения над ними различных операций. Установить утилиты GDAL для Windows можно с помощью OSGeo4W. Удобный визуальный интерфейс для утилит имеется в GDALTools для QGIS. В комплект GDAL входят следующие утилиты:
- gdal_calc - растровый калькулятор, арифметические операции с растрами;
- gdal-config - получить опции необходимые для создания ПО использующего GDAL;
- gdal_contour - получение изолиний по цифровым моделям рельефа (ЦМР);
- gdal_calc - растровый калькулятор, арифметические операции с растрами;
- gdal_edit.py - редактировать информацию о растре;
- gdal_fillnodata.py - заполнение областей имеющих значение NODATA;
- gdal_merge.py - создание мозаик и композитных изображений;
- gdal_polygonize.py - векторизовать растр с получением полигонального слоя;
- gdal_proximity.py - расчитать растр близости;
- gdal_rasterize - растеризация векторных данных;
- gdal_retile.py - создать новый набор тайлов и/или перестроить пирамидные слои;
- 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.tif output.tif
Пакетная конвертация всех JPG в TIF (Windows):
for %i in (*.jpg) do gdal_translate %i %~ni.tif
Стэки, композиты, мозаики
Создание композитного изображения из серии отдельных растров, каждый из которых в своем файле TIF. Разрешение выходного файла устанавливается по первому их растров. Таким образом, если первый канал 15 м, а остальные 30 м, то последние будут пересчитаны на 15 м. Чтобы указать, что производится помещение каждого растра в свой слой, а не мозаицирование, используется ключ -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
Работа с системами координат
Переназначение системы координат
Если вам нужно просто перепрописать систему координат, без пересчета самого растра:
gdal_edit -a_srs "EPSG:4326" input.tif
Конвертация с перепроецированием
gdalwarp позволяет не только конвертировать данные из одного формата в другой, но и одновременно произвести перепроецирование данных из одной системы координат в другую. Для этого используются параметры:
- -a_srs используется для указания системы координат для данных
- -s_srs используется для перезаписи информации о системе координат
- -t_srs перепроецирования данных в требуемую систему координат
Например, перепроецировать из проекции растра в проекцию Альберса можно так:
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.
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. Использовать -geoloc для перепроецирования не нужно.
Импорт данных ASTER L1A:
gdalwarp -overwrite -of GTiff HDF4_EOS:EOS_SWATH:"110601_081441.hdf":VNIR_Band1:ImageData b1.tif
Расчеты
Для разнообразных пересчетов значений пикселей можно использовать 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'
Создание растров уменьшенного разрешения (т.н. 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