GIS-Lab - Вклад [ru]

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

2012-01-15T22:27:05Z

Jekhor: /* Развёрнутое описание */

Снимки [[Каталог_данных_Orbview3|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:
* Узнать, как в том же шаге сделать internal nodata mask, для более аккуратных краёв изображения;

=== Развёрнутое описание ===

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

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

Для тех, кто не хочет читать длинные рассуждения с примерами — краткие выводы:
#gdalwarp без DEM — очень грубая привязка.
#gdalwarp с DEM — получше, но всё равно неточно.
#ENVI EX с DEM — практически идеально, только дорого.

Сравнение результатов ENVI EX + Aster, GDAL + aster и GDAL + SRTM: http://latlon.org/~jek/osm/ov3-ortho-comparison.jpg

Для начала, заглянем в содержимое ZIP-архива с одной сценой, скачанного с сайта USGS:
<pre>
$ 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
</pre>

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

Если загрузить, например, в qgis, TIFF-файл с данными, то никакой привязкой там пахнуть и не будет.
Попробуем разобраться:
<pre>
$ 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
</pre>

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

Для корректной орторектификации требуются параметры RPC и информация о высоте над геоидом. Что-то там про высоту в выводе gdalinfo есть (интересующиеся могут разобраться с этим самостоятельно). Попробуем в лоб:
<pre>
$ 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
</pre>

Ага, появились вполне правдоподобные данные привязки (в географической СК, правда). Загружаем полученный файл в 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:
<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>

Документация на gdalwarp говорит, что для орторектификации можно указать DEM-файл с помощью параметра -to 'RPC_DEM=DEM_merged.tif'. Пробуем:
<pre>
$ 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.
</pre>

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

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

Зато, если скормить исходный растр и DEM программе (увы, платной и под Windows) ENVI EX, она после длительного пережёвывания на выходе даст практически идеально привязанный файл в проекции UTM нужной зоны. Пример можно увидеть здесь: http://latlon.org/~jek/osm/ov3-sample.jpg . Если вы знаете, как достичь таких же результатов с помощью GDAL — поделитесь.

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

2012-01-15T22:11:54Z

Jekhor:

Снимки [[Каталог_данных_Orbview3|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:
* Узнать, как в том же шаге сделать internal nodata mask, для более аккуратных краёв изображения;

=== Развёрнутое описание ===

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

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

Для тех, кто не хочет читать длинные рассуждения с примерами — краткие выводы:
#gdalwarp без DEM — очень грубая привязка.
#gdalwarp с DEM — получше, но всё равно неточно.
#ENVI EX с DEM — практически идеально, только дорого.

Для начала, заглянем в содержимое ZIP-архива с одной сценой, скачанного с сайта USGS:
<pre>
$ 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
</pre>

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

Если загрузить, например, в qgis, TIFF-файл с данными, то никакой привязкой там пахнуть и не будет.
Попробуем разобраться:
<pre>
$ 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
</pre>

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

Для корректной орторектификации требуются параметры RPC и информация о высоте над геоидом. Что-то там про высоту в выводе gdalinfo есть (интересующиеся могут разобраться с этим самостоятельно). Попробуем в лоб:
<pre>
$ 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
</pre>

Ага, появились вполне правдоподобные данные привязки (в географической СК, правда). Загружаем полученный файл в 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:
<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>

Документация на gdalwarp говорит, что для орторектификации можно указать DEM-файл с помощью параметра -to 'RPC_DEM=DEM_merged.tif'. Пробуем:
<pre>
$ 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.
</pre>

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

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

Зато, если скормить исходный растр и DEM программе (увы, платной и под Windows) ENVI EX, она после длительного пережёвывания на выходе даст практически идеально привязанный файл в проекции UTM нужной зоны. Пример можно увидеть здесь: http://latlon.org/~jek/osm/ov3-sample.jpg . Если вы знаете, как достичь таких же результатов с помощью GDAL — поделитесь.

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

2012-01-13T22:23:40Z

Jekhor:

Снимки [[Каталог_данных_Orbview3|Orbview-3]] доступны в исходном виде, без орторектификации. Чтобы привести их в состояние, пригодное для работы, а не только для любования, нужно орторектификацию выполнить.

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

Для тех, кто не хочет читать длинные рассуждения с примерами — краткие выводы:
#gdalwarp без DEM — очень грубая привязка.
#gdalwarp с DEM — получше, но всё равно неточно.
#ENVI EX с DEM — практически идеально, только дорого.

Для начала, заглянем в содержимое ZIP-архива с одной сценой, скачанного с сайта USGS:
<pre>
$ 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
</pre>

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

Если загрузить, например, в qgis, TIFF-файл с данными, то никакой привязкой там пахнуть и не будет.
Попробуем разобраться:
<pre>
$ 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
</pre>

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

Для корректной орторектификации требуются параметры RPC и информация о высоте над геоидом. Что-то там про высоту в выводе gdalinfo есть (интересующиеся могут разобраться с этим самостоятельно). Попробуем в лоб:
<pre>
$ 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
</pre>

Ага, появились вполне правдоподобные данные привязки (в географической СК, правда). Загружаем полученный файл в 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:
<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>

Документация на gdalwarp говорит, что для орторектификации можно указать DEM-файл с помощью параметра -to 'RPC_DEM=DEM_merged.tif'. Пробуем:
<pre>
$ 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.
</pre>

Загружаем в QGIS — совсем другой коленкор! Хотя всё равно получили погрешность порядка 10 метров. Причины слёту не выяснены, да и иметь растр, "перекрученный" в географическую СК, как-то не очень хочется. Как трансформировать с использованием RPC напрямую в, например, UTM zone 35, понять пока не удалось.

Зато, если скормить исходный растр и DEM программе (увы, платной и под Windows) ENVI EX, она после длительного пережёвывания на выходе даст практически идеально привязанный файл в проекции UTM нужной зоны. Пример можно увидеть здесь: http://latlon.org/~jek/osm/ov3-sample.jpg . Если вы знаете, как достичь таких же результатов с помощью GDAL — поделитесь.

Каталог данных OrbView-3

2012-01-13T22:17:42Z

Jekhor:

<p class="ann">Описание индексного каталога данных Orbview-3 свободно доступных через EarthExplorer USGS</p>

9 января USGS [http://www.usgs.gov/newsroom/article.asp?ID=3070 объявило] о том, что стали общедоступны данные OrbView-3 за 2003-2007 гг. Доступ к самим данным дистанционного зондирования осуществляется через инструмент [http://earthexplorer.usgs.gov/ EarthExplorer].

Данные OrbView-3 полученные USGS от компании GeoEye ограничены периодом Сентрябрь 2003 - Март 2007. Платформа OrbView-3 была выведена GeoEye из рабочего состояния 23 апреля 2007 из-за технических проблем на борту. Орбита OrbView-3 постепенно снижалась пока аппарат не сгорел над Тихим Океаном 13 марта 2011. USGS предложило спасти данные OrbView-3 после того, как GeoEye объявила, что она перестанет поддерживать доступ к этим данным в 2011 году. Исследование проведенное USGS показало ценность данных и то, что данные дополняли другие наборы поддерживаемые EROS. В настоящий момент, данные проходят процедуру улучшения, это займет время, но через некоторое время появятся улучшенные продукты на базе данных доступных в текущий момент.

Данные не орторектифицированы. На текущий момент, орторектификация выполняется только для отдельных сцен по запросу (on-demand). Одно из препятствий - отсутствие системных данных по наземной привязке. Единственные данные, которые доступны на данный момент - точки привязки Landsat с точностью +/- 30 метров.

Немного про орторектификацию по имеющимся данным можно найти [[Орторектификация_Orbview3|здесь]].

=== Лицензия ===
По видимому из-за организационной ошибки на данный момент (13.01.2011) к данным прилагается лицензия очень сильно ограничивающая свободу использования данных. Ознакомиться с [http://gis-lab.info/data/orbview-meta/unrestricted_license.txt лицензией] прилагаемой на данный момент. Текущая лицензия является более запретительной чем даже коммерческая лицензия GeoEye ([http://gis-lab.info/data/orbview-meta/geoeye_license_commercial_2007.txt пример] коммерческой лицензии).

Согласно заявлениям по крайней мере двух сотрудников USGS лицензия данных - Public Domain, что означает что данные можно свободно копировать, распространять, использовать для создания производных продуктов и продавать ([http://lists.openstreetmap.org/pipermail/talk/2012-January/061546.html ссылка]).

<pre>-------- Original Message --------
Subject: Re: Fw: OrbView-3 free images?
Date: Wed, 11 Jan 2012 15:18:50 -0600
From: GS-N-EDC Customer Services <custserv at usgs.gov>
CC: Karen A Peterson <kapeterson at usgs.gov>

...
USGS EROS received approximately 180,000 OrbView-3 image "segments'. The data that was delivered by GeoEye was a L1B "Basic Enhanced" radiometrically corrected format.

It is now available to the general public for download at no charge from the EarthExplorer site. http://earthexplorer.usgs.gov

Yes, it is public domain. There is no legal requirement for users to acknowledge or credit USGS as the source for public domain information, but they may wish to do so as a courtesy. If you wish to acknowledge or credit USGS as an information source of data or products, use a line of text as show below. _http://eros.usgs.gov/#/About_Us/Customer_Service/Data_Citation_

"*Data available from the U.S. Geological Survey*"

Please see this page for further information and details: _http://www.usgs.gov/visual-id/credit_usgs.html_
...
</pre>

=== Каталог данных ===
Здесь можно получить слой, который можно использовать в ГИС, чтобы получить быстрый ответ на вопрос, покрыта ли моя территория свободнодоступными снимками OrbView. Слой показывает данные на весь мир.
Каталог данных можно скачать как:
*слой в формате Shape http://gis-lab.info/data/orbview-meta/orbview3-catalog-shp.7z (~17 Mb)
*слой в формате OSM ''еще не готов''
Загрузить скрипт, использовавшийся для создания каталога в формате Shape можно [http://gis-lab.info/programs/python/orbview3_catalog_maker.7z здесь]

[[Файл:Orbview3-coverage.gif|700px|thumb|center]]

=== Метаданные ===
Все метаданные по блокам 20х20 градусов можно скачать здесь: http://gis-lab.info/data/orbview-meta

Все метаданные одним архивом: http://gis-lab.info/data/orbview-meta/orbview3-meta.7z

=== Разграфка блоков ===
Слой каталога был получен посредством импорта метаданных полученных через сервис EarthExplorer. Получение данных осуществлялось блоками 20 на 20 градусов, чтобы избежать ограничения на 10 тыс. сцен на запрос.
[[Файл:Grid20deg.gif|700px|thumb|center]]

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

2012-01-13T21:47:38Z

Jekhor:

Снимки [[Каталог_данных_Orbview3|Orbview-3]] доступны в исходном виде, без орторектификации. Чтобы привести их в состояние, пригодное для работы, а не только для любования, нужно орторектификацию выполнить.

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

Для тех, кто не хочет читать длинные рассуждения с примерами — краткое выводы:
gdalwarp без DEM — очень грубая привязка.
gdalwarp с DEM — получше, но всё равно неточно.
ENVI EX с DEM — практически идеально, только дорого.

Для начала, заглянем в содержимое ZIP-архива с одной сценой, скачанного с сайта USGS:
<pre>
$ 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
</pre>

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

Если загрузить, например, в qgis, TIFF-файл с данными, то никакой привязкой там пахнуть и не будет.
Попробуем разобраться:
<pre>
$ 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
</pre>

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

Для корректной орторектификации требуются параметры RPC и информация о высоте над геоидом. Что-то там про высоту в выводе gdalinfo есть (интересующиеся могут разобраться с этим самостоятельно). Попробуем в лоб:
<pre>
$ 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
</pre>

Ага, появились вполне правдоподобные данные привязки (в географической СК, правда). Загружаем полученный файл в 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:
<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>

Документация на gdalwarp говорит, что для орторектификации можно указать DEM-файл с помощью параметра -to 'RPC_DEM=DEM_merged.tif'. Пробуем:
<pre>
$ 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.
</pre>

Загружаем в QGIS — совсем другой коленкор! Хотя всё равно получили погрешность порядка 10 метров. Причины слёту не выяснены, да и иметь растр, "перекрученный" в географическую СК, как-то не очень хочется. Как трансформировать с использованием RPC напрямую в, например, UTM zone 35, понять пока не удалось.

Зато, если скормить исходный растр и DEM программе (увы, платной и под Windows) ENVI EX, она после длительного пережёвывания на выходе даст практически идеально привязанный файл в проекции UTM нужной зоны. Пример можно увидеть здесь: http://latlon.org/~jek/osm/ov3-sample.jpg . Если вы знаете, как достичь таких же результатов с помощью GDAL — поделитесь.

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

2012-01-13T21:46:58Z

Jekhor: Новая страница: «Снимки Orbview-3 доступны в исходном виде, без орторектификации. Чтобы п...»

Снимки [[Каталог_данных_Orbview|Orbview-3]] доступны в исходном виде, без орторектификации. Чтобы привести их в состояние, пригодное для работы, а не только для любования, нужно орторектификацию выполнить.

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

Для тех, кто не хочет читать длинные рассуждения с примерами — краткое выводы:
gdalwarp без DEM — очень грубая привязка.
gdalwarp с DEM — получше, но всё равно неточно.
ENVI EX с DEM — практически идеально, только дорого.

Для начала, заглянем в содержимое ZIP-архива с одной сценой, скачанного с сайта USGS:
<pre>
$ 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
</pre>

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

Если загрузить, например, в qgis, TIFF-файл с данными, то никакой привязкой там пахнуть и не будет.
Попробуем разобраться:
<pre>
$ 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
</pre>

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

Для корректной орторектификации требуются параметры RPC и информация о высоте над геоидом. Что-то там про высоту в выводе gdalinfo есть (интересующиеся могут разобраться с этим самостоятельно). Попробуем в лоб:
<pre>
$ 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
</pre>

Ага, появились вполне правдоподобные данные привязки (в географической СК, правда). Загружаем полученный файл в 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:
<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>

Документация на gdalwarp говорит, что для орторектификации можно указать DEM-файл с помощью параметра -to 'RPC_DEM=DEM_merged.tif'. Пробуем:
<pre>
$ 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.
</pre>

Загружаем в QGIS — совсем другой коленкор! Хотя всё равно получили погрешность порядка 10 метров. Причины слёту не выяснены, да и иметь растр, "перекрученный" в географическую СК, как-то не очень хочется. Как трансформировать с использованием RPC напрямую в, например, UTM zone 35, понять пока не удалось.

Зато, если скормить исходный растр и DEM программе (увы, платной и под Windows) ENVI EX, она после длительного пережёвывания на выходе даст практически идеально привязанный файл в проекции UTM нужной зоны. Пример можно увидеть здесь: http://latlon.org/~jek/osm/ov3-sample.jpg . Если вы знаете, как достичь таких же результатов с помощью GDAL — поделитесь.