本文原文地址：http://gis-lab.info/qa/orbview3-ortho-gdal.html，借助于Google Translate工具翻译整理了一下。下文中的命令行截图是我本机测试的截图。
一、简介
本文将探讨使用GDAL来对OrbView-3卫星影像进行正射校正。
卫星图像来自免费的OrbView-3航天器，可以通过OrbView-3来了解更多信息。然而，在最原始的数据中，卫星图像是用没有地理位置的Tiff格式存储的。这里就不做详细的介绍了，就是原始的数据不是GeoTIFF格式，就是普通的TIFF格式。但是，下载的原始数据中同时包含了用于正射校正所有必要的元数据。如何利用它进行正射校正，将在稍后进行说明。
二、软件和数据
准备要进行正射校正的卫星图像：
- 免费软件，全部是GDAL的工具，用到的有gdalwarp，gdalbuildvrt和gdal_translate。要获得该软件，你可以去这个去链接下载。
- Orb-View3的卫星图像。可以先从这个链接查找大概的区域。
- 用于正射校正的DEM数据（数字高程模型）。
目前可用于正射校正的DEM数据有：SRTM，ASTER GDEM，这两种DEM的下载地址分别是：
- SRTM（http://gis-lab.info/data/srtm-tif/）
- ASTER GDEM（http://www.gdem.aster.ersdac.or.jp/index.jsp）
为了演示使用GDAL进行正射校正，选择使用的数据集（图像和DEM）是白俄罗斯共和国和库尔斯克地区（点击各自的链接会到各自的示例数据目录）。
库尔斯克地区的数据集，包括：
|
|
文件名 |
说明 |
|
1 |
3V050401P0000692211A520018202082M_001639429.ZIP |
OrbView-3w卫星数据包 |
|
2 |
3v050401p0000692211a520018202082m_001639429_ortho.img |
ERDAS IMAGINE的正射校正结果的格式 |
|
3 |
3v050401p0000692211a520018202082m_001639429_ortho.img.aux.xml |
辅助文件 |
|
4 |
3v050401p0000692211a520018202082m_001639429_ortho.rrd |
金字塔文件 |
|
5 |
dem.tfw |
DEM坐标 |
|
6 |
dem.tif |
DEM |
|
7 |
dem.tif.aux.xml |
DEM辅助数据 |
|
8 |
dem.tif.ovr |
DEM金字塔 |
|
9 |
dem.tif.xml |
DEM辅助数据 |
|
10 |
kursk_ortho.7z |
使用GDAL正射校正结果 |
|
11 |
kursk_ortho_envi_orbview.7z |
使用ENVI正射校正结果 |
白俄罗斯共和国的数据集，包括：
|
|
文件名 |
描述 |
|
1 |
3v050909p0000897861a520004700712m_001631680.zip |
OrbView-3w卫星数据包 |
|
2 |
aster_gdem.tif |
Aster的DEM数据 |
|
3 |
ov3-check.7z |
地面实测的GPS点 |
|
4 |
result_envi.tif |
ENVI-EX正射结果 |
|
5 |
result_gdal.7z |
GDAL正射结果 |
|
6 |
result_gdal_geoid_corrected.7z |
经过大地水准面修正的GDAL正射结果 |
|
7 |
tracks.7z |
地面实测GPS点轨迹 |
三、准备正射校正必要的数据
在这篇文章中，所有的例子中使用的数据是上面白俄罗斯共和国一组数据。
首先，考虑OrbView-3的数据是ZIP压缩格式（3V050909P0000897861A520004700712M_001631680.ZIP）。所以先要对其进行解压：
UNIX系统：
$ ls -13V050909P0000897861A520004700712M_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
在解压之后我们可以看到，里面有Shapefile格式的落图文件以及jpg格式的快视图数据，一个存储实际数据的TIFF文件，卫星参数参数文件，RPC系数文件（scene_rpc.txt）以及其他的一些描述文件。
如果直接用QGIS打开TIFF文件，由于TIFF文件没有投影信息，所以显示的坐标是行列号。可以用gdalinfo工具来查看详细信息：
UNIX系统：
$ gdalinfo 3v050909p0000897861a520004700712m_001631680.tifWindows系统：
>C:\warmerda\bld\bin>gdalinfo.exe3v050909p0000897861a520004700712m_001631680.tif
执行完上面的命令，会输出：
Driver: GTiff/GeoTIFF
Files: E:\GDAL正射\3v050909p0000897861a520004700712m_001631680\3v050909p0000897861a520004700712m_001631680.tif
E:\GDAL正射\3v050909p0000897861a520004700712m_001631680\3v050909p0000897861a520004700712m_001631680_rpc.txt
Size is 8016, 25600
Coordinate System is `'
Metadata:
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=8016x1Type=UInt16, ColorInterp=Gray正如预期的那样，该文件不包含数据的投影和坐标。但是，GDAL读取了数据有理多项式系数（RPC）。如果想知道更多的信息，可以在使用RPC正射卫星影像正射校正的 “ RPC “，以及在相应的论坛主题找到更多的信息【这三个链接的内容比较丰富，尤其是关于RPC的介绍的比较详细，当然也是俄语的】。
如果使用gdalinfo发现输出的元数据中没有RPC信息，请确保你的GDAL版本，至少要GDAL1.8.1以上的版本。【这段谷歌翻译的受不了，一点点都看不懂啥意思，只好自己猜了，囧……】应该可以从文件3v050909p0000897861a520004700712m_001631680.pvl中得到图像的四角经纬度坐标，该坐标系是WGS 84 / UTM区35N或EPSG：32635。
接下来获取必要的数据就是地形数据，下面以SRTM数据作为例子。获取DEM数据的步骤如下：
- 在SRTM数据选择选项填写待校正影像所在区域的行和列的数目。
- 下载http://gis-lab.info/data/srtm-tif选定的压缩文件。
- 解压缩（UNIX系统命令）,Windows下自己肯定会解压吧。
for i in srtm*zip; doyes|unzip $i; done
4. 将所有文件SRTM的GeoTIFF格式转换成一个单一的虚拟文件。使用的工具是gdalbuildvrt。命令是（第一行是UNIX命令，第二行是Windows命令）：
$gdalbuildvrt srtm.vrt SRTMTIF >gdalbuildvrt.exe srtm.vrtSRTM TIF
如果要下载ASTER的DEM文件，访问ASTER GDEM网站，进行搜索查找。选择“选择分块shape文件”，下载文件覆盖scene.shp，下载并解压。一般下载的文件名称类似ASTGTM2_N55E026_dem.tif这样的，对于一个区域，可能需要多块数据。可以使用gdal_merge.py把所有的tif文件合并为一个GeoTIFF格式的文件，命令如下：
$ Gdal_merge.py-ODEM_merged.tif ASTGTM2_N5 [45] E02 [67] / * _dem.tif > python gdal_merge.py -ODEM_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"
点击这个链接，可以查看gdalwarp工具的详细参数和帮助。如果该命令执行失败，首先检查GDAL是否安装成功，然后检查GDAL版本是否支持RPC校正，如果这两个都正常，结果还是失败，那么就是设置的图像输出坐标系有关系。
解决上面的第一个方法就是减少指定的参数。除了指定使用RPC校正之后，其余均使用默认参数，如下：
$ gdalwarp -rpc3v050909p0000897861a520004700712m_001631680.tif test.tif Warning 1:TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered Creating output file that is12925P x 23537L. Warning 1:TIFFReadDirectory:Unknown field with tag 34000 (0x84d0) encountered Processing input file3v050909p0000897861a520004700712m_001631680.tif. 0...10...20...30...40...50...60...70...80...90...100- done.
然后使用gdalinfo查看输出的结果，命令以及输出信息如下：
$ 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["WGS84",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) ( 27d2'55.34"E, 55d 1'14.26"N)
Band 1 Block=12925x1Type=UInt16, ColorInterp=Gray
> gdal_translate -a_srsepsg:32635 -gcp 0.5 0.5 498793 6.11075e+006 173.385 -gcp 8015.5 0.5 5073456.11027e+006 187.386 -gcp 8015.5 25599.5 507347 6.08351e+006 204.881 -gcp 0.5 25599.5 4987586.08385e+006 213.12D:\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. > copyD:\temp\rb\3V050909P0000897861A520004700712M_001631680\3v050909p0000897861a520004700712m_001631680_rpc.txt D:\temp\rb\3v050909p0000897861a520004700712m_001631680.org_rpc.txt
> gdalinfoD:\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 / UTMzone 35N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS84",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 StructureMetadata:
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=8016x1Type=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.tifD:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec.tif
RPC模型使用的坐标系是WGS84和数字高程数据（DEM），使用相对的大地水准面，这个值与真正的大地水准有一定的高差，而在正射校正是需要解决这种差异。如何确定这个高程差，一般使用的是采取图像的中心点的图像坐标（27D 2'550.34“Ë，55D 1'140.26”Ñ），并上传到这个在线资源，或者这个在线资源（此计算可以借助的网上资源，也可以使用的软件，如，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.tifD:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec2.tif
如果要将结果转换到另一个坐标系统，只需将t_srs参数这是为需要的坐标系即可。如果这时gdalwarp执行失败，很有可能就缺少DEM数据导致，下载相邻的DEM数据试试。
五、 结论
为了评估使用GDAL正射校正的进度，这里利用商业软件ENVI EX同GDAL进行对比。如下图，这里是一个图像在同一地点的，绿色的点是GPS轨迹。
在上图中，我们可以看到，当没有使用大地水准面进行正射校正的道路有些偏移。而使用大地水准面的高差进行正射的结果同时用软件ENVI EX的结果是相同的。
使用GDAL进行正射校正会出现下图中的横向锯齿问题，但是使用程序wxGIS处理的结果不会出现这样的情况，wxGIS也是基于GDAL库。为了消除这个问题，在命令行上，你需要添加选项-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.tifD:\temp\rb\3v050909p0000897861a520004700712m_001631680.rec2.tif
下面是另外一组GPS轨迹与处理结果叠加以及GDAL处理的横向锯齿问题。
六、数据下载