Ортокоррекция данных OrbView-3 с помощью GDAL

Эта страница является черновиком статьи.

Снимки OrbView-3 в исходном виде не привязаны и не орторектифицированы. Чтобы привести их в состояние, пригодное для работы, нужно выполнить географическую привязку данных и орторектификацию.

Краткое руководство

Три простых шага для получения компактного восьмибитного растра из исходного:

1. Скачиваем нужные участки рельефа SRTM в GeoTIFF, для орторектификации (распрямления) по ним снимков. Для этого:

• на странице http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp смотрим на картинку и выписываем куда-нибудь номера рядов и столбцов с интересующей нас местностью;
• скачиваем с http://gis-lab.info/data/srtm-tif выбранные zip-файлы, складываем в один каталог;
• распаковываем:

 for i in srtm*zip; do yes|unzip $i; done

2. Для удобства объединяем все файлы SRTM в формате GeoTIFF в единый виртуальный растр:

 gdalbuildvrt srtm.vrt srtm*tif

3. Для получения компактного восьмибитного геотиффа из исходного сырого делаем следующее:

 gdalwarp -co TILED=YES -co JPEG_QUALITY=85 -co COMPRESS=JPEG -multi -dstnodata 0 -srcnodata 0 -overwrite -ot byte -t_srs epsg:3857
 -rpc -to 'RPC_DEM=/home/kom/Downloads/srtm.vrt' 3v050901p0000883761a520004300642m_001631911.tif 3v050901p0000883761a520004300642m_001631911.rec.tif 

TODO:
Bishop проверить с параметром -wo "INIT_DEST=NO_DATA" или -wo "INIT_DEST=0" (для случая -dstnodata 0)
подробнее о параметрах

• Узнать, как в том же шаге сделать internal nodata mask, для более аккуратных краёв изображения;

Развёрнутое руководство

Здесь рассказывается про разные варианты преобразований. Если вы видите, что автор пропустил какие-то очевидные для специалиста вещи — пожалуйста, исправьте.

Использованный пример

Здесь используются иллюстрации из сцены Orbview-3 на территорию Республики Беларусь. Вы можете использовать этот же пример для своих экспериментов.

• Оригинальный снимок 3v050909p0000897861a520004700712m_001631680, скачанный через EarthSat (скачать)
• Фрагмент ASTER GDEM использованный для орторектификации (скачать)
• Треки (GPX) для проверки результата (скачать)
• Результат орторектификации в ENVI EX (скачать)

Процесс

Для тех, кто не хочет читать длинные рассуждения с примерами — краткие выводы:

1. gdalwarp без DEM — очень грубая привязка.
2. gdalwarp с DEM — получше, но всё равно неточно.
3. ENVI EX с DEM — практически идеально, только дорого.

Bishop Проверено опытным путем. Автор при ортокоррекции с использованием ENVI EX применял поправку геоида (22.016 м)
Bishop В GDAL чтобы получить аналогичные результаты нужно добавить:

 gdalwarp -co TILED=YES -co JPEG_QUALITY=85 -co COMPRESS=JPEG -multi -dstnodata 0 -srcnodata 0 -overwrite -ot byte -t_srs epsg:3857
 -rpc -to 'RPC_DEM=/home/kom/Downloads/srtm.vrt' -to 'RPC_HEIGHT=22.016' 3v050901p0000883761a520004300642m_001631911.tif 3v050901p0000883761a520004300642m_001631911.rec.tif

Bishop Расчет может быть сделан онлайн на этом ресурсе. Я взял среднюю точку снимка.
Bishop Кроме того, для устранения потери разрешения необходимо в исходном файле добавить контрольные точки в WGS84/UTM соответствующей зоны (в нашем случае 35N). Для этого в файле номер_сцены.pvl ищем координаты углов в СК WGS84, пересчитываем их в WGS84/UTM и добавляем их к файлу. Ну и выходной файл создавать в той же СК WGS84/UTM.

Сравнение результатов ENVI EX + ASTER GDEM, GDAL + ASTER GDEM и GDAL + SRTM

Для начала, заглянем в содержимое ZIP-архива с одной сценой, скачанного с сайта USGS:

$ ls -1 3v050909p0000897861a520004700712m_001631680*
3v050909p0000897861a520004700712m_001631680_aoi.dbf
3v050909p0000897861a520004700712m_001631680_aoi.prj
3v050909p0000897861a520004700712m_001631680_aoi.shp
3v050909p0000897861a520004700712m_001631680_aoi.shx
3v050909p0000897861a520004700712m_001631680.att
3v050909p0000897861a520004700712m_001631680.dbf
3v050909p0000897861a520004700712m_001631680.eph
3v050909p0000897861a520004700712m_001631680.jgw
3v050909p0000897861a520004700712m_001631680.jpg
3v050909p0000897861a520004700712m_001631680.prj
3v050909p0000897861a520004700712m_001631680.pvl
3v050909p0000897861a520004700712m_001631680_rpc.txt
3v050909p0000897861a520004700712m_001631680.shp
3v050909p0000897861a520004700712m_001631680.shx
3v050909p0000897861a520004700712m_001631680_src.dbf
3v050909p0000897861a520004700712m_001631680_src.prj
3v050909p0000897861a520004700712m_001631680_src.shp
3v050909p0000897861a520004700712m_001631680_src.shx
3v050909p0000897861a520004700712m_001631680.tif

Имеем: шейп-файлы с покрытием, jpg с "превьюшкой" (привязанный!), собственно TIFF с данными, файл с параметрами RPC преобразований (scene_rpc.txt), файл scene.pvl, содержащий некое описание данных и их параметров.

Если загрузить, например, в qgis, TIFF-файл с данными, то никакой привязкой там пахнуть и не будет. Попробуем разобраться:

$ gdalinfo 3v050909p0000897861a520004700712m_001631680.tif 
Warning 1: TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered
Driver: GTiff/GeoTIFF
Files: 3v050909p0000897861a520004700712m_001631680.tif
       3v050909p0000897861a520004700712m_001631680_rpc.txt
Size is 8016, 25600
Coordinate System is `'
Metadata:
  TIFFTAG_MINSAMPLEVALUE=0
  TIFFTAG_MAXSAMPLEVALUE=2047
Image Structure Metadata:
  INTERLEAVE=BAND
RPC Metadata:
  LINE_OFF= +012800.00 pixels
  SAMP_OFF= +004008.00 pixels
  LAT_OFF= +55.02030000 degrees
  LONG_OFF= +027.04780000 degrees
  HEIGHT_OFF= +0179.000 meters
  LINE_SCALE= +012800.00 pixels
  SAMP_SCALE= +004008.00 pixels
  LAT_SCALE= +00.12380000 degrees
  LONG_SCALE= +000.06850000 degrees
  HEIGHT_SCALE= +0300.000 meters
  LINE_NUM_COEFF= -2.104832000000000E-03  -1.642616000000000E-02  -1.027459000000000E+00  +4.182002500000000E-03  -1.902795200000000E-03  +1.614313300000000E-05  +4.786355800000000E-04  -2.127866900000000E-04  +6.958830700000000E-03  -2.260572200000000E-06  -2.225955200000000E-07  -3.746937200000000E-07  +4.648645700000000E-04  -1.801288800000000E-08  +5.140758300000000E-06  +7.566147900000000E-04  -5.452440900000000E-07  +1.394079900000000E-07  -1.828159600000000E-05  +2.421558100000000E-09 
  LINE_DEN_COEFF= +1.000000000000000E+00  -5.006651300000000E-04  -1.457830900000000E-03  +6.037474400000000E-04  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00 
  SAMP_NUM_COEFF= +3.145351000000000E-04  +1.023427000000000E+00  -3.439455200000000E-03  +1.730013100000000E-02  +5.102439600000000E-03  +1.245288300000000E-03  -1.578704500000000E-03  -2.939111200000000E-03  -2.317010900000000E-04  +2.847267800000000E-05  +2.422610700000000E-05  +6.633964600000000E-06  -1.694999000000000E-03  +6.913555700000000E-07  +1.531221300000000E-04  +1.486522300000000E-05  +1.555096200000000E-07  -5.617327200000000E-06  -2.950330800000000E-05  +1.262439900000000E-08 
  SAMP_DEN_COEFF= +1.000000000000000E+00  -6.035266200000000E-04  +6.161647000000000E-03  +6.386015900000000E-04  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00  +0.000000000000000E+00 
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,25600.0)
Upper Right ( 8016.0,    0.0)
Lower Right ( 8016.0,25600.0)
Center      ( 4008.0,12800.0)
Band 1 Block=8016x1 Type=UInt16, ColorInterp=Gray

Как и ожидалось, в файле не содержится данных о привязке. Зато GDAL прочитал данные RPC (rational polynomial coefficients), нужные для корректной привязки и трансформации. Подробнее про орторектификацию с использованием RPC можно прочитать по ссылкам из темы форума: http://gis-lab.info/forum/viewtopic.php?p=27889#p27889. Если команда gdalinfo у вас не вывела метаданные RPC, проверьте версию --- нужен GDAL не менее 1.8.1.

Для корректной орторектификации требуются параметры RPC и информация о высоте над геоидом. Что-то там про высоту в выводе gdalinfo есть (интересующиеся могут разобраться с этим самостоятельно). Попробуем в лоб:

$ gdalwarp -rpc 3v050909p0000897861a520004700712m_001631680.tif test.tif
Warning 1: TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered
Creating output file that is 12925P x 23537L.
Warning 1: TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered
Processing input file 3v050909p0000897861a520004700712m_001631680.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.

Bishop Так делать не рекомендую - есть подозрение, что точки привязки неверно берутся - лучше самостоятельно прописать в СК WGS84/UTM.

Теперь посмотрим еще раз, прописалась ли информация в результирующий файл:

$ gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
       test_rpc.txt
Size is 12925, 23537
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (26.981501010426538,55.143013345911761)
Pixel Size = (0.000010399352347,-0.000010399352347)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  26.9815010,  55.1430133) ( 26d58'53.40"E, 55d 8'34.85"N)
Lower Left  (  26.9815010,  54.8982438) ( 26d58'53.40"E, 54d53'53.68"N)
Upper Right (  27.1159126,  55.1430133) ( 27d 6'57.29"E, 55d 8'34.85"N)
Lower Right (  27.1159126,  54.8982438) ( 27d 6'57.29"E, 54d53'53.68"N)
Center      (  27.0487068,  55.0206286) ( 27d 2'55.34"E, 55d 1'14.26"N)
Band 1 Block=12925x1 Type=UInt16, ColorInterp=Gray

Ага, появились вполне правдоподобные данные привязки (в географической СК, правда). Загружаем полученный файл в QGIS... Гм, в выбранный район, конечно, попали. Но привязка получилась плюс-минус лапоть. Причём лапоть семидесятиметровый.

Гугл рассказал, что для точной привязки нужно использовать данные рельефа DEM (digital elevation model). ОК. Идём на сайт Aster GDEM (http://www.gdem.aster.ersdac.or.jp/index.jsp), регистрируемся, заходим в раздел Search. Выбираем "Select tiles by shapefile", скармливаем файл покрытия scene.shp, скачиваем найденное, распаковываем. Получаем несколько (в данном случае — 4) файлов с данными DEM с именами вида ASTGTM2_N55E026_dem.tif (один тайл соответствует "квадрату" градусной сетки). К сожалению, ни GDAL, ни ENVI EX не умеет работать с несколькими файлами DEM для одной сцены. Поэтому склеиваем их с помощью gdal_merge.py в один GeoTIFF:

$ gdal_merge.py -o DEM_merged.tif  ASTGTM2_N5[45]E02[67]/*_dem.tif 
0...10...20...30...40...50...60...70...80...90...100 - done.

Документация на gdalwarp говорит, что для орторектификации можно указать DEM-файл с помощью параметра -to 'RPC_DEM=DEM_merged.tif'.
Bishop И не забываем про геоид -to 'RPC_HEIGHT=22.016' Пробуем:

$ gdalwarp -rpc -to 'RPC_DEM=/home/ftp/maps/GDEM/DEM_merged.tif'  3v050909p0000897861a520004700712m_001631680.tif test.tif
Warning 1: TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered
Processing input file 3v050909p0000897861a520004700712m_001631680.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.

Загружаем в QGIS — совсем другой коленкор! Хотя всё равно получили погрешность порядка 10 метров. Причины слёту не выяснены.

Для трансформации в другую систему координат просто добавьте -t_srs proj_definition. Если при этом gdalwarp отказывается трансформировать что-либо, скачайте соседние тайлы GDEM — скорее всего, "поля" в трансформированном растре "выскакивают" к соседям.

Зато, если скормить исходный растр и DEM программе (увы, платной и под Windows) ENVI EX, она после длительного пережёвывания на выходе даст практически идеально привязанный файл в проекции UTM нужной зоны. Пример результата:

Если вы знаете, как достичь таких же результатов с помощью GDAL — поделитесь.
Bishop Поделился