使用GDAL库中的RPC校正问题
最近将GDAL库更新至1.11版本之后,发现之前写的RPC像方改正模型校正的结果偏差特别大(更新版本之前结果和PCI处理的结果一致)。所以初步判断是GDAL库的bug,经过各个参数修改发现原来是指定的DEM采样方式导致的。
当指定DEM的采样方式为最邻近时,校正结果偏差很大,当DEM采样方式为双线性采样和三次立方卷积采样时,处理的结果与之前的结果一样。截图如图1所示,红色区域为对比区域,如图2所示。
图1 PCI校正结果全图
图2 图1中的红色区域按实际像素放大显示结果
下面是使用gdalwarp工具外加仿射修正模型进行校正的结果。首先使用DEM为最邻近采样,gdalwarp命令行如下:
gdalwarp -of GTiff -tr 2.5 2.5 -t_srs ESRI::"D:\WGS_1984_UTM_Zone_50N.prj" -rpc -r cubic -wm 1024 -to "RPC_AFFINE=-32.67376693 0.99919599 0.00013941 28.71988840 0.00062236 1.00004356" -to "RPC_DEM=D:\Data\正射纠正基础数据\DEM数据\beijing_all2.img" -to "RPC_DEMINTERPOLATION=near" D:\Data\711_214_26sep2006_p5\SrcData\bandf.tif D:\bandf_rpc1.tif --config GDAL_FILENAME_IS_UTF8 NO处理的结果与PCI结果对比如图3所示。
图3 GDAL指定DEM插值为最邻近处理结果
接下来指定DEM重采样方式为双线性,代码和处理结果如下:
gdalwarp -of GTiff -tr 2.5 2.5 -t_srs ESRI::"D:\WGS_1984_UTM_Zone_50N.prj" -rpc -r cubic -wm 1024 -to "RPC_AFFINE=-32.67376693 0.99919599 0.00013941 28.71988840 0.00062236 1.00004356" -to "RPC_DEM=D:\Data\正射纠正基础数据\DEM数据\beijing_all2.img" -to "RPC_DEMINTERPOLATION=bilinear" D:\Data\711_214_26sep2006_p5\SrcData\bandf.tif D:\bandf_rpc2.tif --config GDAL_FILENAME_IS_UTF8 NO处理的结果与PCI结果对比如图4所示。
图4 GDAL指定DEM插值为双线性处理结果
接下来指定DEM重采样方式为三次立方卷积,代码和处理结果如下:
gdalwarp -of GTiff -tr 2.5 2.5 -t_srs ESRI::"D:\WGS_1984_UTM_Zone_50N.prj" -rpc -r cubic -wm 1024 -to "RPC_AFFINE=-32.67376693 0.99919599 0.00013941 28.71988840 0.00062236 1.00004356" -to "RPC_DEM=D:\Data\正射纠正基础数据\DEM数据\beijing_all2.img" -to "RPC_DEMINTERPOLATION=cubic" D:\Data\711_214_26sep2006_p5\SrcData\bandf.tif D:\bandf_rpc3.tif --config GDAL_FILENAME_IS_UTF8 NO处理的结果与PCI结果对比如图5所示。
图5 GDAL指定DEM插值为三次立方卷积处理结果
修改GDAL源码中的文件gdal_rpc.cpp两处地方,就是将GDT_Int32修改为GDT_Float64。原因很简单就是变量dfDMEH是一个double类型,读取数据的时候也需要按照double类型进行读取,而原来的代码是按照Int32类型读取,导致读取出来的DEM高程值有问题。修改后代码见图6,共有两处,坐标正变换和逆变换都需要修改。
图6 修改后的GDAL代码
修改后重新编译GDAL即可。使用修改后的程序,指定DEM插值为最邻近后处理的结果如图7。
图7 修改后正射的结果对比
已经将该问题反馈至GDAL开发组,具体地址为:http://trac.osgeo.org/gdal/ticket/5553。估计下个版本就会修正。