Ортокоррекция данных OrbView-3 с помощью GDAL: различия между версиями
мНет описания правки |
Bishop (обсуждение | вклад) (Исправление про артефакты) |
||
(не показана 21 промежуточная версия 4 участников) | |||
Строка 1: | Строка 1: | ||
{{Статья| | {{Статья|Опубликована|orbview3-ortho-gdal}} | ||
{{Аннотация|Как произвести ортокоррекцию данных OrbView-3 с помощью GDAL}} | |||
== | == Введение == | ||
В данной статье рассматривается способ ортокоррекции космических снимков с КА OrbView-3 при помощи библиотеки [http://gdal.org GDAL]. | |||
Космические снимки с КА OrbView-3 свободно распространяются и могут быть получены через [[Каталог данных OrbView-3|каталог данных OrbView-3]]. | |||
Однако, в исходном виде, снимки не имеют географической привязки - в формат распространения (TIFF) не внедрены специальные теги для описания привязки, т.е. формат не представляет из себя классического GeoTIFF. Однако, в комплекте поставки имеются все необходимые метаданные для осуществления привязки и ортокоррекции снимков. Как это осуществить будет рассказано далее. | |||
== Программное обеспечение и данные == | |||
Для выполнения ортокоррекции космического снимка необходимо: | |||
# Свободно распространяемое программное обеспечение [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|каталоге]]. | |||
# Рельеф на территорию съемки в формате в виде ЦМР (Цифровой модели рельефа) | |||
В качестве рельефа могут быть использованы: [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] | |||
Для демонстрации операций по ортокоррекции можно воспользоваться подготовленными наборами данных (снимки и рельеф) по районам [http://gis-lab.info/data/orbview-3/samples/kursk/ Курской области] и [http://gis-lab.info/data/orbview-3/samples/rb/ Республике Беларусь]. | |||
Набор данных по Курской области включает в себя: | |||
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 | |||
<pre> | <pre> | ||
3v050909p0000897861a520004700712m_001631680_aoi.dbf | 3v050909p0000897861a520004700712m_001631680_aoi.dbf | ||
3v050909p0000897861a520004700712m_001631680_aoi.prj | 3v050909p0000897861a520004700712m_001631680_aoi.prj | ||
Строка 61: | Строка 77: | ||
</pre> | </pre> | ||
Имеем: шейп-файлы с покрытием, jpg с "превьюшкой" (привязанный!), собственно TIFF с данными, файл с параметрами RPC преобразований (scene_rpc.txt), файл scene.pvl, содержащий некое описание данных и их параметров. | Имеем: шейп-файлы с покрытием, файл jpg с "превьюшкой" (привязанный!), собственно файл TIFF с данными, файл с параметрами RPC преобразований (scene_rpc.txt), файл scene.pvl, содержащий некое описание данных и их параметров. | ||
Если загрузить | Если загрузить TIFF-файл с данными в QGIS, то никакой привязкой там пахнуть и не будет. | ||
Попробуем разобраться: | Попробуем разобраться: | ||
<pre> | <pre> | ||
Строка 102: | Строка 118: | ||
</pre> | </pre> | ||
Как и ожидалось, в файле не содержится данных о привязке. Зато GDAL прочитал данные RPC (rational polynomial coefficients), нужные для корректной привязки и трансформации. Подробнее про | Как и ожидалось, в файле не содержится данных о привязке. Зато 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> | |||
Команда для операционной системы 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 | ||
Строка 113: | Строка 164: | ||
0...10...20...30...40...50...60...70...80...90...100 - done.</pre> | 0...10...20...30...40...50...60...70...80...90...100 - done.</pre> | ||
Проверяем, прописалась ли информация в результирующий файл: | |||
<pre>$ gdalinfo test.tif | <pre>$ gdalinfo test.tif | ||
Строка 143: | Строка 194: | ||
Band 1 Block=12925x1 Type=UInt16, ColorInterp=Gray</pre> | Band 1 Block=12925x1 Type=UInt16, ColorInterp=Gray</pre> | ||
Вариант №2. Прописать систему координат и GCP углов изображения. Для этого в файле scene.pvl ищем координаты углов в СК WGS84, пересчитываем их в WGS84/UTM 35N и добавляем их к файлу. | |||
<pre> | <pre> | ||
> 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> | ||
Проверяем, прописалась ли информация в результирующий файл: | |||
<pre> | <pre> | ||
> gdalinfo D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif | |||
Driver: GTiff/GeoTIFF | |||
Files: D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org.tif | |||
0... | 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> | ||
Теперь можно повторить команду с новым файлом изображения. | |||
> 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]] | |||
[[Файл: | |||
== Ссылки == | |||
* Оригинальный снимок 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. Однако, в комплекте поставки имеются все необходимые метаданные для осуществления привязки и ортокоррекции снимков. Как это осуществить будет рассказано далее.
Программное обеспечение и данные
Для выполнения ортокоррекции космического снимка необходимо:
- Свободно распространяемое программное обеспечение gdalwarp, gdalbuildvrt, gdal_translate из состава библиотеки GDAL. Для получения программного обеспечения можно перейти по ссылке, где перечислены варианты загрузки.
- Космический снимок с КА OrbView-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. Для получения файла рельефа необходимо:
- На странице SRTM Data Selection Options выписываем номера рядов и столбцов с территорией съемки.
- Скачиваем с http://gis-lab.info/data/srtm-tif выбранные zip-файлы и складываем в один каталог.
- Распаковываем (команда для операционной системы UNIX) for i in srtm*zip; do yes|unzip $i; done
- Объединяем все файлы SRTM в формате GeoTIFF в единый виртуальный растр Команда для операционной системы 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 треками.
На изображении видно, что при выполнении ортокоррекции без учета геоида отмечается сдвиг трека в сторону от дороги. Результат выполнения ортокоррекции с учетом геоида полностью идентичен результатам работы ПО 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
Данный ключ отключает излишнюю оптимизацию при интерполяции растра.
Ссылки
- Оригинальный снимок 3v050909p0000897861a520004700712m_001631680, скачанный через EarthSat (скачать)
- Фрагмент ASTER GDEM использованный для ортокоррекции (скачать)
- Треки (GPX) для проверки результата (скачать и скачать)
- Результат ортокоррекции в ENVI EX (скачать)
- Результат ортокоррекции в GDAL (скачать)
- Результат ортокоррекции в GDAL с коррекцией геоида (скачать)