影像坐标重投影代码
程序员文章站
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;
}
推荐阅读