From 23277b911e28b107bc3d50e9d1a7f5a85c52fe01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E9=99=88=E5=A2=9E=E8=BE=89?= <3045316072@qq.com> Date: Fri, 15 Nov 2024 14:21:50 +0800 Subject: [PATCH] =?UTF-8?q?=E6=9B=B4=E6=96=B0=E5=BD=B1=E5=83=8F=E6=8B=BC?= =?UTF-8?q?=E6=8E=A5=E4=BB=A3=E7=A0=81?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- BaseConstVariable.h | 5 +- BaseTool.cpp | 1 - BaseTool.h | 5 +- FileOperator.cpp | 23 ++ FileOperator.h | 1 + ImageOperatorBase.cpp | 834 +++++++++++++++++++++++++++++++++++++++++- ImageOperatorBase.h | 270 +++++--------- LogInfoCls.cpp | 4 + LogInfoCls.h | 4 + RasterToolBase.cpp | 268 ++++++++++++++ RasterToolBase.h | 79 ++++ 11 files changed, 1290 insertions(+), 204 deletions(-) create mode 100644 RasterToolBase.cpp create mode 100644 RasterToolBase.h diff --git a/BaseConstVariable.h b/BaseConstVariable.h index 7842d2e..3d21db4 100644 --- a/BaseConstVariable.h +++ b/BaseConstVariable.h @@ -50,7 +50,10 @@ const double e = sqrt(eSquare); const double earth_Re = 6378136.49; const double earth_Rp = (1 - 1 / f_inverse) * earth_Re; const double earth_We = 0.000072292115; - + +/*********************************************** openMap参数 ********************************************************************/ +static long Paral_num_thread = 6; + /*********************************************** 基础枚举 ********************************************************************/ diff --git a/BaseTool.cpp b/BaseTool.cpp index 602747a..c122680 100644 --- a/BaseTool.cpp +++ b/BaseTool.cpp @@ -115,7 +115,6 @@ std::complex Cubic_kernel_weight(double s) /// double Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12, Landpoint p22) { - return p11.ati * (1 - p0.lon) * (1 - p0.lat) + p12.ati * p0.lon * (1 - p0.lat) + p21.ati * (1 - p0.lon) * p0.lat + diff --git a/BaseTool.h b/BaseTool.h index 86153c5..b2b0220 100644 --- a/BaseTool.h +++ b/BaseTool.h @@ -26,8 +26,6 @@ #include #include #include - - #include "LogInfoCls.h" ///////////////////////////////////// 运行时间打印 @@ -56,8 +54,7 @@ std::complex Cubic_Convolution_interpolation(double u, double v, std::complex Cubic_kernel_weight(double s); -double Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12, - Landpoint p22); +double Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12,Landpoint p22); bool onSegment(Point3 Pi, Point3 Pj, Point3 Q); diff --git a/FileOperator.cpp b/FileOperator.cpp index 5ca06f9..3319db2 100644 --- a/FileOperator.cpp +++ b/FileOperator.cpp @@ -13,6 +13,29 @@ #include #include + + + +QString addMaskToFileName(const QString& filePath,QString _zzname) { + // 获取文件路径和文件名 + QFileInfo fileInfo(filePath); + + // 获取文件名和扩展名 + QString baseName = fileInfo.baseName(); // 不包含扩展名的文件名 + QString extension = fileInfo.suffix(); // 文件扩展名（例如 ".txt", ".jpg"） + + // 拼接新的文件名 + QString newFileName = baseName + _zzname; // 在文件名中增加 "_mask" + if (!extension.isEmpty()) { + newFileName += "." + extension; // 如果有扩展名，添加扩展名 + } + + // 返回新的文件路径 + QString newFilePath = fileInfo.path() + "/" + newFileName; + return newFilePath; +} + + std::vector getFilelist(const QString& folderpath, const QString& filenameExtension, int (*logfun)(QString logtext, int value)) { QString filenameExtensionStr = filenameExtension; diff --git a/FileOperator.h b/FileOperator.h index b6b936c..01f5131 100644 --- a/FileOperator.h +++ b/FileOperator.h @@ -45,5 +45,6 @@ size_t fsize(FILE* fp); QString getParantFromPath(const QString& path); void copyFile(const QString& sourcePath, const QString& destinationPath); +QString addMaskToFileName(const QString& filePath, QString _zzname); // QT FileOperator #endif \ No newline at end of file diff --git a/ImageOperatorBase.cpp b/ImageOperatorBase.cpp index 80abc78..5ee933c 100644 --- a/ImageOperatorBase.cpp +++ b/ImageOperatorBase.cpp @@ -569,24 +569,28 @@ Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, } delete[] temp; } - // else if (gdal_datatype == GDT_UInt64) { - // unsigned long* temp = new unsigned long[rows_count * cols_count]; - // demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, - //rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { for (int j = 0; j < - //cols_count; j++) { datamatrix(i, j) = temp[i * cols_count + j]; - // } - // } - // delete[] temp; - // } - // else if (gdal_datatype == GDT_Int64) { - // long* temp = new long[rows_count * cols_count]; - // demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, - //rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { for (int j = 0; j < - //cols_count; j++) { datamatrix(i, j) = temp[i * cols_count + j]; - // } - // } - // delete[] temp; - // } + else if (gdal_datatype == GDT_UInt64) { + unsigned long* temp = new unsigned long[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < + cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_Int64) { + long* temp = new long[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < + cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } else if(gdal_datatype == GDT_Float32) { float* temp = new float[rows_count * cols_count]; demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, @@ -618,6 +622,134 @@ Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, return datamatrix; } +Eigen::MatrixXi gdalImage::getDatai(int start_row, int start_col, int rows_count, int cols_count, int band_ids) +{ + omp_lock_t lock; + omp_init_lock(&lock); + omp_set_lock(&lock); + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); + GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen( + this->img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 + + GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); + GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids); + + rows_count = start_row + rows_count <= this->height ? rows_count : this->height - start_row; + cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col; + + Eigen::MatrixXi datamatrix(rows_count, cols_count); + + if (gdal_datatype == GDT_Byte) { + char* temp = new char[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); + for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_UInt16) { + unsigned short* temp = new unsigned short[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); + for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_Int16) { + short* temp = new short[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); + for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_UInt32) { + unsigned int* temp = new unsigned int[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); + for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_Int32) { + int* temp = new int[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); + for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_UInt64) { + unsigned long* temp = new unsigned long[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < + cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_Int64) { + long* temp = new long[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < + cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_Float32) { + float* temp = new float[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); + + for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_Float64) { + double* temp = new double[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, + rows_count, gdal_datatype, 0, 0); + + for (int i = 0; i < rows_count; i++) { + for (int j = 0; j < cols_count; j++) { + datamatrix(i, j) = temp[i * cols_count + j]; + } + } + delete[] temp; + } + else { + } + GDALClose((GDALDatasetH)rasterDataset); + omp_unset_lock(&lock); // 锟酵放伙拷斤拷 + omp_destroy_lock(&lock); // 劫伙拷斤拷 + // GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + return datamatrix; +} + Eigen::MatrixXd gdalImage::getGeoTranslation() { return this->gt; @@ -695,6 +827,61 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col omp_destroy_lock(&lock); // 劫伙拷斤拷 } +void gdalImage::saveImage(Eigen::MatrixXi data, int start_row, int start_col, int band_ids) +{ + omp_lock_t lock; + omp_init_lock(&lock); + omp_set_lock(&lock); + if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) { + QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") "; + qDebug() << tip; + throw std::exception(tip.toUtf8().constData()); + } + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); + GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + GDALDataset* poDstDS = nullptr; + if (exists_test(this->img_path)) { + poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update)); + } + else { + poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height, + this->band_num, GDT_Float32, NULL); // 斤拷锟斤拷 + poDstDS->SetProjection(this->projection.toUtf8().constData()); + + double gt_ptr[6]; + for (int i = 0; i < this->gt.rows(); i++) { + for (int j = 0; j < this->gt.cols(); j++) { + gt_ptr[i * 3 + j] = this->gt(i, j); + } + } + poDstDS->SetGeoTransform(gt_ptr); + //delete gt_ptr; + } + + long datarows = data.rows(); + long datacols = data.cols(); + + long* databuffer = new long[datarows * datacols]; // (float*)malloc(datarows * datacols * sizeof(float)); + + for (long i = 0; i < datarows; i++) { + for (long j = 0; j < datacols; j++) { + databuffer[i * datacols + j] = data(i, j); + } + } + // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols, + // datarows, GDT_Float32,band_ids, num,0,0,0); + poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows, + databuffer, datacols, datarows, GDT_Int64, 0, 0); + + GDALFlushCache(poDstDS); + GDALClose((GDALDatasetH)poDstDS); + // GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + delete[] databuffer; + omp_unset_lock(&lock); // 锟酵放伙拷斤拷 + omp_destroy_lock(&lock); // 劫伙拷斤拷 +} + void gdalImage::saveImage() { this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids); @@ -893,6 +1080,42 @@ Eigen::MatrixXd gdalImage::getHist(int bandids) return result; } +RasterExtend gdalImage::getExtend() +{ + RasterExtend extend{ 0,0,0,0 }; + double x1 = this->gt(0, 0); + double y1 = this->gt(1, 0); + + double x2 = this->gt(0, 0) + (this->width - 1) * gt(0, 1) + (0) * gt(0, 2); + double y2 = this->gt(1, 0) + (this->width - 1) * gt(1, 1) + (0) * gt(1, 2); + + double x3 = this->gt(0, 0) + (0) * gt(0, 1) + (this->height - 1) * gt(0, 2); + double y3 = this->gt(1, 0) + (0) * gt(1, 1) + (this->height - 1) * gt(1, 2); + + double x4 = this->gt(0, 0) + (this->width - 1) * gt(0, 1) + (this->height - 1) * gt(0, 2); + double y4 = this->gt(1, 0) + (this->width - 1) * gt(1, 1) + (this->height - 1) * gt(1, 2); + + + extend.min_x = x1 < x2 ? x1 : x2; + extend.max_x = x1 < x2 ? x2 : x1; + extend.min_y = y1 < y2 ? y1 : y2; + extend.max_y = y1 < y2 ? y2 : y1; + + + extend.min_x = extend.min_x < x3 ? extend.min_x : x3; + extend.max_x = extend.max_x > x3 ? extend.max_x : x3; + extend.min_y = extend.min_y < y3 ? extend.min_y : y3; + extend.max_y = extend.max_y > y3 ? extend.max_y : y3; + + + extend.min_x = extend.min_x < x4 ? extend.min_x : x4; + extend.max_x = extend.max_x > x4 ? extend.max_x : x4; + extend.min_y = extend.min_y < y4 ? extend.min_y : y4; + extend.max_y = extend.max_y > y4 ? extend.max_y : y4; + + return extend; +} + gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite) { @@ -930,6 +1153,52 @@ gdalImage CreategdalImage(const QString& img_path, int height, int width, int ba return result_img; } +gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long epsgCode, GDALDataType eType, bool need_gt, bool overwrite, bool isENVI) +{ + if (exists_test(img_path.toUtf8().constData())) { + if (overwrite) { + gdalImage result_img(img_path); + return result_img; + } + else { + throw "file has exist!!!"; + exit(1); + } + } + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7 + GDALDriver* poDriver = isENVI? GetGDALDriverManager()->GetDriverByName("ENVI"): GetGDALDriverManager()->GetDriverByName("GTiff"); + GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, eType, NULL); // 锟斤拷锟斤拷锟斤拷 + if (need_gt) { + OGRSpatialReference oSRS; + + if (oSRS.importFromEPSG(epsgCode) != OGRERR_NONE) { + std::cerr << "Failed to import EPSG code " << epsgCode << std::endl; + throw "Failed to import EPSG code "; + exit(1); + } + char* pszWKT = NULL; + oSRS.exportToWkt(&pszWKT); + std::cout << "WKT of EPSG:"<< epsgCode <<" :\n" << pszWKT << std::endl; + poDstDS->SetProjection(pszWKT); + double gt_ptr[6] = { 0 }; + for (int i = 0; i < gt.rows(); i++) { + for (int j = 0; j < gt.cols(); j++) { + gt_ptr[i * 3 + j] = gt(i, j); + } + } + poDstDS->SetGeoTransform(gt_ptr); + } + for (int i = 1; i <= band_num; i++) { + poDstDS->GetRasterBand(i)->SetNoDataValue(-9999); + } + GDALFlushCache((GDALDatasetH)poDstDS); + GDALClose((GDALDatasetH)poDstDS); + ////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + gdalImage result_img(img_path); + return result_img; +} + gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite) @@ -1275,6 +1544,275 @@ int saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path) return -1; } +ErrorCode MergeRasterProcess(QVector filepaths, QString outfileptah, QString mainString, MERGEMODE mergecode, bool isENVI) +{ + // 参数检查 + if (!isExists(mainString)) { + qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND) )<< "\t" << mainString; + return ErrorCode::FILENOFOUND; + } + else {} + gdalImage mainimg(mainString); + QVector imgdslist(filepaths.count()); + for (long i = 0; i < filepaths.count(); i++) { + if (!isExists(filepaths[i])) { + qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND)) << "\t" << filepaths[i]; + return ErrorCode::FILENOFOUND; + } + else { + imgdslist[i] = gdalImage(filepaths[i]); + if (imgdslist[i].band_num != mainimg.band_num) { + qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTERBAND_NOTEQUAL)) << "\t" << imgdslist[i].band_num <<" != "<< mainimg.band_num; + return ErrorCode::RASTERBAND_NOTEQUAL; + } + } + } + + // 检查坐标系是否统一 + long EPSGCode = GetEPSGFromRasterFile(mainString); + long tempCode = 0; + for (long i = 0; i < filepaths.count(); i++) { + tempCode = GetEPSGFromRasterFile(filepaths[i]); + if (EPSGCode != tempCode) { + qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSAME)) << "\t" << EPSGCode <<"!="<< tempCode; + return ErrorCode::EPSGCODE_NOTSAME; + } + } + + // 检查影像类型是否统一 + GDALDataType mainType = mainimg.getDataType(); + for (long i = 0; i < imgdslist.count(); i++) { + if (mainType != imgdslist[i].getDataType()) { + qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTER_DATETYPE_NOTSAME)) << "\t" << mainType << "!=" << imgdslist[i].getDataType(); + return ErrorCode::RASTER_DATETYPE_NOTSAME; + } + } + + + + Eigen::MatrixXd maingt = mainimg.getGeoTranslation(); + Eigen::MatrixXd rgt = Eigen::MatrixXd::Zero(2,3); + RasterExtend mainExtend = mainimg.getExtend(); + rgt(0, 1) = (mainExtend.max_x - mainExtend.min_x) / (mainimg.width - 1); //dx + rgt(1, 2) =-1*( (mainExtend.max_y - mainExtend.min_y) / (mainimg.height - 1));//dy + QVector extendlist(imgdslist.count()); + for (long i = 0; i < imgdslist.count(); i++) { + extendlist[i] = imgdslist[i].getExtend(); + mainExtend.min_x = mainExtend.min_x < extendlist[i].min_x ? mainExtend.min_x : extendlist[i].min_x; + mainExtend.max_x = mainExtend.max_x > extendlist[i].max_x ? mainExtend.max_x : extendlist[i].max_x; + mainExtend.min_y = mainExtend.min_y < extendlist[i].min_y ? mainExtend.min_y : extendlist[i].min_y; + mainExtend.max_y = mainExtend.max_y > extendlist[i].max_y ? mainExtend.max_y : extendlist[i].max_y; + } + + rgt(0, 0) = mainExtend.min_x; + rgt(1, 0) = mainExtend.max_y; + + // 计算数量 + + long width = std::ceil((mainExtend.max_x - mainExtend.min_x) / rgt(0, 1) + 1); + long height = std::ceil(std::abs((mainExtend.min_y - mainExtend.max_y) / rgt(1, 2)) + 1); + OGRSpatialReference oSRS; + if (oSRS.importFromEPSG(EPSGCode) != OGRERR_NONE) { + qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSUPPORT)) << "\t" << EPSGCode; + return ErrorCode::EPSGCODE_NOTSUPPORT; + } + + gdalImage resultImage = CreategdalImage(outfileptah, height, width, mainimg.band_num, rgt, EPSGCode, mainType, true, true, isENVI); + + QString resultMaskString = addMaskToFileName(outfileptah, QString("_MASK")); + gdalImage maskImage = CreategdalImage(resultMaskString, height, width,1 , rgt, EPSGCode, GDT_UInt16, true, true, isENVI); + + // 初始化 + + // 分块计算 + long resultline = Memory1MB * 100 / 8 / resultImage.width; + resultline = resultline < 100 ? resultline : 100; // 最多100行 + + long bandnum = resultImage.band_num + 1; + long starti = 0; + long rasterCount = imgdslist.count(); + + for (starti = 0; starti < resultImage.height; starti = starti + resultline) { + long blocklines = resultline; + blocklines = starti + blocklines < resultImage.height ? blocklines : resultImage.height - starti; + for (long b = 1; b < bandnum; b++) { + Eigen::MatrixXd data = resultImage.getData(starti, 0, blocklines, resultImage.width, b); + Eigen::MatrixXd maskdata = maskImage.getData(starti, 0, blocklines, resultImage.width, b); + data = data.array() * 0; + maskdata = maskdata.array() * 0; + resultImage.saveImage(data, starti, 0, b); + maskImage.saveImage(maskdata, starti, 0, b); + } + } + + + + + switch (mergecode) + { + case MERGE_GEOCODING: + return MergeRasterInGeoCoding(imgdslist, resultImage, maskImage); + default: + break; + } + + + return ErrorCode::SUCCESS; +} + +ErrorCode MergeRasterInGeoCoding(QVector imgdslist, gdalImage resultimg, gdalImage maskimg) +{ + omp_set_num_threads(Paral_num_thread); + // 逐点合并计算 + QVector extendlist(imgdslist.count()); + for (long i = 0; i < imgdslist.count(); i++) { + extendlist[i] = imgdslist[i].getExtend(); + imgdslist[i].InitInv_gt(); + } + + // 分块计算 + long resultline = Memory1MB * 100 / 8 / resultimg.width; + resultline = resultline < 100 ? resultline : 100; // 最多100行 + + long bandnum = resultimg.band_num+1; + long starti = 0; + long rasterCount = imgdslist.count(); + + long processNumber = 0; + + omp_lock_t lock; + omp_init_lock(&lock); + + + +#pragma omp parallel for + for (starti = 0; starti < resultimg.height; starti = starti + resultline) { + long blocklines = resultline; + blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti; + + long rid = starti; + long cid = 0; + + Landpoint pp = {0,0,0}; + Landpoint lp = { 0,0,0 }; + + for (long ir = 0; ir < rasterCount; ir++) {// 影像 + long minRid = imgdslist[ir].height; + long maxRid = 0; + + + Eigen::MatrixXd ridlist = Eigen::MatrixXd::Zero(blocklines, resultimg.width); + Eigen::MatrixXd cidlist = Eigen::MatrixXd::Zero(blocklines, resultimg.width); + + for (long i = 0; i < blocklines; i++) {// 行号 + rid = starti + i; + for (long j = 0; j < resultimg.width; j++) {// 列号 + cid = j; + resultimg.getLandPoint(i,j,0,pp); + lp = resultimg.getRow_Col(pp.lon, pp.lat); // 获取点坐标 + ridlist(i, j) = lp.lat; + cidlist(i, j) = lp.lon; + } + } + + if (ridlist.maxCoeff() < 0 || ridlist.minCoeff() >= imgdslist[ir].height) { + continue; + } + + if (cidlist.maxCoeff() < 0 || cidlist.minCoeff() >= imgdslist[ir].width) { + continue; + } + + minRid = std::floor(ridlist.minCoeff()); + maxRid = std::ceil(ridlist.maxCoeff()); + minRid = minRid < 0 ? 0 : minRid; + maxRid = maxRid < imgdslist[ir].height ? maxRid : imgdslist[ir].height - 1; + + long rowlen = maxRid - minRid + 1; + if(rowlen <= 0) { + continue; + } + // 获取分配代码 + Landpoint p0{ 0,0,0 }, p11{ 0,0,0 }, p21{ 0,0,0 }, p12{ 0,0,0 }, p22{ 0,0,0 }, p{ 0,0,0 }; + + + long rowcount = 0; + long colcount = 0; + double ridtemp = 0, cidtemp = 0; + + long lastr = 0, nextr = 0; + long lastc = 0, nextc = 0; + + double r0=0, c0 = 0; + + for (long b = 1; b < bandnum; b++) { + Eigen::MatrixXd resultdata = resultimg.getData(starti, 0, blocklines, resultimg.width, b); + Eigen::MatrixXi resultmask = maskimg.getDatai(starti, 0, blocklines, resultimg.width, b); + Eigen::MatrixXd data = imgdslist[ir].getData(minRid, 0, rowlen, imgdslist[ir].width, b); + rowcount = data.rows(); + colcount = data.cols(); + for (long i = 0; i < rowcount; i++) { + for (long j = 0; j < colcount; j++) { + ridtemp = ridlist(i, j); + cidtemp = cidlist(i, j); + + if (ridtemp < 0 || ridtemp >= imgdslist[ir].height||cidtemp<0||cidtemp>=imgdslist[ir].width) { + continue; + } + else {} + + r0 = ridtemp - std::floor(ridtemp); + c0 = cidtemp - std::floor(cidtemp); + + p0 = Landpoint{ c0,r0,0 }; + p11 = Landpoint{ 0,0,data(lastr,lastc) }; + p21 = Landpoint{ 0,1,data(nextr,lastc) }; + p12 = Landpoint{ 1,0,data(lastr,nextc) }; + p22 = Landpoint{ 1,1,data(nextr,nextc) }; + + resultdata(i,j) = resultdata(i, j)+ Bilinear_interpolation(p0, p11, p21, p12, p22); + resultmask(i, j) = resultmask(i, j) + 1; + } + } + resultimg.saveImage(resultdata, starti, 0, b); + maskimg.saveImage(resultmask, starti, 0, b); + } + } + + omp_set_lock(&lock); + processNumber = processNumber + blocklines; + std::cout << "\rprocess bar:\t" << processNumber * 100.0 / resultimg.height << " % " << "\t\t\t"; + omp_unset_lock(&lock); // 锟酵放伙拷斤拷 + } + omp_destroy_lock(&lock); // 劫伙拷斤拷 + std::cout << std::endl; + + // 全体处理 + for (starti = 0; starti < resultimg.height; starti = starti + resultline) { + long blocklines = resultline; + blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti; + for (long b = 1; b < bandnum; b++) { + Eigen::MatrixXd data = resultimg.getData(starti, 0, blocklines, resultimg.width, b); + Eigen::MatrixXd maskdata = maskimg.getData(starti, 0, blocklines, maskimg.width, b); + + for (long i = 0; i < data.rows(); i++) { + for (long j = 0; j < data.cols(); j++) { + if (maskdata(i, j) == 0) { + continue; + } + data(i, j) = data(i, j) / maskdata(i, j); + } + } + + resultimg.saveImage(data, starti, 0, b); + maskimg.saveImage(maskdata, starti, 0, b); + } + } + + + return ErrorCode::SUCCESS; +} + gdalImageComplex::gdalImageComplex(const QString& raster_path) { omp_lock_t lock; @@ -1505,7 +2043,265 @@ void gdalImageComplex::savePreViewImage() - +long getProjectEPSGCodeByLon_Lat(double lon, double lat, ProjectStripDelta stripState) +{ + long EPSGCode = 0; + switch (stripState) { + case ProjectStripDelta::Strip_3: { + break; + }; + case ProjectStripDelta::Strip_6: { + break; + } + default: { + EPSGCode = -1; + break; + } + } + qDebug() << QString(" EPSG code : %1").arg(EPSGCode); + return EPSGCode; +} +long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat) +{ + // EPSG 4534 ~ 4554 3 度带 + // 首先判断是否是在 中国带宽范围 + // 中心经度范围 ：75E ~ 135E 实际范围 73.5E ~ 136.5E, + // 纬度范围 3N ~ 54N，放宽到 0N~ 60N + if (73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内 + long code = trunc((lon - 73.5) / 3) + 4534; + return code; + } + else { // 非中国境内 使用 高斯克吕格 + bool isSouth = lat < 0; // 简单判断南北半球，这里仅为示例，实际应用可能需要更细致的逻辑 + long prefix = isSouth ? 327000 : 326000; + // std::string prefix = isSouth ? "327" : "326"; + lon = fmod(lon + 360.0, 360.0); + long zone = std::floor((lon + 180.0) / 3.0); + prefix = prefix + zone; + return prefix; + } + return 0; +} +long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat) +{ + // EPSG 4502 ~ 4512 6度带 + // 首先判断是否是在 中国带宽范围 + // 中心经度范围 ：75E ~ 135E 实际范围 72.0E ~ 138E, + // 纬度范围 3N ~ 54N，放宽到 0N~ 60N + if (73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内 + long code = trunc((lon - 72.0) / 6) + 4502; + return code; + } + else { // 非中国境内 使用 UTM// 确定带号，6度带从1开始到60，每6度一个带 + int zone = static_cast((lon + 180.0) / 6.0) + 1; + bool isSouth = lon < 0; // 判断是否在南半球 + long epsgCodeBase = isSouth ? 32700 : 32600; // 计算EPSG代码 + long prefix = static_cast(epsgCodeBase + zone); + return prefix; + } + return 0; +} + +QString GetProjectionNameFromEPSG(long epsgCode) +{ + qDebug() << "============= GetProjectionNameFromEPSG ======================"; + OGRSpatialReference oSRS; + + // 设置EPSG代码 + if (oSRS.importFromEPSG(epsgCode) != OGRERR_NONE) { + qDebug() << "epsgcode not recognition"; + return ""; + } + + // 获取并输出坐标系的描述（名称） + const char* pszName = oSRS.GetAttrValue("GEOGCS"); + if (pszName) { + qDebug() << "Coordinate system name for EPSG " + QString::number(epsgCode) + << " is: " + QString::fromStdString(pszName); + return QString::fromStdString(pszName); + } + else { + qDebug() << "Unable to retrieve the name for EPSG " + QString::number(epsgCode); + return ""; + } + + // char* wkt = NULL; + // // 转换为WKT格式 + // oSRS.exportToWkt(&wkt); + // + // qDebug() << wkt; + // + // // 从WKT中解析投影名称，这里简化处理，实际可能需要更复杂的逻辑来准确提取名称 + // std::string wktStr(wkt); + // long start = wktStr.find("PROJCS[\"") + 8; // 找到"PROJCS["后的第一个双引号位置 + // // 从start位置开始找下一个双引号，这之间的内容即为投影名称 + // int end = wktStr.find('\"', start); + // QString projName = QString::fromStdString(wktStr.substr(start, end - start)); + // + // // 释放WKT字符串内存 + // CPLFree(wkt); + + // return projName; +} +long GetEPSGFromRasterFile(QString filepath) +{ + qDebug() << "============= GetEPSGFromRasterFile ======================"; + // QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器 + // QByteArray byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray + //,这个应该会自动释放 const char* charArray = byteArray.constData(); // + // 获取QByteArray的const char*指针 + + { + if (QFile(filepath).exists()) { + qDebug() << "info: the image found.\n"; + } + else { + return -1; + } + + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注册GDAL驱动 + // std::cout<GetProjectionRef(); + + qDebug() << QString::fromUtf8(pszProjection); + + // 创建SpatialReference对象 + OGRSpatialReference oSRS; + if (oSRS.importFromWkt((char**)&pszProjection) != OGRERR_NONE) { + qDebug() << ("Error: Unable to import projection information.\n"); + GDALClose(poDataset); + return -1; + } + + long epsgCode = atoi(oSRS.GetAuthorityCode(nullptr)); // 获取EPSG代码 + + if (epsgCode != 0) { + GDALClose(poDataset); + qDebug() << QString("file %1 :epsg Code %2").arg(filepath).arg(epsgCode); + return epsgCode; + } + else { + qDebug() << "EPSG code could not be determined from the spatial reference."; + GDALClose(poDataset); + return -1; + } + } +} +std::shared_ptr GetCenterPointInRaster(QString filepath) +{ + qDebug() << "============= GetCenterPointInRaster ======================"; + // QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器 + // QByteArray byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray + //,这个应该会自动释放 const char* charArray = byteArray.constData(); // + // 获取QByteArray的const char*指针 + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); + // std::cout<GetGeoTransform(adfGeoTransform) != CE_None) { + qDebug() << "Failed to get GeoTransform"; + return nullptr; + } + + double dfWidth = poDataset->GetRasterXSize(); + double dfHeight = poDataset->GetRasterYSize(); + + // 计算中心点坐标（像素坐标） + double dfCenterX = adfGeoTransform[0] + dfWidth * adfGeoTransform[1] / 2.0 + + dfHeight * adfGeoTransform[2] / 2.0; + double dfCenterY = adfGeoTransform[3] + dfWidth * adfGeoTransform[4] / 2.0 + + dfHeight * adfGeoTransform[5] / 2.0; + + OGRSpatialReference oSRS; + oSRS.importFromWkt(poDataset->GetProjectionRef()); + + if (oSRS.IsGeographic()) { + qDebug() << "Center coords (already in geographic): (" + QString::number(dfCenterX) + + ", " + QString::number(dfCenterY) + ")"; + flag = true; + x = dfCenterX; + y = dfCenterY; + } + else { + // 如果不是地理坐标系，转换到WGS84 + OGRSpatialReference oSRS_WGS84; + oSRS_WGS84.SetWellKnownGeogCS("WGS84"); + + OGRCoordinateTransformation* poCT = + OGRCreateCoordinateTransformation(&oSRS, &oSRS_WGS84); + if (poCT == nullptr) { + qDebug() << "Failed to create coordinate transformation"; + return nullptr; + } + + // double dfLon, dfLat; + if (poCT->Transform(1, &dfCenterX, &dfCenterY)) { + qDebug() << "Center coords (transformed to WGS84): (" + + QString::number(dfCenterX) + ", " + QString::number(dfCenterY) + << ")"; + flag = true; + x = dfCenterX; + y = dfCenterY; + } + else { + qDebug() << "Transformation failed."; + } + OGRCoordinateTransformation::DestroyCT(poCT); + } + } + if (nullptr == poDataset || NULL == poDataset) {} + else { + GDALClose(poDataset); + } + + if (flag) { + std::shared_ptr RasterCenterPoint = std::make_shared(); + RasterCenterPoint->x = x; + RasterCenterPoint->y = y; + RasterCenterPoint->z = 0; + return RasterCenterPoint; + } + else { + return nullptr; + } +} +CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long epsg_code) +{ + OGRSpatialReference oSRS; + if (oSRS.importFromEPSG(epsg_code) == OGRERR_NONE) { + if (oSRS.IsGeographic()) { + return CoordinateSystemType::GeoCoordinateSystem; + } + else if (oSRS.IsProjected()) { + return CoordinateSystemType::ProjectCoordinateSystem; + } + } + else { + return CoordinateSystemType::UNKNOW; + } +} diff --git a/ImageOperatorBase.h b/ImageOperatorBase.h index 7d87945..24efe21 100644 --- a/ImageOperatorBase.h +++ b/ImageOperatorBase.h @@ -24,6 +24,39 @@ #include #include #include // for CPLMalloc() +#include "LogInfoCls.h" + + +enum ProjectStripDelta { + Strip_6, // 6度带 + Strip_3 +}; + +enum CoordinateSystemType { // 坐标系类型 + GeoCoordinateSystem, + ProjectCoordinateSystem, + UNKNOW +}; + +struct PointRaster { // 影像坐标点 + double x; + double y; + double z; +}; + + +struct PointXYZ { + double x, y, z; +}; + +struct PointGeo { + double lon, lat, ati; +}; + +struct PointImage { + double pixel_x, pixel_y; +}; + struct ImageGEOINFO { @@ -38,10 +71,41 @@ enum GDALREADARRCOPYMETHOD { VARIABLEMETHOD // 变量赋值 }; + -// 判断是否需要输出为DLL -#define DLLOUT + +/// 根据经纬度获取 +/// EPSG代码，根据经纬度返回对应投影坐标系统，其中如果在中华人民共和国境内，默认使用 +/// CGCS2000坐标系统(EPSG 4502 ~ 4512 6度带,EPSG 4534 ~ 4554 3度带)，其余地方使用WGS坐标系统， +/// 投影方法 高斯克吕格（国内）， 高斯克吕格 +/// \param long 经度 +/// \param lat 纬度 +/// \return 对应投影坐标系统的 EPSG编码,-1 表示计算错误 +long getProjectEPSGCodeByLon_Lat(double long, double lat,ProjectStripDelta stripState = ProjectStripDelta::Strip_3); + +long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat); + +long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat); + + +QString GetProjectionNameFromEPSG(long epsgCode); + + +long GetEPSGFromRasterFile(QString filepath); + +std::shared_ptr GetCenterPointInRaster(QString filepath); + +CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long EPSGCODE); + + + + + + + + + // 文件打开 // 当指令销毁时，调用GDALClose 销毁类型 std::shared_ptr OpenDataset(const QString& in_path, GDALAccess rwmode= GA_ReadOnly); void CloseDataset(GDALDataset* ptr); @@ -70,11 +134,11 @@ ImageGEOINFO getImageINFO(QString in_path); GDALDataType getGDALDataType(QString fileptah); -struct DemBox { - double min_lat; //纬度 - double min_lon;//经度 - double max_lat;//纬度 - double max_lon;//经度 +struct RasterExtend { + double min_x; //纬度 + double min_y;//经度 + double max_x;//纬度 + double max_y;//经度 }; @@ -97,9 +161,12 @@ public: // 方法 virtual void setTranslationMatrix(Eigen::MatrixXd gt); virtual void setData(Eigen::MatrixXd,int data_band_ids=1); virtual Eigen::MatrixXd getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids); + virtual Eigen::MatrixXi getDatai(int start_row, int start_col, int rows_count, int cols_count, int band_ids); + virtual Eigen::MatrixXd getGeoTranslation(); virtual GDALDataType getDataType(); virtual void saveImage(Eigen::MatrixXd, int start_row, int start_col, int band_ids); + virtual void saveImage(Eigen::MatrixXi, int start_row, int start_col, int band_ids); virtual void saveImage(); virtual void setNoDataValue(double nodatavalue, int band_ids); virtual int InitInv_gt(); @@ -116,6 +183,7 @@ public: // 方法 virtual Eigen::MatrixXd getHist(int bandids); + virtual RasterExtend getExtend(); public: QString img_path; // 图像文件 @@ -154,13 +222,13 @@ public: // 创建影像 gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false); +gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long espgcode, GDALDataType eType=GDT_Float32, bool need_gt = true, bool overwrite = false,bool isENVI=false); gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false); gdalImageComplex CreateEchoComplex(const QString& img_path, int height, int width, int band_num); - int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample); int ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResample = GRIORA_Bilinear); @@ -173,180 +241,24 @@ int saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path); //---------------------------------------------------- +//--------------------- 图像合并流程 ------------------------------ +enum MERGEMODE +{ + MERGE_GEOCODING, +}; + + +ErrorCode MergeRasterProcess(QVector filepath, QString outfileptah, QString mainString, MERGEMODE mergecode = MERGEMODE::MERGE_GEOCODING, bool isENVI = false); + + +ErrorCode MergeRasterInGeoCoding(QVector inimgs, gdalImage resultimg,gdalImage maskimg); + + + template std::shared_ptr readDataArr(gdalImage &img,int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method); -#ifndef DLLOUT - - - - - -#else -//#define DllExport __declspec( dllexport ) -//double __declspec(dllexport) ProcessMGCMathX_MGC(int Xbetaidx, int Xbwidx, double XTao, double satH, char* sigma_path, char* output_path, -// double p1_x, double p1_y, double p2_x, double p2_y, double p3_x, double p3_y, double p4_x, double p4_y) - #endif - -#endif - -template -std::shared_ptr readDataArr(gdalImage& imgds, int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method) -{ - std::shared_ptr result = nullptr; - - - omp_lock_t lock; - omp_init_lock(&lock); - omp_set_lock(&lock); - GDALAllRegister(); - CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); - GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(imgds.img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 - - GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); - GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids); - - rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row; - cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col; - - Eigen::MatrixXd datamatrix(rows_count, cols_count); - - if (gdal_datatype == GDT_Byte) { - char* temp = new char[rows_count * cols_count]; - demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); - result=std::shared_ptr(new T[rows_count * cols_count], delArrPtr); - if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { - std::memcpy(result.get(), temp, rows_count * cols_count); - } - else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { - long count = rows_count * cols_count; - for (long i = 0; i < count; i++) { - result.get()[i] = T(temp[i]); - } - } - - - delete[] temp; - } - else if (gdal_datatype == GDT_UInt16) { - unsigned short* temp = new unsigned short[rows_count * cols_count]; - demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); - result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); - if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { - std::memcpy(result.get(), temp, rows_count * cols_count); - } - else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { - long count = rows_count * cols_count; - for (long i = 0; i < count; i++) { - result.get()[i] = T(temp[i]); - } - } - delete[] temp; - } - else if (gdal_datatype == GDT_Int16) { - short* temp = new short[rows_count * cols_count]; - demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); - result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); - if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { - std::memcpy(result.get(), temp, rows_count * cols_count); - } - else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { - long count = rows_count * cols_count; - for (long i = 0; i < count; i++) { - result.get()[i] = T(temp[i]); - } - } - delete[] temp; - } - else if (gdal_datatype == GDT_UInt32) { - unsigned int* temp = new unsigned int[rows_count * cols_count]; - demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); - result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); - if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { - std::memcpy(result.get(), temp, rows_count * cols_count); - } - else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { - long count = rows_count * cols_count; - for (long i = 0; i < count; i++) { - result.get()[i] = T(temp[i]); - } - } - delete[] temp; - } - else if (gdal_datatype == GDT_Int32) { - int* temp = new int[rows_count * cols_count]; - demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); - result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); - if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { - std::memcpy(result.get(), temp, rows_count * cols_count); - } - else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { - long count = rows_count * cols_count; - for (long i = 0; i < count; i++) { - result.get()[i] = T(temp[i]); - } - } - delete[] temp; - } - // else if (gdal_datatype == GDT_UInt64) { - // unsigned long* temp = new unsigned long[rows_count * cols_count]; - // demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, - //rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { for (int j = 0; j < - //cols_count; j++) { datamatrix(i, j) = temp[i * cols_count + j]; - // } - // } - // delete[] temp; - // } - // else if (gdal_datatype == GDT_Int64) { - // long* temp = new long[rows_count * cols_count]; - // demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, - //rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { for (int j = 0; j < - //cols_count; j++) { datamatrix(i, j) = temp[i * cols_count + j]; - // } - // } - // delete[] temp; - // } - else if (gdal_datatype == GDT_Float32) { - float* temp = new float[rows_count * cols_count]; - demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); - - result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); - if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { - std::memcpy(result.get(), temp, rows_count * cols_count); - } - else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { - long count = rows_count * cols_count; - for (long i = 0; i < count; i++) { - result.get()[i] = T(temp[i]); - } - } - delete[] temp; - } - else if (gdal_datatype == GDT_Float64) { - double* temp = new double[rows_count * cols_count]; - demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); - - result = std::shared_ptr(new T[rows_count * cols_count], delArrPtr); - if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) { - std::memcpy(result.get(), temp, rows_count * cols_count); - } - else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) { - long count = rows_count * cols_count; - for (long i = 0; i < count; i++) { - result.get()[i] = T(temp[i]); - } - } - delete[] temp; - } - else { - } - GDALClose((GDALDatasetH)rasterDataset); - omp_unset_lock(&lock); // 锟酵放伙拷斤拷 - omp_destroy_lock(&lock); // 劫伙拷斤拷 - // GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH - return result; -} - + \ No newline at end of file diff --git a/LogInfoCls.cpp b/LogInfoCls.cpp index 47cba8f..fcdb257 100644 --- a/LogInfoCls.cpp +++ b/LogInfoCls.cpp @@ -18,6 +18,10 @@ std::string errorCode2errInfo(ErrorCode e) _CASE_STR(FILE_NOT_EXIST); _CASE_STR(FIND_ID_ERROR); _CASE_STR(INSERT_ID_ERROR); + _CASE_STR(EPSGCODE_NOTSAME); + _CASE_STR(EPSGCODE_NOTSUPPORT); + _CASE_STR(RASTERBAND_NOTEQUAL); + _CASE_STR(RASTER_DATETYPE_NOTSAME); //GSL 1xx _CASE_STR(Error_GSL_FAILURE ); _CASE_STR(Error_GSL_CONTINUE ); diff --git a/LogInfoCls.h b/LogInfoCls.h index c439b9c..b35520d 100644 --- a/LogInfoCls.h +++ b/LogInfoCls.h @@ -27,6 +27,10 @@ enum ErrorCode { FILE_NOT_EXIST = 9, FIND_ID_ERROR = 10, INSERT_ID_ERROR = 11, + EPSGCODE_NOTSAME = 12, + EPSGCODE_NOTSUPPORT = 13, + RASTERBAND_NOTEQUAL = 14, + RASTER_DATETYPE_NOTSAME = 15, //GSL 1xx Error_GSL_FAILURE = -101, Error_GSL_CONTINUE = -102, /* iteration has not converged */ diff --git a/RasterToolBase.cpp b/RasterToolBase.cpp new file mode 100644 index 0000000..272c185 --- /dev/null +++ b/RasterToolBase.cpp @@ -0,0 +1,268 @@ +/** + * @file RasterProjectBase.cpp + * @brief None + * @author 陈增辉 (3045316072@qq.com) + * @version 2.5.0 + * @date 24-6-4 + * @copyright Copyright (c) Since 2024 中科卫星应用研究院 All rights reserved. + */ + +#include +#include "RasterToolBase.h" +#include "gdal_priv.h" +#include "ogr_spatialref.h" +#include "cpl_conv.h" // for CPLMalloc() +#include +#include +#include + +namespace RasterToolBase { + long getProjectEPSGCodeByLon_Lat(double lon, double lat, ProjectStripDelta stripState) + { + long EPSGCode = 0; + switch(stripState) { + case ProjectStripDelta::Strip_3: { + break; + }; + case ProjectStripDelta::Strip_6: { + break; + } + default: { + EPSGCode = -1; + break; + } + } + qDebug() << QString(" EPSG code : %1").arg(EPSGCode); + return EPSGCode; + } + long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat) + { + // EPSG 4534 ~ 4554 3 度带 + // 首先判断是否是在 中国带宽范围 + // 中心经度范围 ：75E ~ 135E 实际范围 73.5E ~ 136.5E, + // 纬度范围 3N ~ 54N，放宽到 0N~ 60N + if(73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内 + long code = trunc((lon - 73.5) / 3) + 4534; + return code; + } else { // 非中国境内 使用 高斯克吕格 + bool isSouth = lat < 0; // 简单判断南北半球，这里仅为示例，实际应用可能需要更细致的逻辑 + long prefix = isSouth ? 327000 : 326000; + // std::string prefix = isSouth ? "327" : "326"; + lon = fmod(lon + 360.0, 360.0); + long zone = std::floor((lon + 180.0) / 3.0); + prefix = prefix + zone; + return prefix; + } + return 0; + } + long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat) + { + // EPSG 4502 ~ 4512 6度带 + // 首先判断是否是在 中国带宽范围 + // 中心经度范围 ：75E ~ 135E 实际范围 72.0E ~ 138E, + // 纬度范围 3N ~ 54N，放宽到 0N~ 60N + if(73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内 + long code = trunc((lon - 72.0) / 6) + 4502; + return code; + } else { // 非中国境内 使用 UTM// 确定带号，6度带从1开始到60，每6度一个带 + int zone = static_cast((lon + 180.0) / 6.0) + 1; + bool isSouth = lon < 0; // 判断是否在南半球 + long epsgCodeBase = isSouth ? 32700 : 32600; // 计算EPSG代码 + long prefix = static_cast(epsgCodeBase + zone); + return prefix; + } + return 0; + } + + QString GetProjectionNameFromEPSG(long epsgCode) + { + qDebug() << "============= GetProjectionNameFromEPSG ======================"; + OGRSpatialReference oSRS; + + // 设置EPSG代码 + if(oSRS.importFromEPSG(epsgCode) != OGRERR_NONE) { + qDebug() << "epsgcode not recognition"; + return ""; + } + + // 获取并输出坐标系的描述（名称） + const char* pszName = oSRS.GetAttrValue("GEOGCS"); + if(pszName) { + qDebug() << "Coordinate system name for EPSG " + QString::number(epsgCode) + << " is: " + QString::fromStdString(pszName); + return QString::fromStdString(pszName); + } else { + qDebug() << "Unable to retrieve the name for EPSG " + QString::number(epsgCode); + return ""; + } + + // char* wkt = NULL; + // // 转换为WKT格式 + // oSRS.exportToWkt(&wkt); + // + // qDebug() << wkt; + // + // // 从WKT中解析投影名称，这里简化处理，实际可能需要更复杂的逻辑来准确提取名称 + // std::string wktStr(wkt); + // long start = wktStr.find("PROJCS[\"") + 8; // 找到"PROJCS["后的第一个双引号位置 + // // 从start位置开始找下一个双引号，这之间的内容即为投影名称 + // int end = wktStr.find('\"', start); + // QString projName = QString::fromStdString(wktStr.substr(start, end - start)); + // + // // 释放WKT字符串内存 + // CPLFree(wkt); + + // return projName; + } + long GetEPSGFromRasterFile(QString filepath) + { + qDebug() << "============= GetEPSGFromRasterFile ======================"; + // QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器 + // QByteArray byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray + //,这个应该会自动释放 const char* charArray = byteArray.constData(); // + // 获取QByteArray的const char*指针 + + { + if(QFile(filepath).exists()){ + qDebug() << "info: the image found.\n"; + }else{ + return -1; + } + + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注册GDAL驱动 + // std::cout<GetProjectionRef(); + + qDebug() << QString::fromUtf8(pszProjection); + + // 创建SpatialReference对象 + OGRSpatialReference oSRS; + if(oSRS.importFromWkt((char**)&pszProjection) != OGRERR_NONE) { + qDebug() << ("Error: Unable to import projection information.\n"); + GDALClose(poDataset); + return -1; + } + + long epsgCode = atoi(oSRS.GetAuthorityCode(nullptr)); // 获取EPSG代码 + + if(epsgCode != 0) { + GDALClose(poDataset); + qDebug() << QString("file %1 :epsg Code %2").arg(filepath).arg(epsgCode); + return epsgCode; + } else { + qDebug() << "EPSG code could not be determined from the spatial reference."; + GDALClose(poDataset); + return -1; + } + } + } + std::shared_ptr GetCenterPointInRaster(QString filepath) + { + qDebug() << "============= GetCenterPointInRaster ======================"; + // QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器 + // QByteArray byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray + //,这个应该会自动释放 const char* charArray = byteArray.constData(); // + // 获取QByteArray的const char*指针 + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); + // std::cout<GetGeoTransform(adfGeoTransform) != CE_None) { + qDebug() << "Failed to get GeoTransform"; + return nullptr; + } + + double dfWidth = poDataset->GetRasterXSize(); + double dfHeight = poDataset->GetRasterYSize(); + + // 计算中心点坐标（像素坐标） + double dfCenterX = adfGeoTransform[0] + dfWidth * adfGeoTransform[1] / 2.0 + + dfHeight * adfGeoTransform[2] / 2.0; + double dfCenterY = adfGeoTransform[3] + dfWidth * adfGeoTransform[4] / 2.0 + + dfHeight * adfGeoTransform[5] / 2.0; + + OGRSpatialReference oSRS; + oSRS.importFromWkt(poDataset->GetProjectionRef()); + + if(oSRS.IsGeographic()) { + qDebug() << "Center coords (already in geographic): (" + QString::number(dfCenterX) + + ", " + QString::number(dfCenterY) + ")"; + flag = true; + x = dfCenterX; + y = dfCenterY; + } else { + // 如果不是地理坐标系，转换到WGS84 + OGRSpatialReference oSRS_WGS84; + oSRS_WGS84.SetWellKnownGeogCS("WGS84"); + + OGRCoordinateTransformation* poCT = + OGRCreateCoordinateTransformation(&oSRS, &oSRS_WGS84); + if(poCT == nullptr) { + qDebug() << "Failed to create coordinate transformation"; + return nullptr; + } + + // double dfLon, dfLat; + if(poCT->Transform(1, &dfCenterX, &dfCenterY)) { + qDebug() << "Center coords (transformed to WGS84): (" + + QString::number(dfCenterX) + ", " + QString::number(dfCenterY) + << ")"; + flag = true; + x = dfCenterX; + y = dfCenterY; + } else { + qDebug() << "Transformation failed."; + } + OGRCoordinateTransformation::DestroyCT(poCT); + } + } + if(nullptr==poDataset||NULL==poDataset){}else{ + GDALClose(poDataset); + } + + if(flag) { + std::shared_ptr RasterCenterPoint = std::make_shared(); + RasterCenterPoint->x = x; + RasterCenterPoint->y = y; + RasterCenterPoint->z = 0; + return RasterCenterPoint; + } else { + return nullptr; + } + } + CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long epsg_code) + { + OGRSpatialReference oSRS; + if(oSRS.importFromEPSG(epsg_code) == OGRERR_NONE) { + if(oSRS.IsGeographic()) { + return CoordinateSystemType::GeoCoordinateSystem; + } else if(oSRS.IsProjected()) { + return CoordinateSystemType::ProjectCoordinateSystem; + } + } else { + return CoordinateSystemType::UNKNOW; + } + } +} // namespace RasterToolBase \ No newline at end of file diff --git a/RasterToolBase.h b/RasterToolBase.h new file mode 100644 index 0000000..e042beb --- /dev/null +++ b/RasterToolBase.h @@ -0,0 +1,79 @@ +/** + * @file RasterProjectBase.h + * @brief None + * @author 陈增辉 (3045316072@qq.com) + * @version 2.5.0 + * @date 24-6-4 + * @copyright Copyright (c) Since 2024 中科卫星应用研究院 All rights reserved. + */ + +#ifndef LAMPCAE_RASTERTOOLBASE_H +#define LAMPCAE_RASTERTOOLBASE_H +#include "gdal_priv.h" +#include + +namespace RasterToolBase { + + static bool GDALAllRegisterEnable=false; + + + enum ProjectStripDelta{ + Strip_6, // 6度带 + Strip_3 + }; + + enum CoordinateSystemType{ // 坐标系类型 + GeoCoordinateSystem, + ProjectCoordinateSystem, + UNKNOW + }; + + struct PointRaster{ // 影像坐标点 + double x; + double y; + double z; + }; + + + struct PointXYZ{ + double x,y,z; + }; + + struct PointGeo{ + double lon,lat,ati; + }; + + struct PointImage{ + double pixel_x,pixel_y; + }; + + /// 根据经纬度获取 + /// EPSG代码，根据经纬度返回对应投影坐标系统，其中如果在中华人民共和国境内，默认使用 + /// CGCS2000坐标系统(EPSG 4502 ~ 4512 6度带,EPSG 4534 ~ 4554 3度带)，其余地方使用WGS坐标系统， + /// 投影方法 高斯克吕格（国内）， 高斯克吕格 + /// \param long 经度 + /// \param lat 纬度 + /// \return 对应投影坐标系统的 EPSG编码,-1 表示计算错误 + long getProjectEPSGCodeByLon_Lat(double long, double lat, + ProjectStripDelta stripState = ProjectStripDelta::Strip_3); + + long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat); + + long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat); + + + QString GetProjectionNameFromEPSG(long epsgCode) ; + + + long GetEPSGFromRasterFile(QString filepath); + + std::shared_ptr GetCenterPointInRaster(QString filepath); + + CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long EPSGCODE); + + + + +} // namespace RasterProjectConvertor + +#endif // LAMPCAE_RASTERTOOLBASE_H