gdalwarp use rpc & dem bug

[liminlu0314]

Use the following command to run the error,

gdalwarp.exe -rpc E:\GDAL\bug\GF1_PMS1_E116.5_N39.4_20131127_L1A0000117600_PAN.vrt E:\GDAL\bug\out.tif -to RPC_DEM=E:\GDAL\bug\dem.tif -to RPC_DEMINTERPOLATION=bilinear --debug on -t_srs E:\GDAL\bug\prj.txt

The output message is as follows:

...
RPC: DEM (pixel, line) = (-3.23331, 794.046)
RPC: Iteration 1 for (pixel, line) = (727.68, 18000): No elevation value at 116.289614568165 39.2970459265271. Erroring out
RPC: DEM (pixel, line) = (-5.84626, 793.942)
RPC: Iteration 0 for (pixel, line) = (727.68, 18000): No elevation value at 116.287437116176 39.2971332931275. Using elevation 18.8833 at DEM (pixel, line) = (0.5, 793.942) (snapping to boundaries) instead
RPC: DEM (pixel, line) = (-3.23331, 794.046)
RPC: Iteration 1 for (pixel, line) = (727.68, 18000): No elevation value at 116.289614568023 39.2970459265535. Erroring out
GDAL: GDALDriver::Create(GTiff,E:\GDAL\bug\out.tif,37409837,923358668,1,UInt16,0000000000000000)
ERROR 3: Free disk space available is 90299555840 bytes, whereas 69085394524834232 are at least necessary. You can disable this check by defining the CHECK_DISK_FREE_SPACE configuration option to FALSE.
GDAL: GDALClose(E:\GDAL\bug\dem.tif, this=00000193B82511E0)
GDAL: GDALClose(E:\GDAL\bug\GF1_PMS1_E116.5_N39.4_20131127_L1A0000117600_PAN.vrt, this=00000193B82558A0)

The main reason is that in the RPC coordinate transformation, part of the coordinate transformation fails, and then in reprojection (see the code below), the change failed flag is covered to true, resulting in the output image size error.

​https://trac.osgeo.org/gdal/browser/trunk/gdal/alg/gdaltransformer.cpp#L1973

/* -------------------------------------------------------------------- */
/*      Reproject if needed.                                            */
/* -------------------------------------------------------------------- */
    if( psInfo->pReprojectArg )
    {
        if( !psInfo->pReproject( psInfo->pReprojectArg, bDstToSrc,
                                 nPointCount, padfX, padfY, padfZ,
                                 panSuccess ) )
            return FALSE;
    }

​https://trac.osgeo.org/gdal/browser/trunk/gdal/ogr/ogrct.cpp#L1271

/* -------------------------------------------------------------------- */
/*      Establish error information if pabSuccess provided.             */
/* -------------------------------------------------------------------- */
    if( pabSuccess )
    {
        for( int i = 0; i < nCount; i++ )
        {
            if( x[i] == HUGE_VAL || y[i] == HUGE_VAL )
                pabSuccess[i] = FALSE;
            else
                pabSuccess[i] = TRUE;
        }
    }

The attachment is the data used. Since the original tif data is too large to be committed, it is converted to the VRT format.

[Even Rouault]

In 40446:

RPC transformer: set output coordinates to HUGE_VAL when failure occurs, so that a following coordinate transformation can detect the error too (analysis done by liminlu0314, fixes #7090)

[Even Rouault]

In 40447:

RPC transformer: set output coordinates to HUGE_VAL when failure occurs, so that a following coordinate transformation can detect the error too (analysis done by liminlu0314, fixes #7090)