Use RPC in gdal library to correct problems

After updating the gdal database to version 1.11 recently, it was found that the result deviation of the previously written RPC Image Correction Model was extremely large (the results before the version update were consistent with those of the PCI processing ). Therefore, it was initially determined that the gdal database had a bug. After each parameter modification, it was found that it was caused by the specified DEM sampling method.

When the sampling method of the specified DEM is the nearest neighbor, the deviation of the correction result is very large. When the DEM sampling method is bilinear sampling and cubic convolution sampling, the processing result is the same as the previous one. 1, the red area is the contrast area, and 2.

Figure 1 Overview of PCI correction results

Figure 2 The Red Area in Figure 1 is enlarged by actual pixels.

The following is the result of correction using the gdalwarp tool and the affinine correction model. First, use DEM as the nearest sampling. The gdalwarp command line is as follows:

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 "-to" rpc_dem = D: \ data \ normal shot correction basic data \ DEM Data \ 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

The processing result is shown in comparison with PCI result 3.

Figure 3 gdal specifies DEM Interpolation as the nearest processing result

Next, specify the dem sampling method to bilinear. The code and processing result are as follows:

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 "-to" rpc_dem = D: \ data \ normal shot correction basic data \ DEM Data \ 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

The processing result is shown in comparison with PCI Result 4.

Figure 4 gdal-specify DEM Interpolation as bilinear processing result

Next, specify the DEM sampling method as cubic convolution. The code and processing result are as follows:

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 "-to" rpc_dem = D: \ data \ normal shot correction basic data \ DEM Data \ 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

The processing result is shown in comparison with PCI result 5.

Figure 5 gdal-specify DEM Interpolation as cubic convolution processing result

Modify the gdal_rpc.cpp file in the gdal source code to change gdt_int32 to gdt_float64. The simple reason is that the variable dfdmeh is of the double type. When reading data, you also need to read data of the double type. The original Code reads data of the int32 type, the read DEM elevation value is incorrect. After modification, the code is shown in Figure 6. There are two parts in total. Both the coordinate positive transformation and the inverse transformation need to be modified.

Figure 6 modified gdal code

Modify and re-compile gdal. Use the modified program to specify DEM Interpolation as result 7 of the nearest post-processing.

Figure 7 after modification, the comparison of the positive shot result has reported the problem to the gdal development team. The specific address is http://trac.osgeo.org/gdal/ticket/5553. It is estimated that the next version will be corrected.