project logo

Начало работы с GDAL/OGR

Для начала работы вам не понадобится ничего, кроме терминала (командной строки). Если вы захотите визуализировать свои результаты, вы можете использовать одну из настольных ГИС из состава OSGeo-Live, например, QGIS (Quantum GIS (QGIS)).

Это введение разбито на две части: GDAL (растровые данные) и OGR (векторные данные). Начнём с GDAL.

Ниже описывается, как:

GDAL
  • получать информацию о растровых данных с помощью gdalinfo;
  • конвертировать в различные форматы с помощью gdal_translate;
  • перепроецировать данные с помощью gdalwarp;
  • создавать мозаику из растров с помощью gdalwarp или gdal_merge.py;
  • создавать шейп-файл с тайловым индексом растров с помощью gdaltindex.
OGR
  • получать информацию о векторных данных с помощью ogrinfo;
  • использовать ogr2ogr для конвертации данных в другие форматы.

Знакомимся с GDAL

Вы можете найти тестовые данные здесь /usr/local/share/data. В этом «введении» мы будем использовать набор данных Natural Earth. Лучше работать с копией данных, так что первым делом скопируйте их в ваш домашний каталог.

cd /home/user
cp -R /usr/local/share/data/natural_earth/ ./gdal_natural_earth

Возьмём какой-либо растровый файл из состава Natural Earth и world-файл к нему (файл привязки):

cd /home/user/gdal_natural_earth/HYP_50M_SR_W

Tip

Чтобы просмотреть файл, откройте его в какой-либо настольной ГИС, например, QGIS.

Получение информации о растровом файле с помощью gdalinfo

gdalinfo HYP_50M_SR_W.tif
  Driver: GTiff/GeoTIFF
  Files: HYP_50M_SR_W.tif
         HYP_50M_SR_W.tfw
  Size is 10800, 5400
  Coordinate System is `'
  Origin = (-179.999999999999972,90.000000000000000)
  Pixel Size = (0.033333333333330,-0.033333333333330)
  Metadata:
    TIFFTAG_SOFTWARE=Adobe Photoshop CS3 Macintosh
    TIFFTAG_DATETIME=2009:09:19 10:13:17
    TIFFTAG_XRESOLUTION=342.85699
    TIFFTAG_YRESOLUTION=342.85699
    TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  Image Structure Metadata:
    INTERLEAVE=PIXEL
  Corner Coordinates:
  Upper Left  (-180.0000000,  90.0000000)
  Lower Left  (-180.0000000, -90.0000000)
  Upper Right ( 180.0000000,  90.0000000)
  Lower Right ( 180.0000000, -90.0000000)
  Center      (  -0.0000000,   0.0000000)
  Band 1 Block=10800x1 Type=Byte, ColorInterp=Red
  Band 2 Block=10800x1 Type=Byte, ColorInterp=Green
  Band 3 Block=10800x1 Type=Byte, ColorInterp=Blue
Заметим, что растр имеет следующие характеристики:
  • Драйвер “GTiff/GeoTIFF”;
  • размер 10800x5400;
  • 3 канала типа Byte;
  • есть координатная привязка;
  • система координат не задана.

Простая конвертация форматов

Вначале узнаем, какие драйверы поддерживаются. В этом нам поможет флаг –formats утилит GDAL (например, gdal_translate), с его помощью выводится полный список доступных драйверов.

Для каждого формата сообщается, что он поддерживает:
  • только чтение (ro),
  • чтение/запись (rw),
  • чтение/запись/обновление (rw+).
gdal_translate --formats

Флаг –formats можно использовать для получения детального описания конкретного драйвера с указанием опций создания и разрешённых типов данных.

gdalinfo --format jpeg
gdal_translate --format png

Конвертация

Конвертация осуществляется с помощью утилиты gdal_translate. Выходной формат по умолчанию — GeoTIFF. Флаг -of используется для выбора выходного формата, флаг -co — для указания опций создания выходного файла:

gdal_translate -of JPEG -co QUALITY=40 HYP_50M_SR_W.tif HYP_50M_SR_W.jpg

Флаг -ot служит для изменения типа выходного файла:

gdal_translate -ot Int16 HYP_50M_SR_W.tif HYP_50M_SR_W_Int16.tif

Используйте gdalinfo, чтобы проверить тип данных.

Изменение размера и масштабирование данных

Для изменения размера выходного файла может быть использован флаг -outsize.

gdal_translate -outsize 50% 50% HYP_50M_SR_W.tif  HYP_50M_SR_W_small.tif

Используйте gdalinfo, чтобы проверить размер растра.

Для перемасштабирования данных существует флаг -scale. Доступен также прямой контроль за диапазоном входных и выходных данных. Для вывода минимальных/максимальных значений пикселей растра может быть использован флаг gdalinfo -mm.

Теперь разрежем наш растр на две части с помощью флага -srcwin, который делает копию данных на основе положения пикселов исходного растра (xoff yoff xsize ysize). Вы также можете использовать флаг -projwin, чтобы задать границы растра в координатах географической привязки (ulx uly lrx lry).

gdalinfo -mm HYP_50M_SR_W.tif
gdal_translate -srcwin 0 0 5400 5400 HYP_50M_SR_W.tif  west.tif
gdal_translate -srcwin 5400 0 5400 5400 HYP_50M_SR_W.tif  east.tif

Индекс растровых тайлов с помощью gdaltindex

Вы можете создать шейп-файл с индексом растровых тайлов. Для каждого растра будет сгенерирован полигон с границами по охвату растра и с путём к файлу в атрибутах.

gdaltindex index_natural_earth.shp *st.tif

Посмотрим на получившийся шейп-файл в QGIS и ogrinfo (мы ещё рассмотрим ogrinfo ниже).

../../_images/gdal_gdaltindex10.png 
ogrinfo ../HYP_50M_SR_W/ index
INFO: Open of `../HYP_50M_SR_W/'
    using driver `ESRI Shapefile' successful.

Layer name: index
Geometry: Polygon
Feature Count: 2
Extent: (-180.000000, -90.000000) - (180.000000, 90.000000)
Layer SRS WKT: (unknown)
location: String (255.0)
OGRFeature(index):0
  location (String) = east.tif
  POLYGON ((-0.00000000001796 90.0,179.999999999964047 90.0,179.999999999964047 -89.999999999982009,-0.00000000001796 -89.999999999982009,-0.00000000001796 90.0))

OGRFeature(index):1
  location (String) = west.tif
  POLYGON ((-179.999999999999972 90.0,-0.00000000001796 90.0,-0.00000000001796 -89.999999999982009,-179.999999999999972 -89.999999999982009,-179.999999999999972 90.0))

Перепроецирование

Для следующих действий мы предполагаем, что растр HYP_50M_SR_W.tif имеет нужный охват. Как мы выяснили ранее с помощью gdalinfo, у растра не задана система координат, поэтому первым делом мы назначим WGS84 в качестве таковой.

gdal_translate -a_srs WGS84 HYP_50M_SR_W.tif HYP_50M_SR_W_4326.tif

Команда gdalwarp служит для перепроецирования растров. Попробуем перепроецировать наш растр в проекцию Меркатора:

gdalwarp -t_srs '+proj=merc +datum=WGS84' HYP_50M_SR_W_4326.tif mercator.tif

Используйте gdalinfo, чтобы проверить изменения и посмотреть на свойства растра.

../../_images/gdal_mercator10.png

Теперь перепроецируем растр в ортографическую проекцию:

gdalwarp -t_srs '+proj=ortho +datum=WGS84' HYP_50M_SR_W_4326.tif ortho.tif
../../_images/gdal_ortho10.png

Вы обратили внимание, что земные полюса “обрезаны”? gdalwarp не читает данные полностью, поскольку границы возле полюсов нельзя перепроецировать. В качестве одного из возможных решений можно заставить gdalwarp читать дополнительный объём данных вокруг фрагментов растра. Подробнее читайте на странице RasterTutorial http://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial.

Создание мозаик

gdal_merge.py — скрипт на Python, который можно использовать для простых задач создания мозаик. Например, создадим мозаику из двух растров (east.tif и west.tif) в виде единого файла merged.tif:

gdal_merge.py  east.tif west.tif -o merged.tif

Подобная задача может быть решена и с помощью gdalwarp, это утилита имеет ряд преимуществ перед gdal_merge.py, но может медленно работать при сшивке большого количества растров:

gdalwarp east.tif west.tif warpmerged.tif

Знакомимся с OGR

cd /home/user/gdal_natural_earth/

Tip

Чтобы просмотреть данные, откройте шейп-файл в любой настольной ГИС типа QGIS.

Получение информации о векторных данных с помощью ogrinfo

ogrinfo -ro /home/user/gdal_natural_earth
INFO: Open of `/home/user/gdal_natural_earth'
      using driver `ESRI Shapefile' successful.
  1: ne_10m_populated_places (3D Point)
  2: ne_10m_geography_regions_polys (3D Polygon)
  3: ne_10m_admin_1_states_provinces_shp (3D Polygon)
  4: ne_10m_urban_areas (3D Polygon)
  5: ne_10m_geography_marine_polys (3D Polygon)
  6: ne_10m_land (3D Polygon)
  7: ne_10m_geography_regions_elevation_points (3D Point)
  8: ne_10m_admin_0_countries (3D Polygon)
  9: ne_10m_rivers_lake_centerlines (3D Line String)
  10: ne_10m_lakes (3D Polygon)
  11: ne_10m_geography_regions_points (3D Point)
  12: ne_10m_ocean (3D Polygon)

Краткую информацию о векторных данных можно получить с помощью утилиты ogrinfo с флагом -so.

ogrinfo -ro -so ne_10m_admin_0_countries.shp ne_10m_admin_0_countries
INFO: Open of `ne_10m_admin_0_countries.shp'
      using driver `ESRI Shapefile' successful.

Layer name: ne_10m_admin_0_countries
Geometry: 3D Polygon
Feature Count: 254
Extent: (-180.000000, -90.000000) - (180.000000, 83.634101)
Layer SRS WKT:
GEOGCS["GCS_WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS_84",6378137.0,298.257223563]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]]
scalerank: Integer (4.0)
featurecla: String (30.0)
labelrank: Real (16.6)
sovereignt: String (254.0)
sov_a3: String (254.0)
adm0_dif: Real (16.6)
level: Real (16.6)
type: String (254.0)
admin: String (254.0)
adm0_a3: String (254.0)
geou_dif: Real (16.6)
geounit: String (254.0)
gu_a3: String (254.0)
su_dif: Real (16.6)
subunit: String (254.0)
su_a3: String (254.0)
brk_diff: Real (16.6)
name: String (254.0)
name_long: String (254.0)
brk_a3: String (254.0)
brk_name: String (254.0)
brk_group: String (254.0)
abbrev: String (254.0)
postal: String (254.0)
formal_en: String (254.0)
formal_fr: String (254.0)
note_adm0: String (254.0)
note_brk: String (254.0)
name_sort: String (254.0)
name_alt: String (254.0)
mapcolor7: Real (16.6)
mapcolor8: Real (16.6)
mapcolor9: Real (16.6)
mapcolor13: Real (16.6)
pop_est: Real (16.6)
gdp_md_est: Real (16.6)
pop_year: Real (16.6)
lastcensus: Real (16.6)
gdp_year: Real (16.6)
economy: String (254.0)
income_grp: String (254.0)
wikipedia: Real (16.6)
fips_10: String (254.0)
iso_a2: String (254.0)
iso_a3: String (254.0)
iso_n3: String (254.0)
un_a3: String (254.0)
wb_a2: String (254.0)
wb_a3: String (254.0)
woe_id: Real (16.6)
adm0_a3_is: String (254.0)
adm0_a3_us: String (254.0)
adm0_a3_un: Real (16.6)
adm0_a3_wb: Real (16.6)
continent: String (254.0)
region_un: String (254.0)
subregion: String (254.0)
region_wb: String (254.0)
name_len: Real (16.6)
long_len: Real (16.6)
abbrev_len: Real (16.6)
tiny: Real (16.6)
homepart: Real (16.6)

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

ogrinfo -ro ne_10m_admin_0_countries.shp ne_10m_admin_0_countries

Вы можете отфильтровать вывод ogrinfo с помощью стандартной утилиты grep и получить, например, только атрибуты поля COUNTRY.

ogrinfo ne_10m_admin_0_countries.shp ne_10m_admin_0_countries | grep 'admin'

  admin (String) = Aruba
  admin (String) = Afghanistan
  admin (String) = Angola
  admin (String) = Anguilla
  admin (String) = Albania
  admin (String) = Aland
  admin (String) = Andorra
etc.

Вы можете сконвертировать ваши данные в другие форматы. Список поддерживаемых форматов можно получить с помощью флага –formats.

Использование ogr2ogr для конвертации данных между форматами

Вы можете использовать утилиту ogr2ogr для конвертации векторных данных (стандарта simple features) между различными форматами. Полный список форматов OGR с указанием поддержки чтения/записи выводится флагом –formats.

Давайте сконвертируем шейп-файл countries в формат GML.

ogr2ogr --formats
ogr2ogr -f GML countries.xml ne_10m_admin_0_countries.shp

Что ещё попробовать

Есть несколько действий, которые стоит попробовать при работе с GDAL/OGR:

  1. Попробуйте создать растровую мозаику с помощью gdalwarp или gdal_merge.py
  2. Попробуйте создать внутренние слои “пирамид” для оптимизации просмотра растра в низком разрешении с помощью gdaladdo
  3. QGIS использует GDAL/OGR для поддержки большого числа форматов. Эта ГИС также предоставляет плагин GdalTools для работы с растровыми данными, который интегрирует утилиты GDAL в QGIS.
  4. Попробуйте ogr2ogr для импорта/экспорта векторных данных в другие форматы (например, PostGIS). Эта утилита имеет довольно длинный список опций.
  5. Попробуйте конвертацию данных в QGIS через OGR.

Что дальше?

Это «введение» — только первый шаг на дороге освоения GDAL/OGR. На самом деле, доступно гораздо больше функциональности, чем описано здесь.

Официальная страница GDAL:

Всё об OGR:

Руководство по GDAL: