修复坐标系转换的bug

pull/5/head
陈增辉 2025-02-23 11:04:22 +08:00
parent 2bd17af3a6
commit 3292dc63e9
4 changed files with 79 additions and 120 deletions

View File

@ -3566,139 +3566,83 @@ void MergeTiffs(QList<QString> inputFiles, QString outputFile) {
GDALClose(poDstDS);
GDALClose(poFirstDS);
}
bool ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode)
{
void ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode) {
// 注册所有GDAL驱动
GDALAllRegister();
GDALDataset* poDataset = (GDALDataset*)GDALOpen(inRasterPath.toStdString().c_str(), GA_ReadOnly);
if (!poDataset) {
std::cerr << "Failed to open input raster." << std::endl;
return false;
// 打开输入栅格文件
GDALDataset* srcDataset = (GDALDataset*)GDALOpen(inRasterPath.toUtf8().constData(), GA_ReadOnly);
if (!srcDataset) {
// 错误处理:输出文件打开失败
// qDebug() << "无法打开输入文件:" << inRasterPath;
return;
}
// 创建目标空间参考
OGRSpatialReference oSRS;
if (oSRS.importFromEPSG(outepsgcode) != OGRERR_NONE) {
std::cerr << "Failed to import EPSG code." << std::endl;
GDALClose(poDataset);
return false;
// 创建目标坐标系
OGRSpatialReference targetSRS;
if (targetSRS.importFromEPSG(outepsgcode) != OGRERR_NONE) {
GDALClose(srcDataset);
// qDebug() << "无效的EPSG代码" << outepsgcode;
return;
}
char* pszWKT = nullptr;
if (oSRS.exportToWkt(&pszWKT) != OGRERR_NONE) {
std::cerr << "Failed to export SRS to WKT." << std::endl;
GDALClose(poDataset);
return false;
// 获取目标坐标系的WKT表示
char* targetSRSWkt = nullptr;
targetSRS.exportToWkt(&targetSRSWkt);
// 创建重投影后的虚拟数据集Warped VRT
GDALDataset* warpedVRT = (GDALDataset*)GDALAutoCreateWarpedVRT(
srcDataset,
nullptr, // 输入坐标系(默认使用源数据)
targetSRSWkt, // 目标坐标系
GRA_Bilinear, // 重采样方法:双线性插值
0.0, // 最大误差0表示自动计算
nullptr // 其他选项
);
CPLFree(targetSRSWkt); // 释放WKT内存
if (!warpedVRT) {
GDALClose(srcDataset);
// qDebug() << "创建投影转换VRT失败";
return;
}
// 计算输出参数(尺寸和地理变换)
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;
// 获取输出驱动GeoTIFF格式
GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (!driver) {
GDALClose(warpedVRT);
GDALClose(srcDataset);
// qDebug() << "无法获取GeoTIFF驱动";
return;
}
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);
// 创建输出栅格文件
GDALDataset* dstDataset = driver->CreateCopy(
outRasterPath.toUtf8().constData(), // 输出文件路径
warpedVRT, // 输入数据集VRT
false, // 是否严格复制
nullptr, // 创建选项
nullptr, // 进度回调
nullptr // 回调参数
);
// 创建输出数据集
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (!poDriver) {
std::cerr << "GTiff driver not available." << std::endl;
CPLFree(pszWKT);
GDALClose(poDataset);
return false;
if (!dstDataset) {
// qDebug() << "创建输出文件失败:" << outRasterPath;
GDALClose(warpedVRT);
GDALClose(srcDataset);
return;
}
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");
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->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) {
std::cerr << "Warp initialization failed." << std::endl;
GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
GDALDestroyWarpOptions(psWarpOptions);
GDALClose(poDataset);
GDALClose(poOutDataset);
return false;
}
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);
GDALClose(poDataset);
GDALClose(poOutDataset);
return true;
// 释放资源
GDALClose(dstDataset);
GDALClose(warpedVRT);
GDALClose(srcDataset);
}
void ResampleByReferenceRasterB(QString rasterAPath, QString rasterBPath, QString rasterCPath) {
GDALAllRegister(); // 初始化GDAL驱动

View File

@ -270,7 +270,7 @@ int BASECONSTVARIABLEAPI saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString o
void BASECONSTVARIABLEAPI clipRaster(QString inRasterPath, QString outRasterPath, long minRow, long maxRow, long minCol, long maxCol);
// 坐标系转换
bool BASECONSTVARIABLEAPI ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode);
void BASECONSTVARIABLEAPI ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode);
//--------------------- 图像合并流程 ------------------------------
enum MERGEMODE

View File

@ -49,7 +49,7 @@
</attribute>
<layout class="QVBoxLayout" name="verticalLayout_2">
<item>
<layout class="QGridLayout" name="infoLayout"/>
<layout class="QGridLayout" name="mapCanvasLayout"/>
</item>
</layout>
</widget>
@ -85,7 +85,7 @@
<x>0</x>
<y>0</y>
<width>906</width>
<height>23</height>
<height>22</height>
</rect>
</property>
<widget class="QMenu" name="projectMenu">
@ -164,6 +164,9 @@
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="windowTitle">
<string>图层列表</string>
</property>
<attribute name="dockWidgetArea">
<number>1</number>
</attribute>
@ -181,7 +184,7 @@
</sizepolicy>
</property>
<property name="title">
<string>图层列表</string>
<string/>
</property>
<layout class="QHBoxLayout" name="horizontalLayout">
<item>
@ -234,6 +237,9 @@ p, li { white-space: pre-wrap; }
</widget>
</widget>
<widget class="QDockWidget" name="dockWidget_2">
<property name="windowTitle">
<string>下载任务</string>
</property>
<attribute name="dockWidgetArea">
<number>2</number>
</attribute>
@ -357,13 +363,16 @@ p, li { white-space: pre-wrap; }
</widget>
</widget>
<widget class="QDockWidget" name="dockWidget_3">
<property name="windowTitle">
<string>信息窗口</string>
</property>
<attribute name="dockWidgetArea">
<number>1</number>
</attribute>
<widget class="QWidget" name="dockWidgetContents_4">
<layout class="QVBoxLayout" name="verticalLayout_5">
<item>
<layout class="QGridLayout" name="mapCanvasLayout"/>
<layout class="QGridLayout" name="infoLayout"/>
</item>
</layout>
</widget>

View File

@ -50,6 +50,9 @@
<height>30</height>
</size>
</property>
<property name="text">
<string>C:\Users\30453\Desktop\RasterTool\LT1A_S1GBM_20241229\LT1A_S1GBM_20241229WGS84.tif</string>
</property>
</widget>
</item>
<item row="1" column="2">
@ -73,6 +76,9 @@
<height>30</height>
</size>
</property>
<property name="text">
<string>C:\Users\30453\Desktop\RasterTool\LT1A_S1GBM_20241229\LT1A_S1GBM_20241229.tif</string>
</property>
</widget>
</item>
<item row="0" column="0">