使用GDAL工具对OrbView-3数据进行正射校正

最新推荐文章于 2021-11-07 11:46:05 发布

langjitianya43

最新推荐文章于 2021-11-07 11:46:05 发布

阅读量1.1k

点赞数

分类专栏： GDAL

GDAL 专栏收录该内容

12 篇文章 0 订阅

订阅专栏

本文原文地址：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）是白俄罗斯共和国和库尔斯克地区（点击各自的链接会到各自的示例数据目录）。

库尔斯克地区的数据集，包括：

文件名

说明

3V050401P0000692211A520018202082M_001639429.ZIP

OrbView-3w卫星数据包

3v050401p0000692211a520018202082m_001639429_ortho.img

ERDAS IMAGINE的正射校正结果的格式

3v050401p0000692211a520018202082m_001639429_ortho.img.aux.xml

辅助文件

3v050401p0000692211a520018202082m_001639429_ortho.rrd

金字塔文件

dem.tfw

DEM坐标

dem.tif

DEM

dem.tif.aux.xml

DEM辅助数据

dem.tif.ovr

DEM金字塔

dem.tif.xml

DEM辅助数据

kursk_ortho.7z

使用GDAL正射校正结果

kursk_ortho_envi_orbview.7z

使用ENVI正射校正结果

白俄罗斯共和国的数据集，包括：

文件名

描述

3v050909p0000897861a520004700712m_001631680.zip

OrbView-3w卫星数据包

aster_gdem.tif

Aster的DEM数据

ov3-check.7z

地面实测的GPS点

result_envi.tif

ENVI-EX正射结果

result_gdal.7z

GDAL正射结果

result_gdal_geoid_corrected.7z

经过大地水准面修正的GDAL正射结果

tracks.7z

地面实测GPS点轨迹

三、准备正射校正必要的数据

在这篇文章中，所有的例子中使用的数据是上面白俄罗斯共和国一组数据。

首先，考虑OrbView-3的数据是ZIP压缩格式（3V050909P0000897861A520004700712M_001631680.ZIP）。所以先要对其进行解压：

UNIX系统：

[plain] view plaincopyprint?
$ ls -13V050909P0000897861A520004700712M_001631680*
Windows系统：
[plain] view plaincopyprint?

dir * /B
执行完上面的命令，会输出：

[plain] view plaincopyprint?
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系统：

[plain] view plaincopyprint?
$ gdalinfo 3v050909p0000897861a520004700712m_001631680.tif
Windows系统：
[plain] view plaincopyprint?

C:\warmerda\bld\bin>gdalinfo.exe3v050909p0000897861a520004700712m_001631680.tif
执行完上面的命令，会输出：

[plain] view plaincopyprint?
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下自己肯定会解压吧。
[plain] view plaincopyprint?
for i in srtm*zip; doyes|unzip $i; done

将所有文件SRTM的GeoTIFF格式转换成一个单一的虚拟文件。使用的工具是gdalbuildvrt。命令是（第一行是UNIX命令，第二行是Windows命令）：
[plain] view plaincopyprint?
$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格式的文件，命令如下：
[plain] view plaincopyprint?
$ 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卫星影像进行正射校正使用下面的命令：

[plain] view plaincopyprint?

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校正之后，其余均使用默认参数，如下：

[plain] view plaincopyprint?
$ 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查看输出的结果，命令以及输出信息如下：

[plain] view plaincopyprint?
$ 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工具进行转换处理，命令如下：

[plain] view plaincopyprint?

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查看输出文件中的信息：

[plain] view plaincopyprint?

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

现在，您可以重复执行该正射校正命令，得到新的图像文件。
[plain] view plaincopyprint?

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模型上。
对于使用消除大地水准面高差进行正射校正的命令如下：

[plain] view plaincopyprint?

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。示例命令：

[plain] view plaincopyprint?

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调整大地水准面后的正射纠正结果（下载）

转载自：http://blog.csdn.net/liminlu0314/article/details/8193362