#include "ImageOperatorBase.h" #include "GeoOperator.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "FileOperator.h" using namespace std; using namespace Eigen; /** * 输入数据是ENVI格式数据 */ std::shared_ptr OpenDataset(const QString& in_path,GDALAccess rwmode) { GDALAllRegister(); GDALDataset* dataset_ptr = (GDALDataset*)(GDALOpen(in_path.toUtf8().constData(), rwmode)); std::shared_ptr rasterDataset(dataset_ptr, CloseDataset); return rasterDataset; } void CloseDataset(GDALDataset* ptr) { GDALClose(ptr); ptr = NULL; } int TIFF2ENVI(QString in_tiff_path, QString out_envi_path) { std::shared_ptr ds = OpenDataset(in_tiff_path); const char* args[] = { "-of", "ENVI", NULL }; GDALTranslateOptions* psOptions = GDALTranslateOptionsNew((char**)args, NULL); GDALClose(GDALTranslate(out_envi_path.toUtf8().constData(),ds.get(), psOptions,NULL)); GDALTranslateOptionsFree(psOptions); return 0; } int ENVI2TIFF(QString in_envi_path, QString out_tiff_path) { std::shared_ptr ds = OpenDataset(in_envi_path); const char* args[] = { "-of", "Gtiff", NULL }; GDALTranslateOptions* psOptions = GDALTranslateOptionsNew((char**)args, NULL); GDALClose(GDALTranslate(out_tiff_path.toUtf8().constData(), ds.get(), psOptions, NULL)); GDALTranslateOptionsFree(psOptions); return 0; } int CreateDataset(QString new_file_path,int height, int width, int band_num,double* gt, QString projection, GDALDataType gdal_dtype ,bool need_gt) { GDALAllRegister(); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); std::shared_ptr< GDALDataset> poDstDS(poDriver->Create(new_file_path.toUtf8().constData(), width, height, band_num, gdal_dtype, NULL)); if (need_gt) { poDstDS->SetProjection(projection.toUtf8().constData()); poDstDS->SetGeoTransform(gt); } else {} GDALFlushCache((GDALDatasetH)poDstDS.get()); return 0; } int saveDataset(QString new_file_path, int start_line,int start_cols ,int band_ids, int datacols,int datarows,void* databuffer) { GDALAllRegister(); std::shared_ptr poDstDS = OpenDataset(new_file_path,GA_Update); GDALDataType gdal_datatype = poDstDS->GetRasterBand(1)->GetRasterDataType(); poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_cols, start_line, datacols, datarows, databuffer, datacols, datarows, gdal_datatype, 0, 0); GDALFlushCache(poDstDS.get()); return 0; } int block_num_pre_memory(int block_width, int height, GDALDataType gdal_datatype,double memey_size) { // 计算大小 int size_meta = 0; if (gdal_datatype == GDT_Byte) { size_meta = 1; } else if (gdal_datatype == GDT_UInt16) { size_meta = 2; // 只有双通道才能构建 复数矩阵 } else if (gdal_datatype == GDT_UInt16) { size_meta = 2; } else if (gdal_datatype == GDT_Int16) { size_meta = 2; } else if (gdal_datatype == GDT_UInt32) { size_meta = 4; } else if (gdal_datatype == GDT_Int32) { size_meta = 4; } //else if (gdal_datatype == GDT_UInt64) { // size_meta = 8; //} //else if (gdal_datatype == GDT_Int64) { // size_meta = 8; //} else if (gdal_datatype == GDT_Float32) { size_meta = 4; } else if (gdal_datatype == GDT_Float64) { size_meta = 4; } else if (gdal_datatype == GDT_CInt16) { size_meta = 2; } else if (gdal_datatype == GDT_CInt32) { size_meta = 2; } else if (gdal_datatype == GDT_CFloat32) { size_meta = 4; } else if (gdal_datatype == GDT_CFloat64) { size_meta = 8; } else {} int block_num = int(memey_size / (size_meta * block_width)); block_num = block_num > height ? height : block_num; // 行数 block_num = block_num < 1 ? 1 : block_num; return block_num; } Eigen::Matrix ReadComplexMatrixData(int start_line, int width, int line_num, std::shared_ptr rasterDataset, GDALDataType gdal_datatype) { int band_num = rasterDataset->GetRasterCount(); if (gdal_datatype == 0) { return Eigen::Matrix (0, 0); } else if (gdal_datatype < 8) { if (band_num != 2) { return Eigen::Matrix(0, 0); } } else if (gdal_datatype < 12) { if (band_num != 1) { return Eigen::Matrix(0, 0); } } else {} bool _flag = false; Eigen::Matrix data_mat(line_num * width, 2);// 必须强制行优先 if (gdal_datatype == GDT_Byte) { Eigen::MatrixX real_mat(line_num * width, 1); Eigen::MatrixX imag_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag data_mat.col(0) = (real_mat.array().cast()).array(); data_mat.col(1) = (imag_mat.array().cast()).array(); _flag = true; } else if (gdal_datatype == GDT_UInt16) { Eigen::MatrixX real_mat(line_num * width, 1); Eigen::MatrixX imag_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag data_mat.col(0) = (real_mat.array().cast()).array(); data_mat.col(1) = (imag_mat.array().cast()).array(); _flag = true; } else if (gdal_datatype == GDT_Int16) { Eigen::MatrixX real_mat(line_num * width, 1); Eigen::MatrixX imag_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag data_mat.col(0) = (real_mat.array().cast()).array(); data_mat.col(1) = (imag_mat.array().cast()).array(); _flag = true; } else if (gdal_datatype == GDT_UInt32) { Eigen::MatrixX real_mat(line_num * width, 1); Eigen::MatrixX imag_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag data_mat.col(0) = (real_mat.array().cast()).array(); data_mat.col(1) = (imag_mat.array().cast()).array(); _flag = true; } else if (gdal_datatype == GDT_Int32) { Eigen::MatrixX real_mat(line_num * width, 1); Eigen::MatrixX imag_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag data_mat.col(0) = (real_mat.array().cast()).array(); data_mat.col(1) = (imag_mat.array().cast()).array(); _flag = true; } //else if (gdal_datatype == GDT_UInt64) { // Eigen::MatrixX real_mat(line_num * width, 1); // Eigen::MatrixX imag_mat(line_num * width, 1); // rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real // rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag // data_mat.col(0) = (real_mat.array().cast()).array(); // data_mat.col(1) = (imag_mat.array().cast()).array(); // _flag = true; //} //else if (gdal_datatype == GDT_Int64) { // Eigen::MatrixX real_mat(line_num * width, 1); // Eigen::MatrixX imag_mat(line_num * width, 1); // rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real // rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag // data_mat.col(0) = (real_mat.array().cast()).array(); // data_mat.col(1) = (imag_mat.array().cast()).array(); // _flag = true; //} else if (gdal_datatype == GDT_Float32) { Eigen::MatrixX real_mat(line_num * width, 1); Eigen::MatrixX imag_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag data_mat.col(0) = (real_mat.array().cast()).array(); data_mat.col(1) = (imag_mat.array().cast()).array(); _flag = true; } else if (gdal_datatype == GDT_Float64) { Eigen::MatrixX real_mat(line_num * width, 1); Eigen::MatrixX imag_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag data_mat.col(0) = (real_mat.array().cast()).array(); data_mat.col(1) = (imag_mat.array().cast()).array(); _flag = true; } else if (gdal_datatype == GDT_CInt16) { Eigen::MatrixX> complex_short_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, complex_short_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = (complex_short_mat.real().array().cast()).array(); data_mat.col(1) = (complex_short_mat.imag().array().cast()).array(); _flag = true; } else if (gdal_datatype == GDT_CInt32) { Eigen::MatrixX> complex_short_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, complex_short_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = (complex_short_mat.real().array().cast()).array(); data_mat.col(1) = (complex_short_mat.imag().array().cast()).array(); _flag = true; } else if (gdal_datatype == GDT_CFloat32) { Eigen::MatrixX> complex_short_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, complex_short_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = (complex_short_mat.real().array().cast()).array(); data_mat.col(1) = (complex_short_mat.imag().array().cast()).array(); _flag = true; } else if (gdal_datatype == GDT_CFloat64) { Eigen::MatrixX> complex_short_mat(line_num * width, 1); rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, complex_short_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = (complex_short_mat.real().array().cast()).array(); data_mat.col(1) = (complex_short_mat.imag().array().cast()).array(); _flag = true; } else {} // 保存数据 if (_flag) { return data_mat; } else { return Eigen::Matrix(0, 0);// 必须强制行优先; } } Eigen::Matrix ReadMatrixDoubleData(int start_line, int width, int line_num, std::shared_ptr rasterDataset, GDALDataType gdal_datatype, int band_idx) { // 构建矩阵块,使用eigen 进行矩阵计算,加速计算 bool _flag = false; Eigen::Matrix data_mat(line_num * width, 1);// 必须强制行优先 if (gdal_datatype == GDT_Byte) { Eigen::MatrixX real_mat(line_num * width, 1); rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; _flag = true; } else if (gdal_datatype == GDT_UInt16) { Eigen::MatrixX real_mat(line_num * width, 1); rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; _flag = true; } else if (gdal_datatype == GDT_Int16) { Eigen::MatrixX real_mat(line_num * width, 1); rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; _flag = true; } else if (gdal_datatype == GDT_UInt32) { Eigen::MatrixX real_mat(line_num * width, 1); rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; _flag = true; } else if (gdal_datatype == GDT_Int32) { Eigen::MatrixX real_mat(line_num * width, 1); rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; _flag = true; } //else if (gdal_datatype == GDT_UInt64) { // Eigen::MatrixX real_mat(line_num * width, 1); // rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real // data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; // _flag = true; //} //else if (gdal_datatype == GDT_Int64) { // Eigen::MatrixX real_mat(line_num * width, 1); // rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real // data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; // _flag = true; //} else if (gdal_datatype == GDT_Float32) { Eigen::MatrixX real_mat(line_num * width, 1); rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; _flag = true; } else if (gdal_datatype == GDT_Float64) { Eigen::MatrixX real_mat(line_num * width, 1); rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; _flag = true; } else {} return data_mat; } Eigen::MatrixXd getGeoTranslationArray(QString in_path) { return Eigen::MatrixXd(); } ImageGEOINFO getImageINFO(QString in_path) { std::shared_ptr df = OpenDataset(in_path); int width = df->GetRasterXSize(); int heigh = df->GetRasterYSize(); int band_num = df->GetRasterCount(); ImageGEOINFO result; result.width = width; result.height = heigh; result.bandnum = band_num; return result; } GDALDataType getGDALDataType(QString fileptah) { omp_lock_t lock; omp_init_lock(&lock); omp_set_lock(&lock); GDALAllRegister(); GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(fileptah.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); GDALClose((GDALDatasetH)rasterDataset); omp_unset_lock(&lock); //锟酵放伙拷斤拷 omp_destroy_lock(&lock); //劫伙拷斤拷 return gdal_datatype; } gdalImage::gdalImage() { this->height = 0; this->width = 0; this->data_band_ids = 1; this->start_row = 0; this->start_col = 0; } /// /// 斤拷图饺∮帮拷锟?1锟?7 /// /// gdalImage::gdalImage(const QString& raster_path) { omp_lock_t lock; omp_init_lock(&lock); // 锟斤拷始斤拷斤拷 omp_set_lock(&lock); //锟斤拷没斤拷锟?1锟?7 this->img_path = raster_path; GDALAllRegister();// 注绞斤拷斤拷锟?1锟?7 // 锟斤拷DEM影锟斤拷 GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 this->width = rasterDataset->GetRasterXSize(); this->height = rasterDataset->GetRasterYSize(); this->band_num = rasterDataset->GetRasterCount(); double* gt = new double[6]; // 锟斤拷梅斤拷锟斤拷 rasterDataset->GetGeoTransform(gt); this->gt = Eigen::MatrixXd(2, 3); this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]; this->projection = rasterDataset->GetProjectionRef(); // 斤拷斤拷 //double* inv_gt = new double[6];; //GDALInvGeoTransform(gt, inv_gt); // 斤拷斤拷 // 斤拷投影 GDALFlushCache((GDALDatasetH)rasterDataset); GDALClose((GDALDatasetH)rasterDataset); rasterDataset = NULL;// 指矫匡拷 this->InitInv_gt(); delete[] gt; ////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH omp_unset_lock(&lock); //锟酵放伙拷斤拷 omp_destroy_lock(&lock); //劫伙拷斤拷 } gdalImage::~gdalImage() { } void gdalImage::setHeight(int height) { this->height = height; } void gdalImage::setWidth(int width) { this->width = width; } void gdalImage::setTranslationMatrix(Eigen::MatrixXd gt) { this->gt = gt; } void gdalImage::setData(Eigen::MatrixXd, int data_band_ids) { this->data = data; this->data_band_ids = data_band_ids; } Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids = 1) { omp_lock_t lock; omp_init_lock(&lock); omp_set_lock(&lock); GDALAllRegister(); 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; 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); 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; } GDALDataType gdalImage::getDataType() { GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly)); GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); return gdal_datatype; } /// /// 写锟斤拷遥锟斤拷影锟斤拷 /// /// /// /// /// void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0, int band_ids = 1) { 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; qDebug() << tip; throw exception(tip.toUtf8().constData()); } GDALAllRegister(); 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; } int datarows = data.rows(); int datacols = data.cols(); float* databuffer = new float[datarows * datacols];// (float*)malloc(datarows * datacols * sizeof(float)); for (int i = 0; i < datarows; i++) { for (int j = 0; j < datacols; j++) { float temp = float(data(i, j)); databuffer[i * datacols + j] = temp; } } //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_Float32, 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); } void gdalImage::setNoDataValue(double nodatavalue = -9999, int band_ids = 1) { GDALAllRegister();// 注绞斤拷斤拷锟?1锟?7 //GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update)); poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue); GDALFlushCache((GDALDatasetH)poDstDS); GDALClose((GDALDatasetH)poDstDS); } int gdalImage::InitInv_gt() { //1 lon lat = x //1 lon lat = y Eigen::MatrixXd temp=Eigen::MatrixXd::Zero(2, 3); this->inv_gt = temp; double a = this->gt(0, 0); double b = this->gt(0, 1); double c = this->gt(0, 2); double d = this->gt(1, 0); double e = this->gt(1, 1); double f = this->gt(1, 2); double g = 1; double det_gt = b * f - c * e; if (det_gt == 0) { return 0; } this->inv_gt(0, 0) = (c * d - a * f) / det_gt; //2 this->inv_gt(0, 1) = f / det_gt; // lon this->inv_gt(0, 2) = -c / det_gt; // lat this->inv_gt(1, 0) = (a * e - b * d) / det_gt; //1 this->inv_gt(1, 1) = -e / det_gt; // lon this->inv_gt(1, 2) = b / det_gt; // lat return 1; } Landpoint gdalImage::getRow_Col(double lon, double lat) { Landpoint p{ 0,0,0 }; p.lon = this->inv_gt(0, 0) + lon * this->inv_gt(0, 1) + lat * this->inv_gt(0, 2); //x p.lat = this->inv_gt(1, 0) + lon * this->inv_gt(1, 1) + lat * this->inv_gt(1, 2); //y return p; } Landpoint gdalImage::getLandPoint(double row, double col, double ati = 0) { Landpoint p{ 0,0,0 }; p.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); //x p.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); //y p.ati = ati; return p; } double gdalImage::mean(int bandids) { double mean_value = 0; double count = this->height * this->width; int line_invert = 100; int start_ids = 0; do { Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids); mean_value = mean_value + sar_a.sum() / count; start_ids = start_ids + line_invert; } while (start_ids < this->height); return mean_value; } double gdalImage::max(int bandids) { double max_value = 0; bool state_max = true; int line_invert = 100; int start_ids = 0; double temp_max = 0; do { Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids); if (state_max) { state_max = false; max_value = sar_a.maxCoeff(); } else { temp_max = sar_a.maxCoeff(); if (max_value < temp_max) { max_value = temp_max; } } start_ids = start_ids + line_invert; } while (start_ids < this->height); return max_value; } double gdalImage::min(int bandids) { double min_value = 0; bool state_min = true; int line_invert = 100; int start_ids = 0; double temp_min = 0; do { Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids); if (state_min) { state_min = false; min_value = sar_a.minCoeff(); } else { temp_min = sar_a.minCoeff(); if (min_value < temp_min) { min_value = temp_min; } } start_ids = start_ids + line_invert; } while (start_ids < this->height); return min_value; } GDALRPCInfo gdalImage::getRPC() { CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); CPLSetConfigOption("GDAL_DATA", "./data"); GDALAllRegister();//注斤拷锟斤拷 //斤拷锟斤拷 GDALDataset* pDS = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly); //锟斤拷元斤拷锟叫伙拷取RPC锟斤拷息 char** papszRPC = pDS->GetMetadata("RPC"); //斤拷取锟斤拷RPC锟斤拷息斤拷山峁癸拷锟?1锟?7 GDALRPCInfo oInfo; GDALExtractRPCInfo(papszRPC, &oInfo); GDALClose((GDALDatasetH)pDS); return oInfo; } Eigen::MatrixXd gdalImage::getLandPoint(Eigen::MatrixXd points) { if (points.cols() != 3) { throw new exception("the size of points is equit 3!!!"); } Eigen::MatrixXd result(points.rows(), 3); result.col(2) = points.col(2);// 锟竭筹拷 points.col(2) = points.col(2).array() * 0 + 1; points.col(0).swap(points.col(2));// 斤拷 Eigen::MatrixXd gts(3, 2); gts.col(0) = this->gt.row(0); gts.col(1) = this->gt.row(1); result.block(0, 0, points.rows(), 2) = points * gts; return result; } Eigen::MatrixXd gdalImage::getHist(int bandids) { GDALAllRegister();// 注绞斤拷斤拷锟?1锟?7 // 锟斤拷DEM影锟斤拷 GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); GDALRasterBand* xBand = rasterDataset->GetRasterBand(bandids); double dfMin = this->min(bandids); double dfMax = this->max(bandids); int count = int((dfMax - dfMin) / 0.01); count = count > 255 ? count : 255; GUIntBig* panHistogram = new GUIntBig[count]; xBand->GetHistogram(dfMin, dfMax, count, panHistogram, TRUE, FALSE, NULL, NULL); Eigen::MatrixXd result(count, 2); double delta = (dfMax - dfMin) / count; for (int i = 0; i < count; i++) { result(i, 0) = dfMin + i * delta; result(i, 1) = double(panHistogram[i]); } delete[] panHistogram; GDALClose((GDALDatasetH)rasterDataset); return result; } gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt,bool overwrite) { 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();// 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7 GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_Float32, NULL); // 锟斤拷锟斤拷锟斤拷 if (need_gt) { poDstDS->SetProjection(projection.toUtf8().constData()); // 锟斤拷锟斤拷转锟斤拷锟斤拷锟斤拷 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) { if (exists_test(img_path.toUtf8().constData())) { if (overwrite) { gdalImageComplex result_img(img_path); return result_img; } else { throw "file has exist!!!"; exit(1); } } GDALAllRegister(); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_CFloat64, NULL); if (need_gt) { poDstDS->SetProjection(projection.toUtf8().constData()); // 锟斤拷锟斤拷转锟斤拷锟斤拷锟斤拷 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 gdalImageComplex result_img(img_path); return result_img; } int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample) { GDALAllRegister(); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly); if (pDSrc == NULL) { return -1; } GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); if (pDriver == NULL) { GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); return -2; } int width = pDSrc->GetRasterXSize(); int height = pDSrc->GetRasterYSize(); int nBandCount = pDSrc->GetRasterCount(); GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType(); char* pszSrcWKT = NULL; pszSrcWKT = const_cast(pDSrc->GetProjectionRef()); //锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7 if (strlen(pszSrcWKT) <= 0) { OGRSpatialReference oSRS; oSRS.importFromEPSG(4326); //oSRS.SetUTM(50, true); //锟斤拷锟斤拷锟斤拷 锟斤拷锟斤拷120锟斤拷 //oSRS.SetWellKnownGeogCS("WGS84"); oSRS.exportToWkt(&pszSrcWKT); } qDebug() << "GDALCreateGenImgProjTransformer " << endl; void* hTransformArg; hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT, FALSE, 0.0, 1); qDebug() << "no proj " << endl; //(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7) if (hTransformArg == NULL) { GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); return -3; } qDebug() << "has proj " << endl; double dGeoTrans[6] = { 0 }; int nNewWidth = 0, nNewHeight = 0; if (GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg, dGeoTrans, &nNewWidth, &nNewHeight) != CE_None) { GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); return -3; } //GDALDestroyGenImgProjTransformer(hTransformArg); GDALDataset* pDDst = pDriver->Create(pszOutFile, new_width, new_height, nBandCount, dataType, NULL); if (pDDst == NULL) { GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); return -2; } pDDst->SetProjection(pszSrcWKT); pDDst->SetGeoTransform(gt); GDALWarpOptions* psWo = GDALCreateWarpOptions(); //psWo->papszWarpOptions = CSLDuplicate(NULL); psWo->eWorkingDataType = dataType; psWo->eResampleAlg = eResample; psWo->hSrcDS = (GDALDatasetH)pDSrc; psWo->hDstDS = (GDALDatasetH)pDDst; qDebug() << "GDALCreateGenImgProjTransformer" << endl; psWo->pfnTransformer = GDALGenImgProjTransform; psWo->pTransformerArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);; qDebug() << "GDALCreateGenImgProjTransformer has created" << endl; psWo->nBandCount = nBandCount; psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int)); psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int)); for (int i = 0; i < nBandCount; i++) { psWo->panSrcBands[i] = i + 1; psWo->panDstBands[i] = i + 1; } GDALWarpOperation oWo; if (oWo.Initialize(psWo) != CE_None) { GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); GDALClose((GDALDatasetH)(GDALDatasetH)pDDst); return -3; } qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << endl; oWo.ChunkAndWarpMulti(0, 0, new_width, new_height); GDALFlushCache(pDDst); qDebug() << "ChunkAndWarpImage over" << endl; //GDALDestroyGenImgProjTransformer(psWo->pTransformerArg); //GDALDestroyWarpOptions(psWo); GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); GDALClose((GDALDatasetH)(GDALDatasetH)pDDst); ////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH return 0; } int ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResample) { GDALAllRegister(); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_Update); if (pDSrc == NULL) { return -1; } GDALDataType gdal_datatype = pDSrc->GetRasterBand(1)->GetRasterDataType(); GDALRasterBand* demBand = pDSrc->GetRasterBand(band_ids); int width = pDSrc->GetRasterXSize(); int height = pDSrc->GetRasterYSize(); int start_col = 0, start_row = 0, rows_count = 0, cols_count; int row_delta = int(120000000 / width); GDALRasterIOExtraArg psExtraArg; INIT_RASTERIO_EXTRA_ARG(psExtraArg); psExtraArg.eResampleAlg = eResample; do { rows_count = start_row + row_delta < height ? row_delta : height - start_row; cols_count = width; if (gdal_datatype == GDALDataType::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); demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg); delete[] temp; } else if (gdal_datatype == GDALDataType::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); demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg); delete[] temp; } else if (gdal_datatype == GDALDataType::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); demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg); delete[] temp; } start_row = start_row + rows_count; } while (start_row < height); GDALClose((GDALDatasetH)pDSrc); return 0; } int saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path) { int rows = data.rows(); int cols = data.cols(); Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3); gdalImage image_tiff = CreategdalImage(out_tiff_path, rows, cols, 2, gt, "", false, true);// 注意这里保留仿真结果 // 保存二进制文件 Eigen::MatrixXd real_img = data.array().real(); Eigen::MatrixXd imag_img = data.array().imag(); image_tiff.saveImage(real_img, 0, 0, 1); image_tiff.saveImage(imag_img, 0, 0, 2); return -1; } gdalImageComplex::gdalImageComplex(const QString& raster_path) { omp_lock_t lock; omp_init_lock(&lock); omp_set_lock(&lock); this->img_path = raster_path; GDALAllRegister(); GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 this->width = rasterDataset->GetRasterXSize(); this->height = rasterDataset->GetRasterYSize(); this->band_num = rasterDataset->GetRasterCount(); double* gt = new double[6]; rasterDataset->GetGeoTransform(gt); this->gt = Eigen::MatrixXd(2, 3); this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]; double a = this->gt(0, 0); double b = this->gt(0, 1); double c = this->gt(0, 2); double d = this->gt(1, 0); double e = this->gt(1, 1); double f = this->gt(1, 2); this->projection = rasterDataset->GetProjectionRef(); // 斤拷投影 GDALFlushCache((GDALDatasetH)rasterDataset); GDALClose((GDALDatasetH)rasterDataset); rasterDataset = NULL;// 指矫匡拷 this->InitInv_gt(); delete[] gt; ////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH omp_unset_lock(&lock); //锟酵放伙拷斤拷 omp_destroy_lock(&lock); //劫伙拷斤拷 } gdalImageComplex::~gdalImageComplex() { } void gdalImageComplex::setData(Eigen::MatrixXcd data) { this->data = data; } void gdalImageComplex::saveImage(Eigen::MatrixXcd 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; qDebug() << tip; throw exception(tip.toUtf8().constData()); } GDALAllRegister(); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); 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_CFloat64, 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; } int datarows = data.rows(); int datacols = data.cols(); double* databuffer = new double[data.size() * 2]; for (int i = 0; i < data.rows(); i++) { for (int j = 0; j < data.cols(); j++) { databuffer[i * data.cols() * 2 + j * 2] = data(i, j).real(); databuffer[i * data.cols() * 2 + j * 2 + 1] = data(i, j).imag(); } } //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_CFloat64, 0, 0); GDALFlushCache(poDstDS); delete databuffer; GDALClose((GDALDatasetH)poDstDS); //GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH omp_unset_lock(&lock); //锟酵放伙拷斤拷 omp_destroy_lock(&lock); //劫伙拷斤拷 } Eigen::MatrixXcd gdalImageComplex::getDataComplex(int start_row, int start_col, int rows_count, int cols_count, int band_ids) { GDALDataset* poDataset; GDALAllRegister(); // 打开TIFF文件 poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly); if (poDataset == nullptr) { QMessageBox::warning(nullptr, u8"错误", u8"无法打开文件:" + this->img_path); qDebug() << u8"无法打开文件:" + this->img_path; } // 获取数据集的第一个波段 GDALRasterBand* poBand; poBand = poDataset->GetRasterBand(1); // 读取波段信息,假设是复数类型 int nXSize = poBand->GetXSize(); int nYSize = poBand->GetYSize(); double* databuffer = new double[nXSize * nYSize * 2]; poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count, rows_count, GDT_CFloat64, 0, 0); GDALClose((GDALDatasetH)poDataset); Eigen::MatrixXcd rasterData(nYSize, nXSize); // 使用Eigen的MatrixXcd for (size_t i = 0; i < nYSize; i++) { for (size_t j = 0; j < nXSize; j++) { rasterData(i, j) = std::complex(databuffer[i * nXSize * 2 + j * 2], databuffer[i * nXSize * 2 + j * 2 + 1]); } } delete databuffer; return rasterData; } void gdalImageComplex::saveImage() { this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids); }