KJ135
/
BaseLibraryCPP
mirror of http://172.16.0.12:5000/LAMPSARToolSoftware/BaseLibraryCPP.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

更新影像拼接代码

main
陈增辉 2024-11-15 14:21:50 +08:00
parent 45c32c0b40
commit 23277b911e
11 changed files with 1290 additions and 204 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
BaseConstVariable.h
View File

@ -51,6 +51,9 @@ const double earth_Re = 6378136.49;
const double earth_Rp = (1 - 1 / f_inverse) * earth_Re; const double earth_Rp = (1 - 1 / f_inverse) * earth_Re;
const double earth_We = 0.000072292115; const double earth_We = 0.000072292115;
/*********************************************** openMap参数 ********************************************************************/
static long Paral_num_thread = 6;
/*********************************************** 基础枚举 ********************************************************************/ /*********************************************** 基础枚举 ********************************************************************/

1
BaseTool.cpp
View File

@ -115,7 +115,6 @@ std::complex<double> Cubic_kernel_weight(double s)
/// <returns></returns> /// <returns></returns>
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)
{ {
return p11.ati * (1 - p0.lon) * (1 - p0.lat) + return p11.ati * (1 - p0.lon) * (1 - p0.lat) +
p12.ati * p0.lon * (1 - p0.lat) + p12.ati * p0.lon * (1 - p0.lat) +
p21.ati * (1 - p0.lon) * p0.lat + p21.ati * (1 - p0.lon) * p0.lat +

5
BaseTool.h
View File

@ -26,8 +26,6 @@
#include <string> #include <string>
#include <QString> #include <QString>
#include <QStringList> #include <QStringList>
#include "LogInfoCls.h" #include "LogInfoCls.h"
///////////////////////////////////// 运行时间打印 ///////////////////////////////////// 运行时间打印
@ -56,8 +54,7 @@ std::complex<double> Cubic_Convolution_interpolation(double u, double v,
std::complex<double> Cubic_kernel_weight(double s); std::complex<double> Cubic_kernel_weight(double s);
double Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12, double Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12,Landpoint p22);
Landpoint p22);
bool onSegment(Point3 Pi, Point3 Pj, Point3 Q); bool onSegment(Point3 Pi, Point3 Pj, Point3 Q);

23
FileOperator.cpp
View File

@ -13,6 +13,29 @@
#include <QString> #include <QString>
#include <QMessageBox> #include <QMessageBox>
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<QString> getFilelist(const QString& folderpath, const QString& filenameExtension, int (*logfun)(QString logtext, int value)) std::vector<QString> getFilelist(const QString& folderpath, const QString& filenameExtension, int (*logfun)(QString logtext, int value))
{ {
QString filenameExtensionStr = filenameExtension; QString filenameExtensionStr = filenameExtension;

1
FileOperator.h
View File

@ -45,5 +45,6 @@ size_t fsize(FILE* fp);
QString getParantFromPath(const QString& path); QString getParantFromPath(const QString& path);
void copyFile(const QString& sourcePath, const QString& destinationPath); void copyFile(const QString& sourcePath, const QString& destinationPath);
QString addMaskToFileName(const QString& filePath, QString _zzname);
// QT FileOperator // QT FileOperator
#endif #endif

834
ImageOperatorBase.cpp
View File

@ -569,24 +569,28 @@ Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count,
} }
delete[] temp; delete[] temp;
} }
// else if (gdal_datatype == GDT_UInt64) { else if (gdal_datatype == GDT_UInt64) {
// unsigned long* temp = new unsigned long[rows_count * cols_count]; 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, 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 < rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) {
//cols_count; j++) { datamatrix(i, j) = temp[i * cols_count + j]; for (int j = 0; j <
// } cols_count; j++) {
// } datamatrix(i, j) = temp[i * cols_count + j];
// delete[] temp; }
// } }
// else if (gdal_datatype == GDT_Int64) { delete[] temp;
// long* temp = new long[rows_count * cols_count]; }
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, else if (gdal_datatype == GDT_Int64) {
//rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { for (int j = 0; j < long* temp = new long[rows_count * cols_count];
//cols_count; j++) { datamatrix(i, j) = temp[i * cols_count + j]; 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 <
// delete[] temp; cols_count; j++) {
// } datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
else if(gdal_datatype == GDT_Float32) { else if(gdal_datatype == GDT_Float32) {
float* temp = new float[rows_count * cols_count]; float* temp = new float[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, 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; 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() Eigen::MatrixXd gdalImage::getGeoTranslation()
{ {
return this->gt; return this->gt;
@ -695,6 +827,61 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col
omp_destroy_lock(&lock); // 劫伙拷斤拷 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() void gdalImage::saveImage()
{ {
this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids); 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; 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, gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num,
Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite) 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; 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, gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width,
int band_num, Eigen::MatrixXd gt, QString projection, int band_num, Eigen::MatrixXd gt, QString projection,
bool need_gt, bool overwrite) bool need_gt, bool overwrite)
@ -1275,6 +1544,275 @@ int saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path)
return -1; return -1;
} }
ErrorCode MergeRasterProcess(QVector<QString> 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<gdalImage> 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<RasterExtend> 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<gdalImage> imgdslist, gdalImage resultimg, gdalImage maskimg)
{
omp_set_num_threads(Paral_num_thread);
// 逐点合并计算
QVector<RasterExtend> 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) gdalImageComplex::gdalImageComplex(const QString& raster_path)
{ {
omp_lock_t lock; 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<int>((lon + 180.0) / 6.0) + 1;
bool isSouth = lon < 0; // 判断是否在南半球
long epsgCodeBase = isSouth ? 32700 : 32600; // 计算EPSG代码
long prefix = static_cast<int>(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<<filepath.toLocal8Bit().constData()<<std::endl;
// 打开影像文件
GDALDataset* poDataset;
poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
if (NULL == poDataset) {
qDebug() << "Error: Unable to open the image.\n";
return -1;
}
// 获取影像的投影信息
const char* pszProjection = poDataset->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<PointRaster> 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<<filepath.toLocal8Bit().constData()<<std::endl;
GDALDataset* poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
if (nullptr == poDataset || NULL == poDataset) {
qDebug() << "Could not open dataset";
return nullptr; // 表示读取失败
}
double x, y, z;
bool flag = false;
{
double adfGeoTransform[6];
if (poDataset->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<PointRaster> RasterCenterPoint = std::make_shared<PointRaster>();
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;
}
}

268
ImageOperatorBase.h
View File

@ -24,6 +24,39 @@
#include <memory> #include <memory>
#include <QString> #include <QString>
#include <cpl_conv.h> // for CPLMalloc() #include <cpl_conv.h> // 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 { struct ImageGEOINFO {
@ -40,8 +73,39 @@ enum GDALREADARRCOPYMETHOD {
}; };
// 判断是否需要输出为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<PointRaster> GetCenterPointInRaster(QString filepath);
CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long EPSGCODE);
// 文件打开 // 当指令销毁时，调用GDALClose 销毁类型 // 文件打开 // 当指令销毁时，调用GDALClose 销毁类型
std::shared_ptr<GDALDataset> OpenDataset(const QString& in_path, GDALAccess rwmode= GA_ReadOnly); std::shared_ptr<GDALDataset> OpenDataset(const QString& in_path, GDALAccess rwmode= GA_ReadOnly);
void CloseDataset(GDALDataset* ptr); void CloseDataset(GDALDataset* ptr);
@ -70,11 +134,11 @@ ImageGEOINFO getImageINFO(QString in_path);
GDALDataType getGDALDataType(QString fileptah); GDALDataType getGDALDataType(QString fileptah);
struct DemBox { struct RasterExtend {
double min_lat; //纬度 double min_x; //纬度
double min_lon;//经度 double min_y;//经度
double max_lat;//纬度 double max_x;//纬度
double max_lon;//经度 double max_y;//经度
}; };
@ -97,9 +161,12 @@ public: // 方法
virtual void setTranslationMatrix(Eigen::MatrixXd gt); virtual void setTranslationMatrix(Eigen::MatrixXd gt);
virtual void setData(Eigen::MatrixXd,int data_band_ids=1); 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::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 Eigen::MatrixXd getGeoTranslation();
virtual GDALDataType getDataType(); virtual GDALDataType getDataType();
virtual void saveImage(Eigen::MatrixXd, int start_row, int start_col, int band_ids); 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 saveImage();
virtual void setNoDataValue(double nodatavalue, int band_ids); virtual void setNoDataValue(double nodatavalue, int band_ids);
virtual int InitInv_gt(); virtual int InitInv_gt();
@ -116,6 +183,7 @@ public: // 方法
virtual Eigen::MatrixXd getHist(int bandids); virtual Eigen::MatrixXd getHist(int bandids);
virtual RasterExtend getExtend();
public: public:
QString img_path; // 图像文件 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, 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 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); 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 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); 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<QString> filepath, QString outfileptah, QString mainString, MERGEMODE mergecode = MERGEMODE::MERGE_GEOCODING, bool isENVI = false);
ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> inimgs, gdalImage resultimg,gdalImage maskimg);
template<typename T> template<typename T>
std::shared_ptr<T> readDataArr(gdalImage &img,int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method); std::shared_ptr<T> 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
#endif
template<typename T>
std::shared_ptr<T> readDataArr(gdalImage& imgds, int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method)
{
std::shared_ptr<T> 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<T>(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<T>(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<T>(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<T>(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<T>(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<T>(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<T>(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;
}

4
LogInfoCls.cpp
View File

@ -18,6 +18,10 @@ std::string errorCode2errInfo(ErrorCode e)
_CASE_STR(FILE_NOT_EXIST); _CASE_STR(FILE_NOT_EXIST);
_CASE_STR(FIND_ID_ERROR); _CASE_STR(FIND_ID_ERROR);
_CASE_STR(INSERT_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 //GSL 1xx
_CASE_STR(Error_GSL_FAILURE ); _CASE_STR(Error_GSL_FAILURE );
_CASE_STR(Error_GSL_CONTINUE ); _CASE_STR(Error_GSL_CONTINUE );

4
LogInfoCls.h
View File

@ -27,6 +27,10 @@ enum ErrorCode {
FILE_NOT_EXIST = 9, FILE_NOT_EXIST = 9,
FIND_ID_ERROR = 10, FIND_ID_ERROR = 10,
INSERT_ID_ERROR = 11, INSERT_ID_ERROR = 11,
EPSGCODE_NOTSAME = 12,
EPSGCODE_NOTSUPPORT = 13,
RASTERBAND_NOTEQUAL = 14,
RASTER_DATETYPE_NOTSAME = 15,
//GSL 1xx //GSL 1xx
Error_GSL_FAILURE = -101, Error_GSL_FAILURE = -101,
Error_GSL_CONTINUE = -102, /* iteration has not converged */ Error_GSL_CONTINUE = -102, /* iteration has not converged */

268
RasterToolBase.cpp Normal file
View File

@ -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 <QDebug>
#include "RasterToolBase.h"
#include "gdal_priv.h"
#include "ogr_spatialref.h"
#include "cpl_conv.h" // for CPLMalloc()
#include <QTextCodec>
#include <iostream>
#include <QFile>
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<int>((lon + 180.0) / 6.0) + 1;
bool isSouth = lon < 0; // 判断是否在南半球
long epsgCodeBase = isSouth ? 32700 : 32600; // 计算EPSG代码
long prefix = static_cast<int>(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<<filepath.toLocal8Bit().constData()<<std::endl;
// 打开影像文件
GDALDataset* poDataset;
poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
if(NULL==poDataset) {
qDebug() << "Error: Unable to open the image.\n";
return -1;
}
// 获取影像的投影信息
const char* pszProjection = poDataset->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<PointRaster> 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<<filepath.toLocal8Bit().constData()<<std::endl;
GDALDataset* poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
if(nullptr==poDataset||NULL==poDataset) {
qDebug() << "Could not open dataset";
return nullptr; // 表示读取失败
}
double x, y, z;
bool flag = false;
{
double adfGeoTransform[6];
if(poDataset->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<PointRaster> RasterCenterPoint = std::make_shared<PointRaster>();
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

79
RasterToolBase.h Normal file
View File

@ -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 <memory>
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<PointRaster> GetCenterPointInRaster(QString filepath);
CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long EPSGCODE);
} // namespace RasterProjectConvertor
#endif // LAMPCAE_RASTERTOOLBASE_H