使用GDAL進行RPC座標轉換

使用GDAL進行RPC座標轉換

對於高解析度遙感衛星資料而言，目前幾乎都提供了有理函式模型（RFM）來進行影象校正（SPOT系列提供了有理函式模型之外還提供了嚴格軌道模型）。對遙感影像進行校正目前最常用的就是使用有理函式模型來進行校正。

有理函式模型的計算公式網上可以搜到，同時我之前的部落格中也有比較詳細的說明，可以參考這篇部落格http://blog.csdn.net/liminlu0314/article/details/24810593。

GDAL庫從1.3左右就開始提供基於RPC的座標轉換，到目前2.x版本以來，基於RPC的座標轉換一直在持續優化，從1.X版本到2.X版本，計算速度有了較大提升。

下面是使用GDAL基於RPC座標轉換的程式碼：

#include "gdal_priv.h"
#include "gdal_alg.h"

void main()
{
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
    CPLSetConfigOption("GDAL_DATA", "./data");
    GDALAllRegister();//註冊驅動

    const char* pszTif = "F:\\ZY3_TLC_E115.5_N39.8_20150110_L1A0002974893-NAD.tiff";

    //開啟資料
    GDALDataset *pDS = (GDALDataset*)GDALOpen(pszTif, GA_ReadOnly);
    //從元資料中獲取RPC資訊
    char** papszRPC = pDS->GetMetadata("RPC");

    //將獲取的RPC資訊構造成結構體
    GDALRPCInfo oInfo;
    GDALExtractRPCInfo(papszRPC, &oInfo);

    GDALClose((GDALDatasetH)pDS);

    //設定RPC模型中所需的DEM路徑
    char** papszTransOption = NULL;
    papszTransOption = CSLSetNameValue( papszTransOption, "RPC_DEM", "E:\\DEM.img" );   //設定DEM

    //使用RPC資訊，DEM等構造RPC轉換引數
    void *pRPCTransform = GDALCreateRPCTransformer(&oInfo, FALSE, 0, papszTransOption);

    //定義影象的四個角點行列號座標
    double dX[4] = {0, 24525, 24525, 0};
    double dY[4] = {0, 0, 24438, 24438};
    double dZ[4] = {0, 0, 0, 0};
    int nSuccess[4] = {FALSE};

    //輸出原始行列號座標
    printf("row col:\n");
    for(int i=0; i<4; i++)
        printf("%d\tX=%.7f\tY=%.7f\tZ=%.7f\n", i, dX[i], dY[i], dZ[i]);

    //呼叫RPC座標轉換函式進行座標轉換反算（行列號轉經緯度）
    GDALRPCTransform(pRPCTransform, FALSE, 4, dX, dY, dZ, nSuccess);

    //輸出轉換的經緯度
    printf("long lat:\n");
    for(int i=0; i<4; i++)
        printf("%d\tX=%.7f\tY=%.7f\tZ=%.7f\n", i, dX[i], dY[i], dZ[i]);

    //呼叫RPC座標轉換函式進行座標正算（經緯度轉行列號）
    GDALRPCTransform(pRPCTransform, TRUE, 4, dX, dY, dZ, nSuccess);


    //輸出轉換後的行列號座標
    printf("row2 col2:\n");
    for(int i=0; i<4; i++)
        printf("%d\tX=%.7f\tY=%.7f\tZ=%.7f\n", i, dX[i], dY[i], dZ[i]);

    //釋放資源
    GDALDestroyRPCTransformer(pRPCTransform);
    CSLDestroy(papszTransOption);
}