OrbView-3 Orthorectification
Orbview-3 imagery is not georeferences and orthorectified. To make it usable, one needs to do all that.
Instructions
There are three easy steps to get a small 8bit raster images from the source tiffs:
- Download required SRTM regions as GeoTIFF, for imagery orthorectification. To do that:
- on the http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp page look at the picture and write down identificators or rows and columns you need.
- download zip files with those numbers from http://gis-lab.info/data/srtm-tif
- unpack them:
for i in srtm*zip; do yes|unzip $i; done
- For convinience let's join all SRTM GeoTIFF files into a single raster image:
gdalbuildvrt srtm.vrt srtm*tif
- To make a small 8bit GeoTIFF from the source files, run the following command:
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:
- Research how to do an "internal nodata mask", to get more accurate image edges;
Extended instructions
In this section other methods will be discussed.
If you notice the author to be missing some things very obvious for a specialist (the author is amateur), please fix it.
For those who don't like to read long explanations with examples - short summary:
- gdalwarp without DEM - very rough referencing.
- gdalwarp with DEM - better, but still not accurate.
- ENVI EX with DEM - almost perfect, but expensive.
A comparison of results from ENVI EX + Aster, GDAL + aster and GDAL + SRTM: http://latlon.org/~jek/osm/ov3-ortho-comparison.jpg
For a start, let's take a look inside ZIP-archive with a single scene, downloaded from the USGS website:
$ 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
We have: a shape-file with coverage, jpg with "preview" (georeferenced!), TIFF file itself with data, a file with RPC conversion parameters (scene_rpc.txt), and a scene.pvl file, which contains something like a description of the data and its parameters.
If you open that TIFF file in QGIS, you won't find any georeferencing. Let's examine it:
$ 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
As we expected, there are no referencing information in the file. But GDAL has succesfully read RPC (rational polynomial coefficients), which are needed for correct transformation and referencing. You can read in detail about orthorectification using RPC via links from the discussion board thread: http://gis-lab.info/forum/viewtopic.php?p=27889#p27889. If gladinfo command didn't print RPC metadata, check the version: you need GDAL not older than 1.8.1.
For correct orthorectification we need RPC parameters and information about altitude above geoid. There was something about altitude in the gdalinfo output (those who are interested can dig into that). Let's go bruteforce:
$ 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.
Now we'll check once more if the information was written into the resulting file:
$ 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
Well, there are some good-looking georeferencing data (in geodetic CS though). Let's open the file in QGIS... Hm, we hit the target region, that's for sure. But the referencing is very rough: almost hundred meters off.
Гугл рассказал, что для точной привязки нужно использовать данные рельефа DEM (digital elevation model). ОК. Идём на сайт Aster GDEM (http://www.gdem.aster.ersdac.or.jp/index.jsp), регистрируемся, заходим в раздел Search. Выбираем "Select tiles by shapefile", скармливаем файл покрытия scene.shp, скачиваем найденное, распаковываем. Получаем несколько (в данном случае — 4) файлов с данными DEM с именами вида ASTGTM2_N55E026_dem.tif (один тайл соответствует "квадрату" градусной сетки). К сожалению, ни GDAL, ни ENVI EX не умеет работать с несколькими файлами DEM для одной сцены. Поэтому склеиваем их с помощью gdal_merge.py в один GeoTIFF:
$ gdal_merge.py -o DEM_merged.tif ASTGTM2_N5[45]E02[67]/*_dem.tif 0...10...20...30...40...50...60...70...80...90...100 - done.
Документация на gdalwarp говорит, что для орторектификации можно указать DEM-файл с помощью параметра -to 'RPC_DEM=DEM_merged.tif'. Пробуем:
$ gdalwarp -rpc -to 'RPC_DEM=/home/ftp/maps/GDEM/DEM_merged.tif' 3v050909p0000897861a520004700712m_001631680.tif test.tif Warning 1: TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered Processing input file 3v050909p0000897861a520004700712m_001631680.tif. 0...10...20...30...40...50...60...70...80...90...100 - done.
Загружаем в QGIS — совсем другой коленкор! Хотя всё равно получили погрешность порядка 10 метров. Причины слёту не выяснены.
Для трансформации в другую систему координат просто добавьте -t_srs proj_definition. Если при этом gdalwarp отказывается трансформировать что-либо, скачайте соседние тайлы GDEM — скорее всего, "поля" в трансформированном растре "выскакивают" к соседям.
Зато, если скормить исходный растр и DEM программе (увы, платной и под Windows) ENVI EX, она после длительного пережёвывания на выходе даст практически идеально привязанный файл в проекции UTM нужной зоны. Пример можно увидеть здесь: http://latlon.org/~jek/osm/ov3-sample.jpg . Если вы знаете, как достичь таких же результатов с помощью GDAL — поделитесь.
