Ортокоррекция данных OrbView-3 с помощью GDAL: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
(Исправление про артефакты)
 
(не показано 30 промежуточных версий 6 участников)
Строка 1: Строка 1:
Снимки [[Каталог_данных_Orbview3|Orbview-3]] в исходном виде не привязаны и не орторектифицированы. Чтобы привести их в состояние, пригодное для работы, нужно выполнить географическую привязку данных и орторектификацию.
{{Статья|Опубликована|orbview3-ortho-gdal}}
{{Аннотация|Как произвести ортокоррекцию данных OrbView-3 с помощью GDAL}}


=== Краткое руководство ===
== Введение ==


Три простых шага для получения компактного восьмибитного растра из исходного:
В данной статье рассматривается способ ортокоррекции космических снимков с КА OrbView-3 при помощи библиотеки [http://gdal.org GDAL].


1. Скачиваем нужные участки рельефа SRTM в GeoTIFF, для орторектификации (распрямления) по ним снимков. Для этого:
Космические снимки с КА OrbView-3 свободно распространяются и могут быть получены через [[Каталог данных OrbView-3|каталог данных OrbView-3]].
* на странице http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp смотрим на картинку и выписываем куда-нибудь номера рядов и столбцов с интересующей нас местностью;
Однако, в исходном виде, снимки не имеют географической привязки - в формат распространения (TIFF) не внедрены специальные теги для описания привязки, т.е. формат не представляет из себя классического GeoTIFF. Однако, в комплекте поставки имеются все необходимые метаданные для осуществления привязки и ортокоррекции снимков. Как это осуществить будет рассказано далее.
* скачиваем с 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
 
Для выполнения ортокоррекции космического снимка необходимо:
# Свободно распространяемое программное обеспечение [http://gdal.org/gdalwarp.html gdalwarp], [http://gdal.org/gdalbuildvrt.html gdalbuildvrt], [http://gdal.org/gdal_translate.html gdal_translate] из состава библиотеки GDAL. Для получения программного обеспечения можно перейти по [http://trac.osgeo.org/gdal/wiki/DownloadingGdalBinaries ссылке], где перечислены варианты загрузки.
# Космический снимок с КА OrbView-3. Снимки на необходимую территорию ищем в [[Каталог данных OrbView-3|каталоге]].
# Рельеф на территорию съемки в формате в виде ЦМР (Цифровой модели рельефа)


3. Для получения компактного восьмибитного геотиффа из исходного сырого делаем следующее:
В качестве рельефа могут быть использованы: [http://gis-lab.info/qa/srtm.html SRTM], [http://gis-lab.info/qa/aster-gdem.html ASTER GDEM] v.1 и v.2 и др.
Рельеф загружаем отсюда:
* [http://gis-lab.info/data/srtm-tif SRTM]
* [http://www.gdem.aster.ersdac.or.jp/index.jsp ASTER GDEM]


  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
Для демонстрации операций по ортокоррекции можно воспользоваться подготовленными наборами данных (снимки и рельеф) по районам [http://gis-lab.info/data/orbview-3/samples/kursk/ Курской области] и [http://gis-lab.info/data/orbview-3/samples/rb/ Республике Беларусь].


TODO:
Набор данных по Курской области включает в себя:
* Узнать, как в том же шаге сделать internal nodata mask, для более аккуратных краёв изображения;
3V050401P0000692211A520018202082M_001639429.ZIP - космический снимок с КА OrbView-3 (архив)
3v050401p0000692211a520018202082m_001639429_ortho.img - результат ортокоррекции снимка в формате ERDAS Imagine, выполненный в GDAL
3v050401p0000692211a520018202082m_001639429_ortho.img.aux.xml - " -
3v050401p0000692211a520018202082m_001639429_ortho.rrd - " -
dem.tfw - файл рельефа на территорию съемки
dem.tif - " -
dem.tif.aux.xml - " -
dem.tif.ovr - " -
dem.tif.xml - " -
kursk_ortho.7z  - результат ортокоррекции снимка в формате TIFF, выполненный в GDAL (архив)
kursk_ortho_envi_orbview.7z  - результат ортокоррекции снимка в формате TIFF, выполненный в ENVI (архив)


=== Развёрнутое описание ===
Набор данных по Республике Беларусь включает в себя:
3v050909p0000897861a520004700712m_001631680.zip - космический снимок с КА OrbView-3 (архив)
aster_gdem.tif - файл рельефа на территорию съемки
ov3-check.7z - наземные точки снятые с помощью GPS
result_envi.tif  - результат ортокоррекции снимка, выполненный в ENVI EX
result_gdal.7z  - результат ортокоррекции снимка, выполненный в GDAL (архив)
result_gdal_geoid_corrected.7z - результат ортокоррекции снимка с использованием коррекции геоида, выполненный в GDAL (архив)
tracks.7z - наземные точки снятые с помощью GPS


На этой странице рассказано про варианты преобразований.
== Подготовка к ортокорреции ==


Если вы видите, что автор — пропустил какие-то очевидные для специалиста вещи (а автор на самом деле дилетант) — пожалуйста, исправьте.
В данной статье все примеры будут демонстрироваться на наборе данных по Республике Беларусь.


Для тех, кто не хочет читать длинные рассуждения с примерами — краткие выводы:
Для начала рассмотрим содержимое ZIP-архива с материалами космической съемки с КА OrbView-3 (3V050909P0000897861A520004700712M_001631680.ZIP). Для этого выведем листинг директории командой:
#gdalwarp без DEM — очень грубая привязка.
#gdalwarp с DEM — получше, но всё равно неточно.
#ENVI EX с DEM — практически идеально, только дорого.


Для начала, заглянем в содержимое ZIP-архива с одной сценой, скачанного с сайта USGS:
UNIX
$ ls -1 3V050909P0000897861A520004700712M_001631680*
Windows
> dir * /B
<pre>
<pre>
$ ls -1 3v050909p0000897861a520004700712m_001631680*
3v050909p0000897861a520004700712m_001631680_aoi.dbf
3v050909p0000897861a520004700712m_001631680_aoi.dbf
3v050909p0000897861a520004700712m_001631680_aoi.prj
3v050909p0000897861a520004700712m_001631680_aoi.prj
Строка 56: Строка 77:
</pre>
</pre>


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


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


Как и ожидалось, в файле не содержится данных о привязке. Зато GDAL прочитал данные RPC (rational polynomial coefficients), нужные для корректной привязки и трансформации. Подробнее про орторектификацию с использованием RPC можно прочитать по ссылкам из темы форума: http://gis-lab.info/forum/viewtopic.php?p=27889#p27889. Если команда gdalinfo у вас не вывела метаданные RPC, проверьте версию --- нужен GDAL не менее 1.8.1.
Как и ожидалось, в файле не содержится данных о привязке. Зато GDAL прочитал данные RPC (rational polynomial coefficients), нужные для корректной привязки и трансформации. Подробнее про ортокоррекцию с использованием RPC можно прочитать в статье "[[Ортокоррекция космических снимков с использованием RPC]]", а также в [http://gis-lab.info/forum/viewtopic.php?p=27889#p27889 соответствующей теме форума]. Если команда gdalinfo у вас не вывела метаданные RPC, проверьте версию - нужен GDAL не менее 1.8.1.
 
Прежде чем формировать команду для ортокоррекции необходимо определиться с системой координат результирующего изображения. Систему координат нужно подобрать таким образом, чтобы пиксел изображения меньше всего искажался. Для выбора системы координат, исходное изображение было привязано по углам с использованием координат из файла 3v050909p0000897861a520004700712m_001631680.pvl. После этого файл теперь уже в формате GeoTIFF был загружен в QGIS и на нем были испробованы различные системы координат с преобразованием "на лету". При этом, обращалось внимание, чтобы пиксел изображения оставался "квадратным" и вращался на наименьший угол. В результате экспериментов было принято решение об использовании системы координат WGS84/UTM соответствующей зоны.
 
Для нашего случая система координат будет '''WGS 84 / UTM zone 35N''' или '''epsg:32635'''.
 
Получение необходимых данных о рельефе рассмотрим на примере работы с данными '''SRTM'''. Для получения файла рельефа необходимо:
<ol>
<li> На странице [http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp SRTM Data Selection Options] выписываем номера рядов и столбцов с территорией съемки.</li>
<li> Скачиваем с http://gis-lab.info/data/srtm-tif выбранные zip-файлы и складываем в один каталог.</li>
<li> Распаковываем (команда для операционной системы UNIX)</li>
for i in srtm*zip; do yes|unzip $i; done
 
<li>Объединяем все файлы SRTM в формате GeoTIFF в единый виртуальный растр</li>


Для корректной орторектификации требуются параметры RPC и информация о высоте над геоидом. Что-то там про высоту в выводе gdalinfo есть (интересующиеся могут разобраться с этим самостоятельно). Попробуем в лоб:
Команда для операционной системы UNIX:
  gdalbuildvrt srtm.vrt srtm*tif
Команда для операционной системы Windows:
  gdalbuildvrt.exe srtm.vrt srtm*tif
</ol>
 
Для получения файла рельефа '''ASTER GDEM''' переходим на [http://www.gdem.aster.ersdac.or.jp/index.jsp сайт Aster GDEM], регистрируемся и заходим в раздел Search. Выбираем "Select tiles by shapefile", загружаем файл покрытия scene.shp, скачиваем и распаковываем. Получаем несколько (в данном случае — 4) файлов с данными DEM с именами вида ASTGTM2_N55E026_dem.tif (один тайл соответствует "квадрату" градусной сетки). Склеиваем все файлы в один GeoTIFF с помощью gdal_merge.py:
<pre>
$ 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.
</pre>
 
== Ортокорреция ==
Для выполнения ортокоррекции космического снимка с КА OrbView-3 необходимо выполнить следующую команду.
 
  > gdalwarp -dstnodata 0 -srcnodata 0 -overwrite -t_srs epsg:32635 -wo "INIT_DEST=NO_DATA" -rpc -to "RPC_DEM=D:\temp\rb\DEM_merged.tif"
  D:\temp\rb\3V050909P0000897861A520004700712M_001631680\3v050909p0000897861a520004700712m_001631680.tif D:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec.tif
 
Параметры команды gdalwarp можно посмотреть по следующей [http://gdal.org/gdalwarp.html ссылке], подробнее о параметрах -wo можно почитать [http://www.gdal.org/structGDALWarpOptions.html#0ed77f9917bb96c7a9aabd73d4d06e08 здесь].
 
Следует отметить, что в зависимости от версии библиотеки GDAL и операционной системы команда может отработать или завершиться с ошибкой. В случае если команды выполнилась с ошибкой, следует прописать систему координат и экстент исходного изображения.
 
Вариант №1. Попробовать извлечь из метаданных. Минус варианта - выходная система координат будет WGS84, что приводит к потере разрешения (размера пиксела на местности).
<pre>
<pre>
$ gdalwarp -rpc 3v050909p0000897861a520004700712m_001631680.tif test.tif
$ gdalwarp -rpc 3v050909p0000897861a520004700712m_001631680.tif test.tif
Строка 106: Строка 162:
Warning 1: TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered
Warning 1: TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered
Processing input file 3v050909p0000897861a520004700712m_001631680.tif.
Processing input file 3v050909p0000897861a520004700712m_001631680.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.</pre>
$ gdalinfo test.tif
 
Проверяем, прописалась ли информация в результирующий файл:
 
<pre>$ gdalinfo test.tif
Driver: GTiff/GeoTIFF
Driver: GTiff/GeoTIFF
Files: test.tif
Files: test.tif
Строка 133: Строка 192:
Lower Right (  27.1159126,  54.8982438) ( 27d 6'57.29"E, 54d53'53.68"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)
Center      (  27.0487068,  55.0206286) ( 27d 2'55.34"E, 55d 1'14.26"N)
Band 1 Block=12925x1 Type=UInt16, ColorInterp=Gray
Band 1 Block=12925x1 Type=UInt16, ColorInterp=Gray</pre>
</pre>


Ага, появились вполне правдоподобные данные привязки (в географической СК, правда). Загружаем полученный файл в QGIS... Гм, в выбранный район, конечно, попали. Но привязка получилась плюс-минус лапоть. Причём лапоть семидесятиметровый.
Вариант №2. Прописать систему координат и GCP углов изображения. Для этого в файле scene.pvl ищем координаты углов в СК WGS84, пересчитываем их в WGS84/UTM 35N и добавляем их к файлу.
 
Гугл рассказал, что для точной привязки нужно использовать данные рельефа 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:
<pre>
<pre>
$ gdal_merge.py -o DEM_merged.tif ASTGTM2_N5[45]E02[67]/*_dem.tif  
> gdal_translate -a_srs epsg:32635 -gcp 0.5 0.5 498793 6.11075e+006 173.385 -gcp 8015.5 0.5 507345 6.11027e+006 187.386 -gcp 8015.5 25599.5 507347 6.08351e+006 204.881
-gcp 0.5 25599.5 498758 6.08385e+006 213.12 D:\temp\rb\3V050909P0000897861A520004700712M_001631680\3v050909p0000897861a520004700712m_001631680.tif  
D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif
Input file size is 8016, 25600                             
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.
> copy D:\temp\rb\3V050909P0000897861A520004700712M_001631680\3v050909p0000897861a520004700712m_001631680_rpc.txt
D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org_rpc.txt
Скопировано файлов:        1.
</pre>
</pre>


Документация на gdalwarp говорит, что для орторектификации можно указать DEM-файл с помощью параметра -to 'RPC_DEM=DEM_merged.tif'. Пробуем:
Проверяем, прописалась ли информация в результирующий файл:
<pre>
<pre>
$ gdalwarp -rpc -to 'RPC_DEM=/home/ftp/maps/GDEM/DEM_merged.tif' 3v050909p0000897861a520004700712m_001631680.tif test.tif
> gdalinfo D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif                                                                               
Warning 1: TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered
Driver: GTiff/GeoTIFF                                                                           
Processing input file 3v050909p0000897861a520004700712m_001631680.tif.
Files: D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif                           
0...10...20...30...40...50...60...70...80...90...100 - done.
      D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org_rpc.txt                       
Size is 8016, 25600                                                                             
Coordinate System is `'                                                                          
GCP Projection =                                                                                
PROJCS["WGS 84 / UTM zone 35N",                                                                 
    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"]],                                                               
    PROJECTION["Transverse_Mercator"],                                                           
    PARAMETER["latitude_of_origin",0],                                                           
    PARAMETER["central_meridian",27],                                                           
    PARAMETER["scale_factor",0.9996],                                                           
    PARAMETER["false_easting",500000],                                                           
    PARAMETER["false_northing",0],                                                               
    UNIT["metre",1,                                                                             
        AUTHORITY["EPSG","9001"]],                                                               
    AUTHORITY["EPSG","32635"]]                                                                   
GCP[ 0]: Id=1, Info=                                                                         
          (0.5,0.5) -> (498793,6110750,173.385)                                                 
GCP[  1]: Id=2, Info=                                                                           
          (8015.5,0.5) -> (507345,6110270,187.386)                                               
GCP[  2]: Id=3, Info=                                                                           
          (8015.5,25599.5) -> (507347,6083510,204.881)                                           
GCP[  3]: Id=4, Info=                                                                           
          (0.5,25599.5) -> (498758,6083850,213.12)                                              
Metadata:                                                                                       
  AREA_OR_POINT=Area                                                                             
  TIFFTAG_MAXSAMPLEVALUE=2047                                                                   
  TIFFTAG_MINSAMPLEVALUE=0                                                                       
Image Structure Metadata:                                                                       
  INTERLEAVE=BAND                                                                               
RPC Metadata:                                                                                   
  HEIGHT_OFF= +0179.000 meters                                                                   
  HEIGHT_SCALE= +0300.000 meters                                                                 
  LAT_OFF= +55.02030000 degrees                                                                 
  LAT_SCALE= +00.12380000 degrees                                                               
  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                                                                                           
  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_OFF= +012800.00 pixels                                                                   
  LINE_SCALE= +012800.00 pixels                                                       
  LONG_OFF= +027.04780000 degrees                                                               
  LONG_SCALE= +000.06850000 degrees                                                             
  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                                                                                           
  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_OFF= +004008.00 pixels                                                                   
  SAMP_SCALE= +004008.00 pixels                                                                 
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
</pre>
</pre>


Загружаем в QGIS — совсем другой коленкор! Хотя всё равно получили погрешность порядка 10 метров. Причины слёту не выяснены, да и иметь растр, "перекрученный" в географическую СК, как-то не очень хочется. Как трансформировать с использованием RPC напрямую в, например, UTM zone 35, понять пока не удалось.
Теперь можно повторить команду с новым файлом изображения.
  > gdalwarp -dstnodata 0 -srcnodata 0 -overwrite -t_srs epsg:32635 -wo "INIT_DEST=NO_DATA" -rpc -to "RPC_DEM=D:\temp\rb\DEM_merged.tif"
    D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif D:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec.tif
 
Модель RPC использует высотные отметки в системе координат WGS84, а большинство данных о рельефе (DEM) используют относительные отметки геоида. Для устранения этого расхождения необходимо определить превышения геоида над эллипсоидом WGS84 и использовать это значение в команде ортокоррекции. Для определения превышения, возьмем координаты средней точки снимка (27d 2'55.34"E, 55d 1'14.26"N) и загрузим на этот онлайн [http://www.unavco.org/community_science/science-support/geoid/geoid.html ресурс] или на этот онлайн [http://geographiclib.sourceforge.net/cgi-bin/GeoidEval ресурс] (такой расчет может быть выполнен как с помощью онлайн ресурсов, так и с использованием программного обеспечения, например, proj4). В результате получим '''22.0157 м''' для модели EGM96.
 
Для ортокоррекции с учетом геоида выполним следующую команду.
  > gdalwarp -dstnodata 0 -srcnodata 0 -overwrite -t_srs epsg:32635 -wo "INIT_DEST=NO_DATA" -rpc -to "RPC_DEM=D:\temp\rb\DEM_merged.tif" -to "RPC_HEIGHT=22.0157"
    D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif D:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec2.tif
 
Для трансформации в другую систему координат просто добавьте -t_srs proj_definition. Если при этом gdalwarp отказывается трансформировать что-либо, скачайте соседние тайлы GDEM - скорее всего, "поля" в трансформированном растре "выскакивают" к соседям.
 
== Выводы ==
 
Для оценки качества выполнения ортокорреции с использованием GDAL воспользуемся коммерческим программным обеспечением ENVI EX. Выполним в этом программном продукте аналогичные операции и сравним результаты. Приведу изображение одного и того же участка местности с наложенными на него GPS треками.
 
[[Файл:Gdal vs envi.png|740px]]
 
На изображении видно, что при выполнении ортокоррекции без учета геоида отмечается сдвиг трека в сторону от дороги. Результат выполнения ортокоррекции с учетом геоида полностью идентичен результатам работы ПО ENVI EX.
 
Однако, на изображении, полученном в результате ортокорреции при помощи GDAL, имеются артефакты. Причем, ортокорреция, выполненная в программе wxGIS, тоже основанной на библиотеке GDAL, не приводит к появлению артефактов. <strike>Данное необъяснимое поведение библиотеки (версия GDAL 1.9.0, released 2011/12/29) требует отдельного изучения</strike>. Для устранения артефактов в командую строку необходимо добавить ключ -et 0.0. Пример команды:
 
  > gdalwarp -dstnodata 0 -srcnodata 0 -overwrite -t_srs epsg:32635 -wo "INIT_DEST=NO_DATA" -et 0.0 -rpc -to "RPC_DEM=D:\temp\rb\DEM_merged.tif" -to "RPC_HEIGHT=22.0157"
    D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif D:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec2.tif
 
Данный ключ отключает излишнюю оптимизацию при интерполяции растра.
 
[[Файл:Gdal ortho paradox.png|740px]]


Зато, если скормить исходный растр и DEM программе (увы, платной и под Windows) ENVI EX, она после длительного пережёвывания на выходе даст практически идеально привязанный файл в проекции UTM нужной зоны. Пример можно увидеть здесь: http://latlon.org/~jek/osm/ov3-sample.jpg . Если вы знаете, как достичь таких же результатов с помощью GDAL — поделитесь.
== Ссылки ==
* Оригинальный снимок 3v050909p0000897861a520004700712m_001631680, скачанный через EarthSat ([http://gis-lab.info/data/orbview-3/samples/rb/3v050909p0000897861a520004700712m_001631680.zip скачать])
* Фрагмент ASTER GDEM использованный для ортокоррекции ([http://gis-lab.info/data/orbview-3/samples/rb/DEM_merged.tif скачать])
* Треки (GPX) для проверки результата ([http://gis-lab.info/data/orbview-3/samples/rb/tracks.7z скачать] и [http://gis-lab.info/data/orbview-3/samples/rb/ov3-check.7z скачать])
* Результат ортокоррекции в ENVI EX ([http://gis-lab.info/data/orbview-3/samples/rb/result_envi.tif скачать])
* Результат ортокоррекции в GDAL ([http://gis-lab.info/data/orbview-3/samples/rb/result_gdal.7z скачать])
* Результат ортокоррекции в GDAL с коррекцией геоида ([http://gis-lab.info/data/orbview-3/samples/rb/result_gdal_geoid_corrected.7z скачать])

Текущая версия от 17:48, 17 июня 2012

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


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

Введение

В данной статье рассматривается способ ортокоррекции космических снимков с КА OrbView-3 при помощи библиотеки GDAL.

Космические снимки с КА OrbView-3 свободно распространяются и могут быть получены через каталог данных OrbView-3. Однако, в исходном виде, снимки не имеют географической привязки - в формат распространения (TIFF) не внедрены специальные теги для описания привязки, т.е. формат не представляет из себя классического GeoTIFF. Однако, в комплекте поставки имеются все необходимые метаданные для осуществления привязки и ортокоррекции снимков. Как это осуществить будет рассказано далее.

Программное обеспечение и данные

Для выполнения ортокоррекции космического снимка необходимо:

  1. Свободно распространяемое программное обеспечение gdalwarp, gdalbuildvrt, gdal_translate из состава библиотеки GDAL. Для получения программного обеспечения можно перейти по ссылке, где перечислены варианты загрузки.
  2. Космический снимок с КА OrbView-3. Снимки на необходимую территорию ищем в каталоге.
  3. Рельеф на территорию съемки в формате в виде ЦМР (Цифровой модели рельефа)

В качестве рельефа могут быть использованы: SRTM, ASTER GDEM v.1 и v.2 и др. Рельеф загружаем отсюда:

Для демонстрации операций по ортокоррекции можно воспользоваться подготовленными наборами данных (снимки и рельеф) по районам Курской области и Республике Беларусь.

Набор данных по Курской области включает в себя:

3V050401P0000692211A520018202082M_001639429.ZIP - космический снимок с КА OrbView-3 (архив)
3v050401p0000692211a520018202082m_001639429_ortho.img - результат ортокоррекции снимка в формате ERDAS Imagine, выполненный в GDAL
3v050401p0000692211a520018202082m_001639429_ortho.img.aux.xml - " -
3v050401p0000692211a520018202082m_001639429_ortho.rrd - " -
dem.tfw - файл рельефа на территорию съемки
dem.tif - " -
dem.tif.aux.xml - " -
dem.tif.ovr - " -
dem.tif.xml - " -
kursk_ortho.7z  - результат ортокоррекции снимка в формате TIFF, выполненный в GDAL (архив)
kursk_ortho_envi_orbview.7z  - результат ортокоррекции снимка в формате TIFF, выполненный в ENVI (архив)

Набор данных по Республике Беларусь включает в себя:

3v050909p0000897861a520004700712m_001631680.zip - космический снимок с КА OrbView-3 (архив)
aster_gdem.tif - файл рельефа на территорию съемки
ov3-check.7z - наземные точки снятые с помощью GPS
result_envi.tif  - результат ортокоррекции снимка, выполненный в ENVI EX
result_gdal.7z  - результат ортокоррекции снимка, выполненный в GDAL (архив)
result_gdal_geoid_corrected.7z - результат ортокоррекции снимка с использованием коррекции геоида, выполненный в GDAL (архив)
tracks.7z - наземные точки снятые с помощью GPS

Подготовка к ортокорреции

В данной статье все примеры будут демонстрироваться на наборе данных по Республике Беларусь.

Для начала рассмотрим содержимое ZIP-архива с материалами космической съемки с КА OrbView-3 (3V050909P0000897861A520004700712M_001631680.ZIP). Для этого выведем листинг директории командой:

UNIX

$ ls -1 3V050909P0000897861A520004700712M_001631680*

Windows

> dir * /B
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, содержащий некое описание данных и их параметров.

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

$ 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 можно прочитать в статье "Ортокоррекция космических снимков с использованием RPC", а также в соответствующей теме форума. Если команда gdalinfo у вас не вывела метаданные RPC, проверьте версию - нужен GDAL не менее 1.8.1.

Прежде чем формировать команду для ортокоррекции необходимо определиться с системой координат результирующего изображения. Систему координат нужно подобрать таким образом, чтобы пиксел изображения меньше всего искажался. Для выбора системы координат, исходное изображение было привязано по углам с использованием координат из файла 3v050909p0000897861a520004700712m_001631680.pvl. После этого файл теперь уже в формате GeoTIFF был загружен в QGIS и на нем были испробованы различные системы координат с преобразованием "на лету". При этом, обращалось внимание, чтобы пиксел изображения оставался "квадратным" и вращался на наименьший угол. В результате экспериментов было принято решение об использовании системы координат WGS84/UTM соответствующей зоны.

Для нашего случая система координат будет WGS 84 / UTM zone 35N или epsg:32635.

Получение необходимых данных о рельефе рассмотрим на примере работы с данными SRTM. Для получения файла рельефа необходимо:

  1. На странице SRTM Data Selection Options выписываем номера рядов и столбцов с территорией съемки.
  2. Скачиваем с http://gis-lab.info/data/srtm-tif выбранные zip-файлы и складываем в один каталог.
  3. Распаковываем (команда для операционной системы UNIX)
  4. for i in srtm*zip; do yes|unzip $i; done
  5. Объединяем все файлы SRTM в формате GeoTIFF в единый виртуальный растр
  6. Команда для операционной системы UNIX: gdalbuildvrt srtm.vrt srtm*tif Команда для операционной системы Windows: gdalbuildvrt.exe srtm.vrt srtm*tif

Для получения файла рельефа ASTER GDEM переходим на сайт Aster GDEM, регистрируемся и заходим в раздел Search. Выбираем "Select tiles by shapefile", загружаем файл покрытия scene.shp, скачиваем и распаковываем. Получаем несколько (в данном случае — 4) файлов с данными DEM с именами вида ASTGTM2_N55E026_dem.tif (один тайл соответствует "квадрату" градусной сетки). Склеиваем все файлы в один GeoTIFF с помощью gdal_merge.py:

$ 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.

Ортокорреция

Для выполнения ортокоррекции космического снимка с КА OrbView-3 необходимо выполнить следующую команду.

 > gdalwarp -dstnodata 0 -srcnodata 0 -overwrite -t_srs epsg:32635 -wo "INIT_DEST=NO_DATA" -rpc -to "RPC_DEM=D:\temp\rb\DEM_merged.tif"
 D:\temp\rb\3V050909P0000897861A520004700712M_001631680\3v050909p0000897861a520004700712m_001631680.tif D:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec.tif 

Параметры команды gdalwarp можно посмотреть по следующей ссылке, подробнее о параметрах -wo можно почитать здесь.

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

Вариант №1. Попробовать извлечь из метаданных. Минус варианта - выходная система координат будет WGS84, что приводит к потере разрешения (размера пиксела на местности).

$ 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.

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

$ 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

Вариант №2. Прописать систему координат и GCP углов изображения. Для этого в файле scene.pvl ищем координаты углов в СК WGS84, пересчитываем их в WGS84/UTM 35N и добавляем их к файлу.

> gdal_translate -a_srs epsg:32635 -gcp 0.5 0.5 498793 6.11075e+006 173.385 -gcp 8015.5 0.5 507345 6.11027e+006 187.386 -gcp 8015.5 25599.5 507347 6.08351e+006 204.881 
-gcp 0.5 25599.5 498758 6.08385e+006 213.12 D:\temp\rb\3V050909P0000897861A520004700712M_001631680\3v050909p0000897861a520004700712m_001631680.tif 
D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif
Input file size is 8016, 25600                              
0...10...20...30...40...50...60...70...80...90...100 - done.

> copy D:\temp\rb\3V050909P0000897861A520004700712M_001631680\3v050909p0000897861a520004700712m_001631680_rpc.txt 
D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org_rpc.txt
Скопировано файлов:         1. 

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

> gdalinfo D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif                                                                                 
Driver: GTiff/GeoTIFF                                                                             
Files: D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif                             
       D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org_rpc.txt                         
Size is 8016, 25600                                                                               
Coordinate System is `'                                                                           
GCP Projection =                                                                                  
PROJCS["WGS 84 / UTM zone 35N",                                                                   
    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"]],                                                                
    PROJECTION["Transverse_Mercator"],                                                            
    PARAMETER["latitude_of_origin",0],                                                            
    PARAMETER["central_meridian",27],                                                             
    PARAMETER["scale_factor",0.9996],                                                             
    PARAMETER["false_easting",500000],                                                            
    PARAMETER["false_northing",0],                                                                
    UNIT["metre",1,                                                                               
        AUTHORITY["EPSG","9001"]],                                                                
    AUTHORITY["EPSG","32635"]]                                                                    
GCP[  0]: Id=1, Info=                                                                          
          (0.5,0.5) -> (498793,6110750,173.385)                                                   
GCP[  1]: Id=2, Info=                                                                             
          (8015.5,0.5) -> (507345,6110270,187.386)                                                
GCP[  2]: Id=3, Info=                                                                             
          (8015.5,25599.5) -> (507347,6083510,204.881)                                            
GCP[  3]: Id=4, Info=                                                                             
          (0.5,25599.5) -> (498758,6083850,213.12)                                                
Metadata:                                                                                         
  AREA_OR_POINT=Area                                                                              
  TIFFTAG_MAXSAMPLEVALUE=2047                                                                     
  TIFFTAG_MINSAMPLEVALUE=0                                                                        
Image Structure Metadata:                                                                         
  INTERLEAVE=BAND                                                                                 
RPC Metadata:                                                                                     
  HEIGHT_OFF= +0179.000 meters                                                                    
  HEIGHT_SCALE= +0300.000 meters                                                                  
  LAT_OFF= +55.02030000 degrees                                                                   
  LAT_SCALE= +00.12380000 degrees                                                                 
  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                                                                                            
  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_OFF= +012800.00 pixels                                                                     
  LINE_SCALE= +012800.00 pixels                                                         
  LONG_OFF= +027.04780000 degrees                                                                 
  LONG_SCALE= +000.06850000 degrees                                                               
  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                                                                                            
  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_OFF= +004008.00 pixels                                                                     
  SAMP_SCALE= +004008.00 pixels                                                                   
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 

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

 > gdalwarp -dstnodata 0 -srcnodata 0 -overwrite -t_srs epsg:32635 -wo "INIT_DEST=NO_DATA" -rpc -to "RPC_DEM=D:\temp\rb\DEM_merged.tif"
   D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif D:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec.tif 

Модель RPC использует высотные отметки в системе координат WGS84, а большинство данных о рельефе (DEM) используют относительные отметки геоида. Для устранения этого расхождения необходимо определить превышения геоида над эллипсоидом WGS84 и использовать это значение в команде ортокоррекции. Для определения превышения, возьмем координаты средней точки снимка (27d 2'55.34"E, 55d 1'14.26"N) и загрузим на этот онлайн ресурс или на этот онлайн ресурс (такой расчет может быть выполнен как с помощью онлайн ресурсов, так и с использованием программного обеспечения, например, proj4). В результате получим 22.0157 м для модели EGM96.

Для ортокоррекции с учетом геоида выполним следующую команду.

 > gdalwarp -dstnodata 0 -srcnodata 0 -overwrite -t_srs epsg:32635 -wo "INIT_DEST=NO_DATA" -rpc -to "RPC_DEM=D:\temp\rb\DEM_merged.tif" -to "RPC_HEIGHT=22.0157"
   D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif D:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec2.tif

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

Выводы

Для оценки качества выполнения ортокорреции с использованием GDAL воспользуемся коммерческим программным обеспечением ENVI EX. Выполним в этом программном продукте аналогичные операции и сравним результаты. Приведу изображение одного и того же участка местности с наложенными на него GPS треками.

Gdal vs envi.png

На изображении видно, что при выполнении ортокоррекции без учета геоида отмечается сдвиг трека в сторону от дороги. Результат выполнения ортокоррекции с учетом геоида полностью идентичен результатам работы ПО ENVI EX.

Однако, на изображении, полученном в результате ортокорреции при помощи GDAL, имеются артефакты. Причем, ортокорреция, выполненная в программе wxGIS, тоже основанной на библиотеке GDAL, не приводит к появлению артефактов. Данное необъяснимое поведение библиотеки (версия GDAL 1.9.0, released 2011/12/29) требует отдельного изучения. Для устранения артефактов в командую строку необходимо добавить ключ -et 0.0. Пример команды:

 > gdalwarp -dstnodata 0 -srcnodata 0 -overwrite -t_srs epsg:32635 -wo "INIT_DEST=NO_DATA" -et 0.0 -rpc -to "RPC_DEM=D:\temp\rb\DEM_merged.tif" -to "RPC_HEIGHT=22.0157"
   D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif D:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec2.tif

Данный ключ отключает излишнюю оптимизацию при интерполяции растра.

Gdal ortho paradox.png

Ссылки

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