gdal_rpc.cpp GDALCreateRPCTransformer incorrect linear approximator

[jaredrubinaz]

When I used the corner coordinates returned by gdal for an image with RPC coefficients and reproject from ground to image coordinates, the image coordinates are way off. Based upon walking through the simple linear approximator in GDALCreateRPCTranformer method, it looks like values [2] and [4] of adfGTFromLL are swapped. When I reproject the corner coordinates they are within 15-20 image coordinates..which I would expect.

In addition, when building the transform, the psRPCInfo->dfHEIGHT_OFF should be used for elevation instead of 0.0.

So the following update is needed, in addition to using psRPCInfo->dfHEIGHT_OF instead of 0.0

adfGTFromLL[2] = (dfRefLineDelta - dfRefLine) / dfLLDelta;

should be adfGTFromLL[4] = (dfRefLineDelta - dfRefLine) / dfLLDelta;

and

adfGTFromLL[4] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta; should be adfGTFromLL[2] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;

thanks Jared

Below is how I reimplemented the code block

void *GDALCreateRPCTransformer( GDALRPCInfo *psRPCInfo, int bReversed, 
                                double dfPixErrThreshold )

{
    GDALRPCTransformInfo *psTransform;

/* -------------------------------------------------------------------- */
/*      Initialize core info.                                           */
/* -------------------------------------------------------------------- */
    psTransform = (GDALRPCTransformInfo *) 
        CPLCalloc(sizeof(GDALRPCTransformInfo),1);

    memcpy( &(psTransform->sRPC), psRPCInfo, sizeof(GDALRPCInfo) );
    psTransform->bReversed = bReversed;
    psTransform->dfPixErrThreshold = dfPixErrThreshold;

    strcpy( psTransform->sTI.szSignature, "GTI" );
    psTransform->sTI.pszClassName = "GDALRPCTransformer";
    psTransform->sTI.pfnTransform = GDALRPCTransform;
    psTransform->sTI.pfnCleanup = GDALDestroyRPCTransformer;
    psTransform->sTI.pfnSerialize = NULL;

/* -------------------------------------------------------------------- */
/*      Establish a reference point for calcualating an affine          */
/*      geotransform approximate transformation.                        */
/* -------------------------------------------------------------------- */
    double adfGTFromLL[6], dfRefPixel, dfRefLine;

    double dfRefLong = (psRPCInfo->dfMIN_LONG + psRPCInfo->dfMAX_LONG) * 0.5;
    double dfRefLat  = (psRPCInfo->dfMIN_LAT  + psRPCInfo->dfMAX_LAT ) * 0.5;

	 //    RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat, 0.0, 
	 //                       &dfRefPixel, &dfRefLine );
	
	 // AFA - Use the height offset instead
    RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat, psRPCInfo->dfHEIGHT_OFF, 
                       &dfRefPixel, &dfRefLine );

/* -------------------------------------------------------------------- */
/*      Transform nearby locations to establish affine direction        */
/*      vectors.                                                        */
/* -------------------------------------------------------------------- */
    double dfRefPixelDelta, dfRefLineDelta, dfLLDelta = 0.0001;
    
	 //  RPCTransformPoint( psRPCInfo, dfRefLong+dfLLDelta, dfRefLat, 0.0, 
	 //                       &dfRefPixelDelta, &dfRefLineDelta );

	 // AFA - Use the height offset instead
    RPCTransformPoint( psRPCInfo, dfRefLong+dfLLDelta, dfRefLat, psRPCInfo->dfHEIGHT_OFF, 
                       &dfRefPixelDelta, &dfRefLineDelta );

    adfGTFromLL[1] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
	 //    adfGTFromLL[2] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
	 // AFA Correction
    adfGTFromLL[4] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
    
	 //    RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat+dfLLDelta, 0.0, 
	 //                       &dfRefPixelDelta, &dfRefLineDelta );

	 // AFA - Use the height offset instead
    RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat+dfLLDelta, psRPCInfo->dfHEIGHT_OFF, 
                       &dfRefPixelDelta, &dfRefLineDelta );

	 //    adfGTFromLL[4] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
	 // AFA Correction
    adfGTFromLL[2] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
    adfGTFromLL[5] = (dfRefLineDelta - dfRefLine) / dfLLDelta;

    adfGTFromLL[0] = dfRefPixel 
        - adfGTFromLL[1] * dfRefLong - adfGTFromLL[2] * dfRefLat;
    adfGTFromLL[3] = dfRefLine 
        - adfGTFromLL[4] * dfRefLong - adfGTFromLL[5] * dfRefLat;

    GDALInvGeoTransform( adfGTFromLL, psTransform->adfPLToLatLongGeoTransform);
    
    return psTransform;
}

Jared

[Even Rouault]

Hum, it looks like it has just been independantly fixed by r26963, #5395