修复 ConvertCoordinateSystem 函数bug

pull/5/head
陈增辉 2025-02-21 15:19:00 +08:00
parent ff3259d0cd
commit 3844933fac
1 changed files with 70 additions and 45 deletions

View File

@ -3571,59 +3571,84 @@ void MergeTiffs(QList<QString> inputFiles, QString outputFile) {
bool ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode)
{
GDALAllRegister(); // 注册所有GDAL驱动
// Initialize GDAL
GDALAllRegister();
// 打开输入栅格
GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inRasterPath.toUtf8().constData(), GA_ReadOnly);
if (!poSrcDS) return false;
// 创建目标坐标系
OGRSpatialReference oTargetSRS;
oTargetSRS.importFromEPSG(outepsgcode);
char* pszTargetSRS = nullptr;
oTargetSRS.exportToWkt(&pszTargetSRS);
// 设置Warp选项
GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->hSrcDS = poSrcDS;
psWarpOptions->hDstDS = nullptr;
psWarpOptions->nBandCount = poSrcDS->GetRasterCount();
psWarpOptions->panSrcBands = (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
psWarpOptions->panDstBands = (int*)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
for (int i = 0; i < psWarpOptions->nBandCount; ++i) {
psWarpOptions->panSrcBands[i] = i + 1;
psWarpOptions->panDstBands[i] = i + 1;
// Open the input raster file
GDALDataset* poDataset = (GDALDataset*)GDALOpen(inRasterPath.toStdString().c_str(), GA_ReadOnly);
if (poDataset == nullptr) {
std::cerr << "Failed to open input raster file." << std::endl;
return false;
}
// 设置坐标转换参数
psWarpOptions->papszWarpOptions = CSLSetNameValue(nullptr, "INIT_DEST", "NO_DATA");
psWarpOptions->dfWarpMemoryLimit = 8000; // 内存限制8000MB
// Create the output raster driver
const char* pszDriver = "GTiff"; // You can change this depending on your desired output format
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName(pszDriver);
if (poDriver == nullptr) {
std::cerr << "Failed to get output driver." << std::endl;
GDALClose(poDataset);
return false;
}
// 创建坐标转换器
GDALDataType datetype = poDataset->GetRasterBand(1)->GetRasterDataType();
// Create the output raster dataset
GDALDataset* poOutDataset = poDriver->Create(outRasterPath.toStdString().c_str(),
poDataset->GetRasterXSize(),
poDataset->GetRasterYSize(),
poDataset->GetRasterCount(),
poDataset->GetRasterBand(1)->GetRasterDataType(),
nullptr);
// Create the target spatial reference using the EPSG code
OGRSpatialReference oSRS;
oSRS.importFromEPSG(outepsgcode);
char* pszWKT = NULL;
oSRS.exportToWkt(&pszWKT);
std::cout << "WKT of EPSG:" << outepsgcode << " :\n" << pszWKT << std::endl;
poOutDataset->SetProjection(pszWKT);
GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
CPLSetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); // 使用所有可用的CPU核心
CPLSetConfigOption("GDAL_CACHEMAX", "16000"); // 设置缓存大小为500MB
psWarpOptions->eWorkingDataType = datetype;
if (datetype == GDT_Int16 || datetype == GDT_Int32 || datetype == GDT_Int64 || datetype == GDT_Int8 ||
datetype == GDT_UInt16 || datetype == GDT_UInt32 || datetype == GDT_UInt64) {
psWarpOptions->eResampleAlg = GRA_NearestNeighbour;
}
else {
psWarpOptions->eResampleAlg = GRA_Bilinear;
}
psWarpOptions->hSrcDS = (GDALDatasetH)poDataset;
psWarpOptions->hDstDS = (GDALDatasetH)poOutDataset;
psWarpOptions->nBandCount = poDataset->GetRasterCount();
psWarpOptions->pTransformerArg =
GDALCreateGenImgProjTransformer(poSrcDS, nullptr, nullptr, pszTargetSRS,
FALSE, 0, 1);
GDALCreateGenImgProjTransformer(poDataset,
GDALGetProjectionRef(poDataset),
poOutDataset,
GDALGetProjectionRef(poOutDataset),
FALSE, 0.0, 1);
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
// 创建输出文件
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = poDriver->Create(
outRasterPath.toUtf8().constData(),
0, 0, 0, GDT_Unknown, nullptr);
// 执行坐标转换
GDALWarpOperation oOperation;
oOperation.Initialize(psWarpOptions);
oOperation.ChunkAndWarpImage(0, 0,
poSrcDS->GetRasterXSize(),
poSrcDS->GetRasterYSize());
// 清理资源
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
GDALDestroyWarpOptions(psWarpOptions);
CPLFree(pszTargetSRS);
GDALClose(poSrcDS);
GDALClose(poDstDS);
// 初始化并且执行变换操作
GDALWarpOperation oWo;
if (oWo.Initialize(psWarpOptions) != CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)poDataset);
GDALClose((GDALDatasetH)(GDALDatasetH)poOutDataset);
return false;
}
oWo.ChunkAndWarpMulti(0, 0, poDataset->GetRasterYSize(), poDataset->GetRasterXSize());
GDALFlushCache(poOutDataset);
qDebug() << "ChunkAndWarpImage over" << Qt::endl;
// GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
// GDALDestroyWarpOptions(psWo);
GDALClose((GDALDatasetH)(GDALDatasetH)poDataset);
GDALClose((GDALDatasetH)(GDALDatasetH)poOutDataset);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return true;
}