使用GDAL工具对OrbView3数据进行正射校正

本文原文地址：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.tif

Windows系统：

>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

第二个方法是指定图像的坐标系和四个角点的坐标。查看图像的四个角点坐标，可以在文件scene.pvl中进行查找。然后使用gdal_translate工具进行转换处理，命令如下：

> 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

然后再使用gdalinfo查看输出文件中的信息：

> 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处理的横向锯齿问题。

六、数据下载

原始数据3v050909p0000897861a520004700712m_001631680可以通过EarthSat（下载）
ASTER GDEM数据用于正射校正（下载）
GPS点轨迹数据（GPX）检查结果（下载和下载）
ENVI EX正射结果（下载）
GDAL的正射纠正结果（下载）
GDAL调整大地水准面后的正射纠正结果（下载）