欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

影像坐标重投影代码

程序员文章站 2022-04-01 10:06:25
...

源码如下,其中pInput为输入影像路径,pOut为输出影像路径,gcs为输出影像的地理坐标系

int DoReproject(string pInput, string pOut,string gcs)
{
    GDALAllRegister();
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
    GDALDataset* pDTIn = (GDALDataset*)GDALOpen(pInput.c_str(), GA_ReadOnly);
    if (pDTIn == NULL)
    {
        printf("文件打开失败: %s\n", pInput.c_str());
        return -1;
    }
    ///判断是否与目标坐标系统一直
    OGRSpatialReference obj_OSRS;
    obj_OSRS.SetWellKnownGeogCS(gcs.c_str());
    char *pszWKT = NULL;
    obj_OSRS.exportToWkt(&pszWKT);

    ///创建从UTM到84坐标系的坐标转换关系
    void *  hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDTIn, pDTIn->GetProjectionRef(),
        NULL, pszWKT, false, 0, 1);
    if (hTransformArg == NULL)
    {
        printf("无法创建转换方程!\n");
        GDALClose((GDALDatasetH)pDTIn);
        return -1;
    }

    //计算输出影像范围
    int nWidth = 0;
    int nHeight = 0;
    double dOutTransform[6];
    ///hTransformArg:坐标转换关系
    CPLErr eErr = GDALSuggestedWarpOutput((GDALDatasetH)pDTIn,
        GDALGenImgProjTransform, hTransformArg,
        dOutTransform, &nWidth, &nHeight);
    /*if (eErr == CE_None)
    {
        GDALDestroyGenImgProjTransformer(hTransformArg);
        return -1;
    }*/

    //创建输出影像
    GDALDataType eType = pDTIn->GetRasterBand(1)->GetRasterDataType();
    int nBandCnt = pDTIn->GetRasterCount();
    const char* format = "GTiff";
    GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName(format);
    GDALDataset* pDTReproj = (GDALDataset*)GDALCreate((GDALDriverH)pDriver,
        pOut.c_str(), nWidth, nHeight, nBandCnt, eType, NULL);
    if (pDTReproj == NULL)
    {
        printf("创建失败: %s\n", pOut.c_str());
        return -1;
    }
    pDTReproj->SetProjection(pszWKT);
    pDTReproj->SetGeoTransform(dOutTransform);

    //psWarpOptions set
    GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->hDstDS = pDTReproj;///输出图像指针
    psWarpOptions->hSrcDS = pDTIn;///原始图像指针
    psWarpOptions->nBandCount = 0;///处理的波段数据,0为处理所有波段
    psWarpOptions->eResampleAlg = GRA_Bilinear;///重采样方式
    psWarpOptions->panSrcBands = NULL;
    psWarpOptions->panDstBands = NULL;
    psWarpOptions->pfnProgress = GDALTermProgress;
    psWarpOptions->pTransformerArg =
        GDALCreateGenImgProjTransformer((GDALDatasetH)pDTIn,
        GDALGetProjectionRef(pDTIn),
        (GDALDatasetH)pDTReproj,
        GDALGetProjectionRef(pDTReproj),
        FALSE, 0.0, 1);   //创建从原始到输出结果的投影变换参数
    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;  //坐标变换函数指针


    //真正的重投影
    GDALWarpOperation oOperation;
    oOperation.Initialize(psWarpOptions);
    oOperation.ChunkAndWarpImage(0, 0, nWidth, nHeight);

    //销毁转换
    GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
    GDALDestroyWarpOptions(psWarpOptions);

    GDALClose(pDTIn);
    GDALClose(pDTReproj);

    return true;
}
相关标签: 影像处理 GDAL