pull/5/head
parent
3844933fac
commit
2bd17af3a6
|
|
@ -3571,84 +3571,126 @@ void MergeTiffs(QList<QString> inputFiles, QString outputFile) {
|
|||
|
||||
bool ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode)
|
||||
{
|
||||
// Initialize GDAL
|
||||
GDALAllRegister();
|
||||
|
||||
// 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;
|
||||
if (!poDataset) {
|
||||
std::cerr << "Failed to open input raster." << std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
// 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;
|
||||
// 创建目标空间参考
|
||||
OGRSpatialReference oSRS;
|
||||
if (oSRS.importFromEPSG(outepsgcode) != OGRERR_NONE) {
|
||||
std::cerr << "Failed to import EPSG code." << std::endl;
|
||||
GDALClose(poDataset);
|
||||
return false;
|
||||
}
|
||||
|
||||
GDALDataType datetype = poDataset->GetRasterBand(1)->GetRasterDataType();
|
||||
char* pszWKT = nullptr;
|
||||
if (oSRS.exportToWkt(&pszWKT) != OGRERR_NONE) {
|
||||
std::cerr << "Failed to export SRS to WKT." << std::endl;
|
||||
GDALClose(poDataset);
|
||||
return false;
|
||||
}
|
||||
|
||||
// Create the output raster dataset
|
||||
GDALDataset* poOutDataset = poDriver->Create(outRasterPath.toStdString().c_str(),
|
||||
poDataset->GetRasterXSize(),
|
||||
poDataset->GetRasterYSize(),
|
||||
poDataset->GetRasterCount(),
|
||||
poDataset->GetRasterBand(1)->GetRasterDataType(),
|
||||
nullptr);
|
||||
// 计算输出参数（尺寸和地理变换）
|
||||
void* hTransformArg = GDALCreateGenImgProjTransformer(poDataset, GDALGetProjectionRef(poDataset), nullptr, pszWKT, FALSE, 0, 1);
|
||||
if (!hTransformArg) {
|
||||
std::cerr << "Failed to create transformer for warp output calculation." << std::endl;
|
||||
CPLFree(pszWKT);
|
||||
GDALClose(poDataset);
|
||||
return false;
|
||||
}
|
||||
|
||||
// Create the target spatial reference using the EPSG code
|
||||
OGRSpatialReference oSRS;
|
||||
oSRS.importFromEPSG(outepsgcode);
|
||||
double adfGeoTransform[6];
|
||||
int nPixels = 0, nLines = 0;
|
||||
if (GDALSuggestedWarpOutput(poDataset, GDALGenImgProjTransform, hTransformArg, adfGeoTransform, &nPixels, &nLines) != CE_None) {
|
||||
std::cerr << "Failed to calculate output dimensions." << std::endl;
|
||||
GDALDestroyGenImgProjTransformer(hTransformArg);
|
||||
CPLFree(pszWKT);
|
||||
GDALClose(poDataset);
|
||||
return false;
|
||||
}
|
||||
GDALDestroyGenImgProjTransformer(hTransformArg);
|
||||
|
||||
char* pszWKT = NULL;
|
||||
oSRS.exportToWkt(&pszWKT);
|
||||
std::cout << "WKT of EPSG:" << outepsgcode << " :\n" << pszWKT << std::endl;
|
||||
// 创建输出数据集
|
||||
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
|
||||
if (!poDriver) {
|
||||
std::cerr << "GTiff driver not available." << std::endl;
|
||||
CPLFree(pszWKT);
|
||||
GDALClose(poDataset);
|
||||
return false;
|
||||
}
|
||||
|
||||
GDALDataType dataType = poDataset->GetRasterBand(1)->GetRasterDataType();
|
||||
GDALDataset* poOutDataset = poDriver->Create(outRasterPath.toStdString().c_str(), nPixels, nLines,
|
||||
poDataset->GetRasterCount(), dataType, nullptr);
|
||||
if (!poOutDataset) {
|
||||
std::cerr << "Failed to create output dataset." << std::endl;
|
||||
CPLFree(pszWKT);
|
||||
GDALClose(poDataset);
|
||||
return false;
|
||||
}
|
||||
|
||||
// 设置输出数据集的地理参考
|
||||
poOutDataset->SetGeoTransform(adfGeoTransform);
|
||||
poOutDataset->SetProjection(pszWKT);
|
||||
CPLFree(pszWKT);
|
||||
|
||||
// 配置GDAL Warp选项
|
||||
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;
|
||||
}
|
||||
|
||||
CPLSetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS");
|
||||
CPLSetConfigOption("GDAL_CACHEMAX", "16000");
|
||||
|
||||
psWarpOptions->eWorkingDataType = dataType;
|
||||
psWarpOptions->eResampleAlg = (dataType == GDT_Int16 || dataType == GDT_Int32 || dataType == GDT_Int64 || dataType == GDT_Int8 ||
|
||||
dataType == GDT_UInt16 || dataType == GDT_UInt32 || dataType == GDT_UInt64) ? GRA_NearestNeighbour : GRA_Bilinear;
|
||||
|
||||
psWarpOptions->hSrcDS = (GDALDatasetH)poDataset;
|
||||
psWarpOptions->hDstDS = (GDALDatasetH)poOutDataset;
|
||||
psWarpOptions->nBandCount = poDataset->GetRasterCount();
|
||||
psWarpOptions->pTransformerArg =
|
||||
GDALCreateGenImgProjTransformer(poDataset,
|
||||
GDALGetProjectionRef(poDataset),
|
||||
poOutDataset,
|
||||
GDALGetProjectionRef(poOutDataset),
|
||||
FALSE, 0.0, 1);
|
||||
|
||||
// 设置波段映射
|
||||
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;
|
||||
}
|
||||
|
||||
// 创建坐标转换器
|
||||
psWarpOptions->pTransformerArg = GDALCreateGenImgProjTransformer(poDataset, GDALGetProjectionRef(poDataset),
|
||||
poOutDataset, GDALGetProjectionRef(poOutDataset),
|
||||
FALSE, 0.0, 1);
|
||||
psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
|
||||
|
||||
|
||||
// 初始化并且执行变换操作
|
||||
// 执行坐标变换
|
||||
GDALWarpOperation oWo;
|
||||
if (oWo.Initialize(psWarpOptions) != CE_None) {
|
||||
GDALClose((GDALDatasetH)(GDALDatasetH)poDataset);
|
||||
GDALClose((GDALDatasetH)(GDALDatasetH)poOutDataset);
|
||||
std::cerr << "Warp initialization failed." << std::endl;
|
||||
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
|
||||
GDALDestroyWarpOptions(psWarpOptions);
|
||||
GDALClose(poDataset);
|
||||
GDALClose(poOutDataset);
|
||||
return false;
|
||||
}
|
||||
|
||||
oWo.ChunkAndWarpMulti(0, 0, poDataset->GetRasterYSize(), poDataset->GetRasterXSize());
|
||||
|
||||
if (oWo.ChunkAndWarpMulti(0, 0, nLines, nPixels) != CE_None) {
|
||||
std::cerr << "Warp operation failed." << std::endl;
|
||||
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
|
||||
GDALDestroyWarpOptions(psWarpOptions);
|
||||
GDALClose(poDataset);
|
||||
GDALClose(poOutDataset);
|
||||
return false;
|
||||
}
|
||||
|
||||
// 清理资源
|
||||
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
|
||||
GDALDestroyWarpOptions(psWarpOptions);
|
||||
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
|
||||
GDALClose(poDataset);
|
||||
GDALClose(poOutDataset);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue