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

2427 lines
83 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

#include "stdafx.h"
#include "ImageOperatorBase.h"
#include "BaseTool.h"
#include "GeoOperator.h"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <omp.h>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <gdal.h>
#include <gdal_utils.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <proj.h>
#include <string.h>
#include <memory.h>
#include <memory>
#include <iostream>
#include "FileOperator.h"
#include <opencv2/opencv.hpp>
#include <QMessageBox>
#include <QDir>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <QProgressDialog>
/**
* 输入数据是ENVI格式数据
*/
std::shared_ptr<GDALDataset> OpenDataset(const QString& in_path, GDALAccess rwmode)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* dataset_ptr = (GDALDataset*)(GDALOpen(in_path.toUtf8().constData(), rwmode));
std::shared_ptr<GDALDataset> 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<GDALDataset> 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<GDALDataset> 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();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
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();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
std::shared_ptr<GDALDataset> 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<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
ReadComplexMatrixData(int start_line, int width, int line_num,
std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype)
{
int band_num = rasterDataset->GetRasterCount();
if(gdal_datatype == 0) {
return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(0, 0);
} else if(gdal_datatype < 8) {
if(band_num != 2) {
return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(0, 0);
}
} else if(gdal_datatype < 12) {
if(band_num != 1) {
return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(0, 0);
}
} else {
}
bool _flag = false;
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> data_mat(
line_num * width, 2); // 必须强制行优先
if(gdal_datatype == GDT_Byte) {
Eigen::MatrixX<char> real_mat(line_num * width, 1);
Eigen::MatrixX<char> 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<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_UInt16) {
Eigen::MatrixX<unsigned short> real_mat(line_num * width, 1);
Eigen::MatrixX<unsigned short> 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<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_Int16) {
Eigen::MatrixX<short> real_mat(line_num * width, 1);
Eigen::MatrixX<short> 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<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_UInt32) {
Eigen::MatrixX<unsigned int> real_mat(line_num * width, 1);
Eigen::MatrixX<unsigned int> 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<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_Int32) {
Eigen::MatrixX<int> real_mat(line_num * width, 1);
Eigen::MatrixX<int> 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<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
}
// else if (gdal_datatype == GDT_UInt64) {
// Eigen::MatrixX<unsigned long> real_mat(line_num * width, 1);
// Eigen::MatrixX<unsigned long> 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<double>()).array(); data_mat.col(1) =
//(imag_mat.array().cast<double>()).array(); _flag = true;
// }
// else if (gdal_datatype == GDT_Int64) {
// Eigen::MatrixX<long> real_mat(line_num * width, 1);
// Eigen::MatrixX<long> 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<double>()).array(); data_mat.col(1) =
//(imag_mat.array().cast<double>()).array(); _flag = true;
// }
else if(gdal_datatype == GDT_Float32) {
Eigen::MatrixX<float> real_mat(line_num * width, 1);
Eigen::MatrixX<float> 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<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_Float64) {
Eigen::MatrixX<double> real_mat(line_num * width, 1);
Eigen::MatrixX<double> 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<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_CInt16) {
Eigen::MatrixX<std::complex<short>> 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<double>()).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_CInt32) {
Eigen::MatrixX<std::complex<int>> 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<double>()).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_CFloat32) {
Eigen::MatrixX<std::complex<float>> 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<double>()).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_CFloat64) {
Eigen::MatrixX<std::complex<double>> 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<double>()).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
_flag = true;
} else {
}
// 保存数据
if(_flag) {
return data_mat;
} else {
return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
0, 0); // 必须强制行优先;
}
}
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
ReadMatrixDoubleData(int start_line, int width, int line_num,
std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype,
int band_idx)
{
// 构建矩阵块，使用eigen 进行矩阵计算，加速计算
bool _flag = false;
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> data_mat(
line_num * width, 1); // 必须强制行优先
if(gdal_datatype == GDT_Byte) {
Eigen::MatrixX<char> 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<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_UInt16) {
Eigen::MatrixX<unsigned short> 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<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_Int16) {
Eigen::MatrixX<short> 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<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_UInt32) {
Eigen::MatrixX<unsigned int> 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<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_Int32) {
Eigen::MatrixX<int> 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<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
}
// else if (gdal_datatype == GDT_UInt64) {
// Eigen::MatrixX<unsigned long> 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<double>()).array().pow(2)).log10() * 10.0; _flag = true;
// }
// else if (gdal_datatype == GDT_Int64) {
// Eigen::MatrixX<long> 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<double>()).array().pow(2)).log10() * 10.0; _flag = true;
// }
else if(gdal_datatype == GDT_Float32) {
Eigen::MatrixX<float> 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<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_Float64) {
Eigen::MatrixX<double> 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<double>()).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<GDALDataset> 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();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
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;
}
/// <summary>
/// 斤拷图饺∮帮拷锟?1锟?7
/// </summary>
/// <param name="dem_path"></param>
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();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
// 锟斤拷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();
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::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::MatrixXi gdalImage::getDatai(int start_row, int start_col, int rows_count, int cols_count, int band_ids)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
this->img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids);
rows_count = start_row + rows_count <= this->height ? rows_count : this->height - start_row;
cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
Eigen::MatrixXi datamatrix(rows_count, cols_count);
if (gdal_datatype == GDT_Byte) {
char* temp = new char[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0);
for (int i = 0; i < rows_count; i++) {
for (int j = 0; j < cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_UInt16) {
unsigned short* temp = new unsigned short[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0);
for (int i = 0; i < rows_count; i++) {
for (int j = 0; j < cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_Int16) {
short* temp = new short[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0);
for (int i = 0; i < rows_count; i++) {
for (int j = 0; j < cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_UInt32) {
unsigned int* temp = new unsigned int[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0);
for (int i = 0; i < rows_count; i++) {
for (int j = 0; j < cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_Int32) {
int* temp = new int[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0);
for (int i = 0; i < rows_count; i++) {
for (int j = 0; j < cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
//else if (gdal_datatype == GDT_UInt64) {
// unsigned long* temp = new unsigned long[rows_count * cols_count];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
// rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) {
// for (int j = 0; j <
// cols_count; j++) {
// datamatrix(i, j) = temp[i * cols_count + j];
// }
// }
// delete[] temp;
//}
//else if (gdal_datatype == GDT_Int64) {
// long* temp = new long[rows_count * cols_count];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
// for (int i = 0; i < rows_count; i++) {
// for (int j = 0; j < cols_count; j++) {
// datamatrix(i, j) = temp[i * cols_count + j];
// }
// }
// delete[] temp;
//}
else if (gdal_datatype == GDT_Float32) {
float* temp = new float[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0);
for (int i = 0; i < rows_count; i++) {
for (int j = 0; j < cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_Float64) {
double* temp = new double[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0);
for (int i = 0; i < rows_count; i++) {
for (int j = 0; j < cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
else {
}
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return datamatrix;
}
Eigen::MatrixXd gdalImage::getGeoTranslation()
{
return this->gt;
}
GDALDataType gdalImage::getDataType()
{
GDALDataset* rasterDataset =
(GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
return gdal_datatype;
}
/// <summary>
///
/// </summary>
/// <param name="data"></param>
/// <param name="start_row"></param>
/// <param name="start_col"></param>
/// <param name="band_ids"></param>
void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0,
int band_ids = 1)
{
GDALDataType datetype = this->getDataType();
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, datetype, NULL); // 斤拷锟斤拷
if (nullptr == poDstDS) {
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());
return;
}
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, datetype, 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(Eigen::MatrixXi data, int start_row, int start_col, int band_ids)
{
GDALDataType datetype=this->getDataType();
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, datetype, 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();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?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);
}
void gdalImage::setNoDataValuei(int nodatavalue, int band_ids)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?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);
}
double gdalImage::getNoDataValue(int band_ids)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
// GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
double v= poDstDS->GetRasterBand(band_ids)->GetNoDataValue( );
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
return v;
}
int gdalImage::getNoDataValuei(int band_ids)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
// GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
int v= poDstDS->GetRasterBand(band_ids)->GetNoDataValue( );
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
return v;
}
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;
}
void gdalImage::getLandPoint(double row, double col, double ati, Landpoint& Lp)
{
Lp.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); // x
Lp.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); // y
Lp.ati = ati;
}
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::BandmaxValue(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::BandminValue(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();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注斤拷锟斤拷
// 斤拷锟斤拷
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 std::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();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?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->BandminValue(bandids);
double dfMax = this->BandmaxValue(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;
}
RasterExtend gdalImage::getExtend()
{
RasterExtend extend{ 0,0,0,0 };
double x1 = this->gt(0, 0);
double y1 = this->gt(1, 0);
double x2 = this->gt(0, 0) + (this->width - 1) * gt(0, 1) + (0) * gt(0, 2);
double y2 = this->gt(1, 0) + (this->width - 1) * gt(1, 1) + (0) * gt(1, 2);
double x3 = this->gt(0, 0) + (0) * gt(0, 1) + (this->height - 1) * gt(0, 2);
double y3 = this->gt(1, 0) + (0) * gt(1, 1) + (this->height - 1) * gt(1, 2);
double x4 = this->gt(0, 0) + (this->width - 1) * gt(0, 1) + (this->height - 1) * gt(0, 2);
double y4 = this->gt(1, 0) + (this->width - 1) * gt(1, 1) + (this->height - 1) * gt(1, 2);
extend.min_x = x1 < x2 ? x1 : x2;
extend.max_x = x1 < x2 ? x2 : x1;
extend.min_y = y1 < y2 ? y1 : y2;
extend.max_y = y1 < y2 ? y2 : y1;
extend.min_x = extend.min_x < x3 ? extend.min_x : x3;
extend.max_x = extend.max_x > x3 ? extend.max_x : x3;
extend.min_y = extend.min_y < y3 ? extend.min_y : y3;
extend.max_y = extend.max_y > y3 ? extend.max_y : y3;
extend.min_x = extend.min_x < x4 ? extend.min_x : x4;
extend.max_x = extend.max_x > x4 ? extend.max_x : x4;
extend.min_y = extend.min_y < y4 ? extend.min_y : y4;
extend.max_y = extend.max_y > y4 ? extend.max_y : y4;
return extend;
}
gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num,
Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite)
{
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 = 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;
}
gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long epsgCode, GDALDataType eType, bool need_gt, bool overwrite, bool isENVI)
{
if (exists_test(img_path.toUtf8().constData())) {
if (overwrite) {
gdalImage result_img(img_path);
return result_img;
}
else {
throw "file has exist!!!";
exit(1);
}
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
GDALDriver* poDriver = isENVI? GetGDALDriverManager()->GetDriverByName("ENVI"): GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, eType, NULL); // 锟斤拷锟斤拷锟斤拷
if (need_gt) {
OGRSpatialReference oSRS;
if (oSRS.importFromEPSG(epsgCode) != OGRERR_NONE) {
std::cerr << "Failed to import EPSG code " << epsgCode << std::endl;
throw "Failed to import EPSG code ";
exit(1);
}
char* pszWKT = NULL;
oSRS.exportToWkt(&pszWKT);
std::cout << "WKT of EPSG:"<< epsgCode <<" :\n" << pszWKT << std::endl;
poDstDS->SetProjection(pszWKT);
double gt_ptr[6] = { 0 };
for (int i = 0; i < gt.rows(); i++) {
for (int j = 0; j < gt.cols(); j++) {
gt_ptr[i * 3 + j] = gt(i, j);
}
}
poDstDS->SetGeoTransform(gt_ptr);
}
for (int i = 1; i <= band_num; i++) {
poDstDS->GetRasterBand(i)->SetNoDataValue(-9999);
}
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
gdalImage result_img(img_path);
return result_img;
}
gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width,
int band_num, Eigen::MatrixXd gt, QString projection,
bool need_gt, bool overwrite)
{
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();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
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(0); // 回波部分
//}
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
gdalImageComplex result_img(img_path);
return result_img;
}
gdalImageComplex CreateEchoComplex(const QString& img_path, int height, int width, int band_num)
{
// 创建图像
Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
//Xgeo = GeoTransform[0] + Xpixel * GeoTransform[1] + Ypixel * GeoTransform[2]
//Ygeo = GeoTransform[3] + Xpixel * GeoTransform[4] + Ypixel * GeoTransform[5]
// X
gt(0, 0) = 0; gt(0, 2) = 1; gt(0, 2) = 0;
gt(1, 0) = 0; gt(1, 1) = 0; gt(1, 2) = 1;
// Y
QString projection = "";
gdalImageComplex echodata = CreategdalImageComplex(img_path, height, width, 1, gt, projection, false, true);
return echodata;
}
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width,
int new_height, GDALResampleAlg eResample)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
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<char*>(pDSrc->GetProjectionRef());
// 锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
if(strlen(pszSrcWKT) <= 0) {
OGRSpatialReference oSRS;
oSRS.importFromEPSG(4326);
// oSRS.SetUTM(50, true); //锟斤拷锟斤拷锟斤拷 锟斤拷锟斤拷120锟斤拷
// oSRS.SetWellKnownGeogCS("WGS84");
oSRS.exportToWkt(&pszSrcWKT);
}
qDebug() << "GDALCreateGenImgProjTransformer " << Qt::endl;
void* hTransformArg;
hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT,
FALSE, 0.0, 1);
qDebug() << "no proj " << Qt::endl;
//(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7)
if(hTransformArg == NULL) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -3;
}
qDebug() << "has proj " << Qt::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" << Qt::endl;
psWo->pfnTransformer = GDALGenImgProjTransform;
psWo->pTransformerArg = GDALCreateGenImgProjTransformer(
(GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);
;
qDebug() << "GDALCreateGenImgProjTransformer has created" << Qt::endl;
psWo->nBandCount = nBandCount;
psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));
for(int i = 0; i < nBandCount; i++) {
psWo->panSrcBands[i] = i + 1;
psWo->panDstBands[i] = i + 1;
}
GDALWarpOperation oWo;
if(oWo.Initialize(psWo) != CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
return -3;
}
qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << Qt::endl;
oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
GDALFlushCache(pDDst);
qDebug() << "ChunkAndWarpImage over" << Qt::endl;
// GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
// GDALDestroyWarpOptions(psWo);
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
////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", "YES");
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 alignRaster(QString inputPath, QString referencePath, QString outputPath, GDALResampleAlg eResample)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* pDSrc = (GDALDataset*)GDALOpen(inputPath.toUtf8().constData(), GA_ReadOnly);
if (pDSrc == NULL) {
return -1;
}
GDALDataset* tDSrc = (GDALDataset*)GDALOpen(referencePath.toUtf8().constData(), GA_ReadOnly);
if (tDSrc == NULL) {
return -1;
}
long new_width = tDSrc->GetRasterXSize();
long new_height = tDSrc->GetRasterYSize();
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<char*>(pDSrc->GetProjectionRef());
if (strlen(pszSrcWKT) <= 0) {
OGRSpatialReference oSRS;
oSRS.importFromEPSG(4326);
oSRS.exportToWkt(&pszSrcWKT);
}
qDebug() << "GDALCreateGenImgProjTransformer " << Qt::endl;
void* hTransformArg;
hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT,FALSE, 0.0, 1);
qDebug() << "no proj " << Qt::endl;
if (hTransformArg == NULL) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -3;
}
qDebug() << "has proj " << Qt::endl;
std::shared_ptr<double> dGeoTrans(new double[6], delArrPtr);
int nNewWidth = 0,
nNewHeight = 0;
if (GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg,dGeoTrans.get(), &nNewWidth, &nNewHeight) != CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -3;
}
GDALDataset* pDDst = pDriver->Create(outputPath.toUtf8().constData(), new_width, new_height, pDSrc->GetRasterCount(), dataType, NULL);
if (pDDst == NULL) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -2;
}
std::shared_ptr<double> gt(new double[6], delArrPtr);
tDSrc->GetGeoTransform(gt.get());
pDDst->SetProjection(pszSrcWKT);
pDDst->SetGeoTransform(gt.get());
GDALWarpOptions* psWo = GDALCreateWarpOptions();
psWo->eWorkingDataType = dataType;
psWo->eResampleAlg = eResample;
psWo->hSrcDS = (GDALDatasetH)pDSrc;
psWo->hDstDS = (GDALDatasetH)pDDst;
psWo->pfnTransformer = GDALGenImgProjTransform;
psWo->pTransformerArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);
qDebug() << "GDALCreateGenImgProjTransformer has created" << Qt::endl;
psWo->nBandCount = nBandCount;
psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));
for (int i = 0; i < nBandCount; i++) {
psWo->panSrcBands[i] = i + 1;
psWo->panDstBands[i] = i + 1;
}
GDALWarpOperation oWo;
if (oWo.Initialize(psWo) != CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)tDSrc);
return -3;
}
qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << Qt::endl;
oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
GDALFlushCache(pDDst);
qDebug() << "ChunkAndWarpImage over" << Qt::endl;
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)tDSrc);
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;
}
ErrorCode MergeRasterProcess(QVector<QString> filepaths, QString outfileptah, QString mainString, MERGEMODE mergecode, bool isENVI, ShowProessAbstract* dia )
{
// 参数检查
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_Int32, true, true, isENVI);
// 初始化
long resultline = Memory1MB * 500 / 8 / resultImage.width;
resultline = resultline < 10000 ? resultline : 10000; // 最多100行
resultline = resultline > 0 ? resultline : 2;
long bandnum = resultImage.band_num + 1;
long starti = 0;
long rasterCount = imgdslist.count();
QProgressDialog progressDialog(u8"初始化影像", u8"终止", 0, resultImage.height);
progressDialog.setWindowTitle(u8"初始化影像");
progressDialog.setWindowModality(Qt::WindowModal);
progressDialog.setAutoClose(true);
progressDialog.setValue(0);
progressDialog.setMaximum(resultImage.height);
progressDialog.setMinimum(0);
progressDialog.show();
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::MatrixXi maskdata = maskImage.getDatai(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);
}
if (nullptr != dia) {
dia->showProcess(starti * 1.0 / resultImage.height, u8"初始化影像数据");
}
progressDialog.setValue(starti+ blocklines);
}
progressDialog.close();
switch (mergecode)
{
case MERGE_GEOCODING:
return MergeRasterInGeoCoding(imgdslist, resultImage, maskImage, dia);
default:
break;
}
return ErrorCode::SUCCESS;
}
ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resultimg, gdalImage maskimg, ShowProessAbstract* dia)
{
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 * 300 / 8 / resultimg.width;
resultline = resultline < 1000 ? resultline : 1000; // 最多100行
long bandnum = resultimg.band_num+1;
long starti = 0;
long rasterCount = imgdslist.count();
long processNumber = 0;
QProgressDialog progressDialog(u8"合并影像", u8"终止", 0, resultimg.height);
progressDialog.setWindowTitle(u8"合并影像");
progressDialog.setWindowModality(Qt::WindowModal);
progressDialog.setAutoClose(true);
progressDialog.setValue(0);
progressDialog.setMaximum(resultimg.height);
progressDialog.setMinimum(0);
progressDialog.show();
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 = resultimg.getData(starti, 0, blocklines, resultimg.width, 1);
ridlist = ridlist.array() * 0;
Eigen::MatrixXd cidlist = ridlist.array() * 0;
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);
double nodata = imgdslist[ir].getNoDataValue(b);
for (long ii = 0; ii < data.rows(); ii++) {
for (long jj = 0; jj < data.cols(); jj++) {
if (std::abs(data(ii, jj) - nodata) < 1e-6) {
data(ii, jj) = 0;
}
}
}
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);
lastr = std::floor(ridtemp);
nextr = std::ceil(ridtemp);
lastc = std::floor(cidtemp);
lastc = std::ceil(cidtemp);
lastr = lastr - minRid;
nextr = nextr - minRid;
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";
if (nullptr != dia) {
dia->showProcess(processNumber * 1.0 / resultimg.height, u8"合并图像");
}
if (progressDialog.maximum() <= processNumber) {
processNumber = progressDialog.maximum() - 1;
}
progressDialog.setValue(processNumber);
omp_unset_lock(&lock);
}
omp_destroy_lock(&lock);
std::cout << std::endl;
progressDialog.setWindowTitle(u8"影像掩膜");
progressDialog.setLabelText(u8"影像掩膜");
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) {
data(i, j) = -9999;
continue;
}
data(i, j) = data(i, j) / maskdata(i, j);
}
}
resultimg.saveImage(data, starti, 0, b);
maskimg.saveImage(maskdata, starti, 0, b);
}
if (nullptr != dia) {
dia->showProcess((starti + blocklines) * 1.0 / resultimg.height, u8"影像掩膜");
}
progressDialog.setValue(starti + blocklines);
}
resultimg.setNoDataValue(-9999);
progressDialog.close();
return ErrorCode::SUCCESS;
}
gdalImageComplex::gdalImageComplex(const QString& raster_path)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
this->img_path = raster_path;
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
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 std::exception(tip.toUtf8().constData());
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
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();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
// 打开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<double>(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);
}
void gdalImageComplex::savePreViewImage()
{
qDebug()<<"void gdalImageComplex::savePreViewImage()";
Eigen::MatrixXd data_abs = Eigen::MatrixXd::Zero(this->height, this->width);
data_abs = (this->data.array().real().pow(2) + this->data.array().imag().pow(2))
.array()
.log10()*10.0; // 计算振幅
double min_abs = data_abs.minCoeff(); // 最大值
double max_abs = data_abs.maxCoeff(); // 最小值
double delta = (max_abs - min_abs) / 1000; // 1000分位档
Eigen::MatrixX<size_t> data_idx =
((data_abs.array() - min_abs).array() / delta).array().floor().cast<size_t>();
// 初始化
double hist[1001];
for(size_t i = 0; i < 1001; i++) {
hist[i] = 0; // 初始化
}
for(size_t i = 0; i < this->height; i++) {
for(size_t j = 0; j < this->width; j++) {
hist[data_idx(i, j)]++;
}
}
// 统计
size_t count = this->height * this->width;
double precent = 0;
size_t curCount = 0;
double pre2 = 0;
bool findprec_2 = true;
double pre98 = 0;
bool findprec_98 = true;
for(size_t i = 0; i < 1001; i++) {
precent = precent + hist[i];
if(findprec_2 & precent / count > 0.02) {
pre2 = i * delta + min_abs;
findprec_2 = false;
}
if(findprec_98 & precent / count > 0.98) {
pre98 = (i-1) * delta + min_abs;
findprec_98 = false;
}
}
// 拉伸
delta = (pre98-pre2)/200;
data_idx=
((data_abs.array() - pre2).array() / delta).array().floor().cast<size_t>();
for(size_t i = 0; i < this->height; i++) {
for(size_t j = 0; j < this->width; j++) {
if(data_idx(i,j)<0){
data_idx(i,j)=0;
}
else if(data_idx(i,j)>255){
data_idx(i,j)=255;
}else{
}
}
}
// 赋值
QString filePath = this->img_path;
QFile file(filePath);
QFileInfo fileInfo(file);
QString directory = fileInfo.absolutePath();
qDebug() << "文件所在目录：" << directory;
QString baseName = fileInfo.completeBaseName();
qDebug() << "无后缀文件名：" << baseName;
// 创建文件路径
QString previewImagePath = JoinPath(directory, baseName+"_preview.png");
cv::Mat img(this->height, this->width, CV_8U ,cv::Scalar(0));
for(size_t i = 0; i < this->height; i++) {
for(size_t j = 0; j < this->width; j++) {
img.at<uchar>(i,j)= (uchar)(data_idx(i,j));
}
}
//std::string outimgpath=previewImagePath.toUtf8().data();
cv::imwrite(previewImagePath.toUtf8().data(), img);
}
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;
}
std::cout << pszProjection << std::endl;
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;
}
}
void ShowProessAbstract::showProcess(double precent, QString tip)
{
}
void ShowProessAbstract::showToolInfo(QString tip)
{
}
Reference in New Issue View Git Blame Copy Permalink