修复bug

pull/5/head
chenzenghui 2025-02-24 11:41:41 +08:00
parent 3292dc63e9
commit 32be3d8f48
4 changed files with 159 additions and 88 deletions

View File

@ -3606,12 +3606,14 @@ void ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long o
if (!warpedVRT) {
GDALClose(srcDataset);
// qDebug() << "创建投影转换VRT失败";
qDebug() << u8"创建投影转换VRT失败";
return;
}
// 获取输出驱动GeoTIFF格式
GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
QString filesuffer = getFileExtension(outRasterPath).toLower();
bool isTiff = filesuffer.contains("tif");
GDALDriver* driver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
if (!driver) {
GDALClose(warpedVRT);
GDALClose(srcDataset);
@ -3643,103 +3645,154 @@ void ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long o
}
void ResampleByReferenceRasterB(QString rasterAPath, QString rasterBPath, QString rasterCPath) {
GDALAllRegister(); // 初始化GDAL驱动
void ResampleByReferenceRasterB(QString pszSrcFile, QString RefrasterBPath, QString pszOutFile, GDALResampleAlg eResample) {
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile.toLocal8Bit().constData(), GA_ReadOnly);
if (pDSrc == NULL) {
qDebug() << u8"do not open In Raster file: " << pszSrcFile;
return ;
}
// 阶段1读取参考栅格B的元数据
GDALDatasetH hRefDS = GDALOpen(rasterBPath.toUtf8().constData(), GA_ReadOnly);
if (!hRefDS) {
qDebug() << u8"参考栅格打开失败:" << rasterBPath;
GDALDataset* pDRef = (GDALDataset*)GDALOpen(RefrasterBPath.toLocal8Bit().constData(), GA_ReadOnly);
if (pDRef == NULL) {
qDebug() << u8"do not open Ref Raster file: " << RefrasterBPath;
return;
}
// 获取地理变换参数和尺寸
double geotransform[6];
int refWidth = GDALGetRasterXSize(hRefDS);
int refHeight = GDALGetRasterYSize(hRefDS);
if (GDALGetGeoTransform(hRefDS, geotransform) != CE_None) {
GDALClose(hRefDS);
qDebug() << u8"参考栅格缺少地理变换参数:" << rasterBPath;
return;
QString filesuffer = getFileExtension(pszOutFile).toLower();
bool isTiff = filesuffer.contains("tif");
GDALDriver* pDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
if (pDriver == NULL) {
qDebug() << "not open driver";
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return ;
}
int width = pDSrc->GetRasterXSize();
int height = pDSrc->GetRasterYSize();
int nBandCount = pDSrc->GetRasterCount();
GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();
char* pszSrcWKT = NULL;
pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef());
// 锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
if (strlen(pszSrcWKT) <= 0) {
OGRSpatialReference oSRS;
oSRS.importFromEPSG(4326);
// oSRS.SetUTM(50, true); //锟斤拷锟斤拷锟斤拷 锟斤拷锟斤拷120锟斤拷
// oSRS.SetWellKnownGeogCS("WGS84");
oSRS.exportToWkt(&pszSrcWKT);
}
// 计算目标地理范围
const double xmin = geotransform[0];
const double ymax = geotransform[3];
const double xmax = xmin + geotransform[1] * refWidth;
const double ymin = ymax + geotransform[5] * refHeight;
char* pdstSrcWKT = NULL;
pdstSrcWKT = const_cast<char*>(pDRef->GetProjectionRef());
// 获取目标坐标系
const char* projRef = GDALGetProjectionRef(hRefDS);
OGRSpatialReference targetSRS;
if (projRef[0] == '\0' || targetSRS.importFromWkt(projRef) != OGRERR_NONE) {
GDALClose(hRefDS);
qDebug() << u8"参考栅格坐标系无效:" << rasterBPath;
return;
}
GDALClose(hRefDS); // 尽早释放参考数据集
// 阶段2配置Warp参数
char** warpOptions = nullptr;
warpOptions = CSLAddString(warpOptions, "-t_srs"); // 目标坐标系
char* targetSRSWkt = nullptr;
targetSRS.exportToWkt(&targetSRSWkt);
warpOptions = CSLAddString(warpOptions, targetSRSWkt);
warpOptions = CSLAddString(warpOptions, "-te"); // 目标范围
warpOptions = CSLAddString(warpOptions, QString::number(xmin, 'f', 10).toUtf8().constData());
warpOptions = CSLAddString(warpOptions, QString::number(ymin, 'f', 10).toUtf8().constData());
warpOptions = CSLAddString(warpOptions, QString::number(xmax, 'f', 10).toUtf8().constData());
warpOptions = CSLAddString(warpOptions, QString::number(ymax, 'f', 10).toUtf8().constData());
warpOptions = CSLAddString(warpOptions, "-ts"); // 目标尺寸
warpOptions = CSLAddString(warpOptions, QString::number(refWidth).toUtf8().constData());
warpOptions = CSLAddString(warpOptions, QString::number(refHeight).toUtf8().constData());
warpOptions = CSLAddString(warpOptions, "-r"); // 双线性重采样
warpOptions = CSLAddString(warpOptions, "bilinear");
warpOptions = CSLAddString(warpOptions, "-of"); // 输出为GTiff
warpOptions = CSLAddString(warpOptions, "GTiff");
warpOptions = CSLAddString(warpOptions, "-wo"); // 启用多线程
warpOptions = CSLAddString(warpOptions, "NUM_THREADS=ALL_CPUS");
warpOptions = CSLAddString(warpOptions, "-co"); // 压缩优化
warpOptions = CSLAddString(warpOptions, "COMPRESS=LZW");
// 阶段3执行重采样操作
GDALWarpAppOptions* warpAppOptions = GDALWarpAppOptionsNew(warpOptions, nullptr);
GDALDatasetH hSrcDS = GDALOpen(rasterAPath.toUtf8().constData(), GA_ReadOnly);
if (!hSrcDS) {
GDALWarpAppOptionsFree(warpAppOptions);
CSLDestroy(warpOptions);
CPLFree(targetSRSWkt);
qDebug() << u8"源栅格打开失败:" << rasterAPath;
return;
// 锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
if (strlen(pdstSrcWKT) <= 0)
{
OGRSpatialReference oSRS;
oSRS.importFromEPSG(4326);
// oSRS.SetUTM(50, true); //锟斤拷锟斤拷锟斤拷 锟斤拷锟斤拷120锟斤拷
// oSRS.SetWellKnownGeogCS("WGS84");
oSRS.exportToWkt(&pdstSrcWKT);
}
int bUsageError = FALSE;
GDALDatasetH hDstDS = GDALWarp(
rasterCPath.toUtf8().constData(), // 输出路径
nullptr, // 自动选择驱动
1, &hSrcDS, // 输入数据集
warpAppOptions, // 转换参数
&bUsageError // 错误状态
);
// 阶段4资源清理与状态检查
GDALClose(hSrcDS);
if (hDstDS) GDALClose(hDstDS);
GDALWarpAppOptionsFree(warpAppOptions);
CSLDestroy(warpOptions);
CPLFree(targetSRSWkt);
int new_width = pDRef->GetRasterXSize();
int new_height = pDRef->GetRasterYSize();
double gt[6] ;
pDRef->GetGeoTransform(gt);
if (bUsageError || !hDstDS) {
qDebug() << u8"重采样过程发生错误";
// GDALDestroyGenImgProjTransformer(hTransformArg);
qDebug() << "create init GDALDataset ";
GDALDataset* pDDst =
pDriver->Create(pszOutFile.toLocal8Bit().constData(), new_width, new_height, nBandCount, dataType, NULL);
if (pDDst == NULL) {
qDebug() << "not create init GDALDataset ";
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
return ;
}
pDDst->SetProjection(pdstSrcWKT);
pDDst->SetGeoTransform(gt);
qDebug() << "GDALCreateGenImgProjTransformer " << Qt::endl;
void* hTransformArg;
hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT,
FALSE, 0.0, 1);
qDebug() << "no proj ";
//(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7)
if (hTransformArg == NULL) {
qDebug() << "hTransformArg create failure";
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
return;
}
qDebug() << u8"重采样成功完成:" << rasterCPath;
qDebug() << "has proj ";
double dGeoTrans[6] = { 0 };
int nNewWidth = 0, nNewHeight = 0;
if (GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg,
dGeoTrans, &nNewWidth, &nNewHeight)
!= CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
return ;
}
GDALWarpOptions* psWo = GDALCreateWarpOptions();
CPLSetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); // 使用所有可用的CPU核心
CPLSetConfigOption("GDAL_CACHEMAX", "16000"); // 设置缓存大小为500MB
// psWo->papszWarpOptions = CSLDuplicate(NULL);
psWo->eWorkingDataType = dataType;
psWo->eResampleAlg = eResample;
psWo->hSrcDS = (GDALDatasetH)pDSrc;
psWo->hDstDS = (GDALDatasetH)pDDst;
qDebug() << "GDALCreateGenImgProjTransformer" ;
psWo->pfnTransformer = GDALGenImgProjTransform;
psWo->pTransformerArg = GDALCreateGenImgProjTransformer(
(GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);
;
qDebug() << "GDALCreateGenImgProjTransformer has created" << Qt::endl;
psWo->nBandCount = nBandCount;
psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));
for (int i = 0; i < nBandCount; i++) {
psWo->panSrcBands[i] = i + 1;
psWo->panDstBands[i] = i + 1;
}
GDALWarpOperation oWo;
if (oWo.Initialize(psWo) != CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
return ;
}
qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << Qt::endl;
oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
GDALFlushCache(pDDst);
qDebug() << "ChunkAndWarpImage over" << Qt::endl;
// GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
// GDALDestroyWarpOptions(psWo);
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return ;
}

View File

@ -261,6 +261,9 @@ ErrorCode BASECONSTVARIABLEAPI transformCoordinate(double x, double y, int sou
int BASECONSTVARIABLEAPI alignRaster(QString inputPath, QString referencePath, QString outputPath, GDALResampleAlg eResample);
void BASECONSTVARIABLEAPI ResampleByReferenceRasterB(QString pszSrcFile, QString RefrasterBPath, QString pszOutFile, GDALResampleAlg eResample);
//--------------------- 保存文博 -------------------------------
int BASECONSTVARIABLEAPI saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path);

View File

@ -2,6 +2,7 @@
#include "ui_QResampleRefrenceRaster.h"
#include <QMessageBox>
#include <QFileDialog>
#include "ImageOperatorBase.h"
QResampleRefrenceRaster::QResampleRefrenceRaster(QWidget *parent)
: QDialog(parent),ui(new Ui::QResampleRefrenceRasterClass)
@ -83,12 +84,26 @@ void QResampleRefrenceRaster::ondialogBtnaccepted()
QString RefRasterPath = this->ui->lineEditRefRaster->text();
QString OutRasterPath = this->ui->lineEditOutRaster->text();
std::shared_ptr<double> gt(new double[6], delArrPtr);
gdalImage refimage(RefRasterPath);
gt.get()[0] = refimage.gt(0, 0);
gt.get()[1] = refimage.gt(0, 1);
gt.get()[2] = refimage.gt(0, 2);
gt.get()[3] = refimage.gt(1, 0);
gt.get()[4] = refimage.gt(1, 1);
gt.get()[5] = refimage.gt(1, 2);
ResampleGDAL(inRasterPath.toLocal8Bit().constData(), OutRasterPath.toLocal8Bit().constData(),
gt.get(), refimage.width, refimage.height,
GDALResampleAlg::GRA_Bilinear);
//alignRaster(inRasterPath, RefRasterPath, OutRasterPath,GDALResampleAlg::GRA_Bilinear);
QMessageBox::information(this, tr(u8"Ìáʾ"), tr(u8"completed!!!"));
}
void QResampleRefrenceRaster::ondialogBtnrejected()

View File

@ -120,7 +120,7 @@
</size>
</property>
<property name="text">
<string>输影像:</string>
<string>输影像:</string>
</property>
</widget>
</item>