LAMPMeasureDataProcessTeam
/
BaseCommonLibrary
forked from KJ135/BaseCommonLibrary
1
0
Fork 0
Code Pull Requests Activity
BaseCommonLibrary/BaseTool/ImageOperatorBase.cpp

1794 lines
62 KiB
C++
Raw Permalink 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>
#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>
#include <gdal_priv.h>
#include <ogr_spatialref.h> // OGRSpatialReference 用于空间参考转换
#include <gdal_alg.h> // 用于 GDALWarp 操作
/**
* 输入数据是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<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(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;
}
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);
qDebug() << "create init GDALDataset " ;
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();
CPLSetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); // 使用所有可用的CPU核心
CPLSetConfigOption("GDAL_CACHEMAX", "16000"); // 设置缓存大小为500MB
// psWo->papszWarpOptions = CSLDuplicate(NULL);
psWo->eWorkingDataType = dataType;
psWo->eResampleAlg = eResample;
psWo->hSrcDS = (GDALDatasetH)pDSrc;
psWo->hDstDS = (GDALDatasetH)pDDst;
qDebug() << "GDALCreateGenImgProjTransformer" << 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;
}
void clipRaster(QString inRasterPath, QString outRasterPath, long minRow, long maxRow, long minCol, long maxCol)
{
long rownum = maxRow - minRow + 1;
long colnum = maxCol - minCol + 1;
gdalImage inimg(inRasterPath);
Eigen::MatrixXd gt = inimg.gt;
Landpoint lp = inimg.getLandPoint(minRow, minCol, 0);
gt(0, 0) = lp.lon;
gt(1, 0) = lp.lat;
gdalImage outimg= CreategdalImageDouble(outRasterPath, rownum, colnum, inimg.band_num, gt, inimg.projection, true, true, true);
for (long bi = 1; bi < inimg.band_num + 1; bi++) {
Eigen::MatrixXd brasterData = inimg.getData(minRow, minCol, rownum, colnum, bi);
outimg.saveImage(brasterData, 0, 0, bi);
qDebug() << "writer raster band : " << bi;
}
qDebug() << "writer raster overring";
}
long getProjectEPSGCodeByLon_Lat(double long, 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驱动
// qDebug()<<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;
}
qDebug() << pszProjection ;
const char* epscodestr = oSRS.GetAuthorityCode(nullptr);
if (NULL == epscodestr || nullptr == epscodestr) {
qDebug() << "EPSG code string could not be determined from the spatial reference.";
GDALClose(poDataset);
return -1;
}
long epsgCode = atoi(epscodestr); // 获取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");
// qDebug()<<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;
}
}
else {
return CoordinateSystemType::UNKNOW;
}
}
void resampleRaster(const char* inputRaster, const char* outputRaster, double targetPixelSizeX, double targetPixelSizeY) {
// 初始化GDAL
GDALAllRegister();
// 打开输入栅格文件
GDALDataset* poDataset = (GDALDataset*)GDALOpen(inputRaster, GA_ReadOnly);
if (poDataset == nullptr) {
std::cerr << "Failed to open raster file." << std::endl;
return;
}
// 获取原始栅格的空间参考
double adfGeoTransform[6];
if (poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
std::cerr << "Failed to get geotransform." << std::endl;
GDALClose(poDataset);
return;
}
// 获取原始栅格的尺寸
int nXSize = poDataset->GetRasterXSize();
int nYSize = poDataset->GetRasterYSize();
// 计算目标栅格的尺寸
double targetXSize = (adfGeoTransform[1] * nXSize) / targetPixelSizeX;
double targetYSize = (adfGeoTransform[5] * nYSize) / targetPixelSizeY;
// 创建目标数据集（输出栅格）
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (poDriver == nullptr) {
std::cerr << "Failed to get GTiff driver." << std::endl;
GDALClose(poDataset);
return;
}
// 创建输出数据集
GDALDataset* poOutDataset = poDriver->Create(outputRaster, (int)targetXSize, (int)targetYSize, poDataset->GetRasterCount(), GDT_Float32, nullptr);
if (poOutDataset == nullptr) {
std::cerr << "Failed to create output raster." << std::endl;
GDALClose(poDataset);
return;
}
// 设置输出数据集的地理变换（坐标系）
double targetGeoTransform[6] = {
adfGeoTransform[0], targetPixelSizeX, 0, adfGeoTransform[3], 0, -targetPixelSizeY
};
poOutDataset->SetGeoTransform(targetGeoTransform);
// 设置输出数据集的投影信息
poOutDataset->SetProjection(poDataset->GetProjectionRef());
// 进行重采样
for (int i = 0; i < poDataset->GetRasterCount(); i++) {
GDALRasterBand* poBand = poDataset->GetRasterBand(i + 1);
GDALRasterBand* poOutBand = poOutDataset->GetRasterBand(i + 1);
// 使用GDAL的重采样方法，选择一个适当的重采样方式
poOutBand->RasterIO(GF_Write, 0, 0, (int)targetXSize, (int)targetYSize,
nullptr, (int)targetXSize, (int)targetYSize,
poBand->GetRasterDataType(), 0, 0, nullptr);
}
// 关闭数据集
GDALClose(poDataset);
GDALClose(poOutDataset);
qDebug() << "Resampling completed." ;
}
void transformRaster(const char* inputFile, const char* outputFile, int sourceEPSG, int targetEPSG) {
// 初始化 GDAL 库
GDALAllRegister();
// 打开源栅格文件
GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile, GA_ReadOnly);
if (poSrcDS == nullptr) {
qDebug() << "Failed to open input file:" << inputFile;
return;
}
// 获取源栅格的基本信息
int nXSize = poSrcDS->GetRasterXSize();
int nYSize = poSrcDS->GetRasterYSize();
int nBands = poSrcDS->GetRasterCount();
GDALDataType eDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();
// 创建目标栅格文件
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (poDriver == nullptr) {
qDebug() << "GTiff driver not available.";
GDALClose(poSrcDS);
return;
}
GDALDataset* poDstDS = poDriver->Create(outputFile, nXSize, nYSize, nBands, eDT, nullptr);
if (poDstDS == nullptr) {
qDebug() << "Failed to create output file:" << outputFile;
GDALClose(poSrcDS);
return;
}
// 设置目标栅格的空间参考系统
OGRSpatialReference oSRS;
oSRS.importFromEPSG(targetEPSG);
char* pszWKT = nullptr;
oSRS.exportToWkt(&pszWKT);
poDstDS->SetProjection(pszWKT);
CPLFree(pszWKT);
// 复制元数据
poDstDS->SetMetadata(poSrcDS->GetMetadata());
// 复制每个波段的数据
for (int i = 1; i <= nBands; ++i) {
GDALRasterBand* poSrcBand = poSrcDS->GetRasterBand(i);
GDALRasterBand* poDstBand = poDstDS->GetRasterBand(i);
float* pafScanline = (float*)CPLMalloc(sizeof(float) * nXSize);
for (int j = 0; j < nYSize; ++j) {
poSrcBand->RasterIO(GF_Read, 0, j, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
poDstBand->RasterIO(GF_Write, 0, j, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
}
CPLFree(pafScanline);
}
// 关闭数据集
GDALClose(poSrcDS);
GDALClose(poDstDS);
qDebug() << "Raster transformation completed successfully.";
}
ErrorCode transformCoordinate(double x, double y, int sourceEPSG, int targetEPSG, Point2 &p) {
// 创建源坐标系（原始坐标系）
OGRSpatialReference sourceSRS;
sourceSRS.importFromEPSG(sourceEPSG); // 使用 EPSG 编码来指定坐标系
// 创建目标坐标系（目标坐标系）
OGRSpatialReference targetSRS;
targetSRS.importFromEPSG(targetEPSG); // WGS84 坐标系 EPSG:4326
// 创建坐标变换对象
OGRCoordinateTransformation* transform = OGRCreateCoordinateTransformation(&sourceSRS, &targetSRS);
if (transform == nullptr) {
qDebug() << "Failed to create coordinate transformation." ;
return ErrorCode::FAIL;
}
// 转换坐标
double transformedX = x;
double transformedY = y;
if (transform->Transform(1, &transformedX, &transformedY)) {
qDebug() << "Original Coordinates: (" << x << ", " << y << ")" ;
qDebug() << "Transformed Coordinates (EPSG:" << targetEPSG << "): (" << transformedX << ", " << transformedY << ")" ;
}
else {
qDebug() << "Coordinate transformation failed." ;
}
// 清理
delete transform;
p.x = transformedX;
p.y = transformedY;
return ErrorCode::SUCCESS;
}
void cropRasterByLatLon(const char* inputFile, const char* outputFile,double minLon, double maxLon, double minLat, double maxLat) {
// 初始化 GDAL 库
GDALAllRegister();
// 打开栅格数据集
GDALDataset* poDataset = (GDALDataset*)GDALOpen(inputFile, GA_ReadOnly);
if (poDataset == nullptr) {
std::cerr << "Failed to open input raster." << std::endl;
return;
}
// 获取栅格数据的地理参考信息
double adfGeoTransform[6];
if (poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
std::cerr << "Failed to get geotransform." << std::endl;
GDALClose(poDataset);
return;
}
// 获取输入影像的投影信息
const char* projection = poDataset->GetProjectionRef();
// 根据经纬度计算出裁剪区域对应的栅格像素坐标
int xMin = (int)((minLon - adfGeoTransform[0]) / adfGeoTransform[1]);
int xMax = (int)((maxLon - adfGeoTransform[0]) / adfGeoTransform[1]);
int yMin = (int)((maxLat - adfGeoTransform[3]) / adfGeoTransform[5]);
int yMax = (int)((minLat - adfGeoTransform[3]) / adfGeoTransform[5]);
// 创建裁剪区域的目标栅格数据集
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (poDriver == nullptr) {
std::cerr << "Failed to get GTiff driver." << std::endl;
GDALClose(poDataset);
return;
}
// 创建输出栅格数据集，指定尺寸
int width = xMax - xMin;
int height = yMax - yMin;
GDALDataset* poOutDataset = poDriver->Create(outputFile, width, height, poDataset->GetRasterCount(), GDT_Float32, nullptr);
if (poOutDataset == nullptr) {
std::cerr << "Failed to create output raster." << std::endl;
GDALClose(poDataset);
return;
}
// 设置输出栅格的投影信息和地理变换
poOutDataset->SetProjection(projection);
double newGeoTransform[6] = { adfGeoTransform[0] + xMin * adfGeoTransform[1], adfGeoTransform[1], 0.0, adfGeoTransform[3] + yMin * adfGeoTransform[5], 0.0, adfGeoTransform[5] };
poOutDataset->SetGeoTransform(newGeoTransform);
// 循环读取源数据并写入目标数据集
for (int i = 0; i < poDataset->GetRasterCount(); ++i) {
GDALRasterBand* poBand = poDataset->GetRasterBand(i + 1);
GDALRasterBand* poOutBand = poOutDataset->GetRasterBand(i + 1);
// 读取源数据
int* pData = new int[width * height];
poBand->RasterIO(GF_Read, xMin, yMin, width, height, pData, width, height, GDT_Int32, 0, 0);
// 写入目标数据
poOutBand->RasterIO(GF_Write, 0, 0, width, height, pData, width, height, GDT_Int32, 0, 0);
delete[] pData;
}
qDebug() << "Raster cropped and saved to: " << outputFile ;
// 清理
GDALClose(poDataset);
GDALClose(poOutDataset);
}
ErrorCode DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QString demsloperPath)
{
gdalImage demds(dempath);
gdalImage demxyz = CreategdalImage(demxyzpath, demds.height, demds.width, 3, demds.gt, demds.projection, true, true);// X,Y,Z
// 分块计算并转换为XYZ
Eigen::MatrixXd demArr = demds.getData(0, 0, demds.height, demds.width, 1);
Eigen::MatrixXd demR = demArr;
Landpoint LandP{ 0,0,0 };
Point3 GERpoint{ 0,0,0 };
double R = 0;
double dem_row = 0, dem_col = 0, dem_alt = 0;
long line_invert = 1000;
double rowidx = 0;
double colidx = 0;
for (int max_rows_ids = 0; max_rows_ids < demds.height; max_rows_ids = max_rows_ids + line_invert) {
Eigen::MatrixXd demdata = demds.getData(max_rows_ids, 0, line_invert, demds.width, 1);
Eigen::MatrixXd xyzdata_x = demdata.array() * 0;
Eigen::MatrixXd xyzdata_y = demdata.array() * 0;
Eigen::MatrixXd xyzdata_z = demdata.array() * 0;
int datarows = demdata.rows();
int datacols = demdata.cols();
for (int i = 0; i < datarows; i++) {
for (int j = 0; j < datacols; j++) {
rowidx = i + max_rows_ids;
colidx = j;
demds.getLandPoint(rowidx, colidx, demdata(i, j), LandP); // 获取地理坐标
LLA2XYZ(LandP, GERpoint); // 经纬度转换为地心坐标系
xyzdata_x(i, j) = GERpoint.x;
xyzdata_y(i, j) = GERpoint.y;
xyzdata_z(i, j) = GERpoint.z;
}
}
demxyz.saveImage(xyzdata_x, max_rows_ids, 0, 1);
demxyz.saveImage(xyzdata_y, max_rows_ids, 0, 2);
demxyz.saveImage(xyzdata_z, max_rows_ids, 0, 3);
}
// 计算坡向角
gdalImage demsloperxyz = CreategdalImage(demsloperPath, demds.height, demds.width, 4, demds.gt, demds.projection, true, true);// X,Y,Z,cosangle
line_invert = 1000;
long start_ids = 0;
long dem_rows = 0, dem_cols = 0;
for (start_ids = 1; start_ids < demds.height; start_ids = start_ids + line_invert) {
Eigen::MatrixXd demdata = demds.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
long startlineid = start_ids;
Eigen::MatrixXd demsloper_x = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
Eigen::MatrixXd demsloper_y = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 2);
Eigen::MatrixXd demsloper_z = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 3);
Eigen::MatrixXd demsloper_angle = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 4);
Landpoint p0, p1, p2, p3, p4, pslopeVector, pp;
Vector3D slopeVector;
dem_rows = demsloper_y.rows();
dem_cols = demsloper_y.cols();
double sloperAngle = 0;
Vector3D Zaxis = { 0,0,1 };
double rowidx = 0, colidx = 0;
for (long i = 1; i < dem_rows - 1; i++) {
for (long j = 1; j < dem_cols - 1; j++) {
rowidx = i + startlineid;
colidx = j;
demds.getLandPoint(rowidx, colidx, demdata(i, j), p0);
demds.getLandPoint(rowidx - 1, colidx, demdata(i - 1, j), p1);
demds.getLandPoint(rowidx, colidx - 1, demdata(i, j - 1), p2);
demds.getLandPoint(rowidx + 1, colidx, demdata(i + 1, j), p3);
demds.getLandPoint(rowidx, colidx + 1, demdata(i, j + 1), p4);
pslopeVector = getSlopeVector(p0, p1, p2, p3, p4); // 地面坡向矢量
slopeVector = { pslopeVector.lon,pslopeVector.lat,pslopeVector.ati };
pp = LLA2XYZ(p0);
Zaxis.x = pp.lon;
Zaxis.y = pp.lat;
Zaxis.z = pp.ati;
sloperAngle = getCosAngle(slopeVector, Zaxis); // 地面坡向角
demsloper_x(i, j) = slopeVector.x;
demsloper_y(i, j) = slopeVector.y;
demsloper_z(i, j) = slopeVector.z;
demsloper_angle(i, j) = sloperAngle;
}
}
demsloperxyz.saveImage(demsloper_x, start_ids - 1, 0, 1);
demsloperxyz.saveImage(demsloper_y, start_ids - 1, 0, 2);
demsloperxyz.saveImage(demsloper_z, start_ids - 1, 0, 3);
demsloperxyz.saveImage(demsloper_angle, start_ids - 1, 0, 4);
}
return ErrorCode::SUCCESS;
}
void ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode) {
// 注册所有GDAL驱动
GDALAllRegister();
// 打开输入栅格文件
GDALDataset* srcDataset = (GDALDataset*)GDALOpen(inRasterPath.toUtf8().constData(), GA_ReadOnly);
if (!srcDataset) {
// 错误处理：输出文件打开失败
// qDebug() << "无法打开输入文件：" << inRasterPath;
return;
}
// 创建目标坐标系
OGRSpatialReference targetSRS;
if (targetSRS.importFromEPSG(outepsgcode) != OGRERR_NONE) {
GDALClose(srcDataset);
// qDebug() << "无效的EPSG代码：" << outepsgcode;
return;
}
GDALDataType datetype = srcDataset->GetRasterBand(1)->GetRasterDataType();
// 获取目标坐标系的WKT表示
char* targetSRSWkt = nullptr;
targetSRS.exportToWkt(&targetSRSWkt);
bool flag = (datetype == GDT_Byte || datetype == GDT_Int8 || datetype == GDT_Int16 ||datetype == GDT_UInt16 || datetype == GDT_Int32 || datetype == GDT_UInt32 || datetype == GDT_Int64 || datetype == GDT_UInt64);
// 创建重投影后的虚拟数据集（Warped VRT）
GDALDataset* warpedVRT = flag? (GDALDataset*)GDALAutoCreateWarpedVRT(
srcDataset,
nullptr, // 输入坐标系（默认使用源数据）
targetSRSWkt, // 目标坐标系
GRA_NearestNeighbour, // 重采样方法：双线性插值
0.0, // 最大误差（0表示自动计算）
nullptr // 其他选项
) :(GDALDataset*)GDALAutoCreateWarpedVRT(
srcDataset,
nullptr, // 输入坐标系（默认使用源数据）
targetSRSWkt, // 目标坐标系
GRA_Bilinear, // 重采样方法：双线性插值
0.0, // 最大误差（0表示自动计算）
nullptr // 其他选项
);
CPLFree(targetSRSWkt); // 释放WKT内存
if (!warpedVRT) {
GDALClose(srcDataset);
qDebug() << u8"创建投影转换VRT失败";
return;
}
// 获取输出驱动（GeoTIFF格式）
QString filesuffer = getFileExtension(outRasterPath).toLower();
bool isTiff = filesuffer.contains("tif");
GDALDriver* driver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
if (!driver) {
GDALClose(warpedVRT);
GDALClose(srcDataset);
// qDebug() << "无法获取GeoTIFF驱动";
return;
}
// 创建输出栅格文件
GDALDataset* dstDataset = driver->CreateCopy(
outRasterPath.toUtf8().constData(), // 输出文件路径
warpedVRT, // 输入数据集（VRT）
false, // 是否严格复制
nullptr, // 创建选项
nullptr, // 进度回调
nullptr // 回调参数
);
if (!dstDataset) {
// qDebug() << "创建输出文件失败：" << outRasterPath;
GDALClose(warpedVRT);
GDALClose(srcDataset);
return;
}
// 释放资源
GDALClose(dstDataset);
GDALClose(warpedVRT);
GDALClose(srcDataset);
}
void ResampleByReferenceRasterB(QString pszSrcFile, QString RefrasterBPath, QString pszOutFile, GDALResampleAlg eResample) {
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile.toLocal8Bit().constData(), GA_ReadOnly);
if (pDSrc == NULL) {
qDebug() << u8"do not open In Raster file: " << pszSrcFile;
return ;
}
GDALDataset* pDRef = (GDALDataset*)GDALOpen(RefrasterBPath.toLocal8Bit().constData(), GA_ReadOnly);
if (pDRef == NULL) {
qDebug() << u8"do not open Ref Raster file: " << RefrasterBPath;
return;
}
QString filesuffer = getFileExtension(pszOutFile).toLower();
bool isTiff = filesuffer.contains("tif");
GDALDriver* pDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
if (pDriver == NULL) {
qDebug() << "not open driver";
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return ;
}
int width = pDSrc->GetRasterXSize();
int height = pDSrc->GetRasterYSize();
int nBandCount = pDSrc->GetRasterCount();
GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();
char* pszSrcWKT = NULL;
pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef());
// 锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
if (strlen(pszSrcWKT) <= 0) {
OGRSpatialReference oSRS;
oSRS.importFromEPSG(4326);
// oSRS.SetUTM(50, true); //锟斤拷锟斤拷锟斤拷 锟斤拷锟斤拷120锟斤拷
// oSRS.SetWellKnownGeogCS("WGS84");
oSRS.exportToWkt(&pszSrcWKT);
}
char* pdstSrcWKT = NULL;
pdstSrcWKT = const_cast<char*>(pDRef->GetProjectionRef());
// 锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
if (strlen(pdstSrcWKT) <= 0)
{
OGRSpatialReference oSRS;
oSRS.importFromEPSG(4326);
// oSRS.SetUTM(50, true); //锟斤拷锟斤拷锟斤拷 锟斤拷锟斤拷120锟斤拷
// oSRS.SetWellKnownGeogCS("WGS84");
oSRS.exportToWkt(&pdstSrcWKT);
}
int new_width = pDRef->GetRasterXSize();
int new_height = pDRef->GetRasterYSize();
double gt[6] ;
pDRef->GetGeoTransform(gt);
// GDALDestroyGenImgProjTransformer(hTransformArg);
qDebug() << "create init GDALDataset ";
GDALDataset* pDDst =
pDriver->Create(pszOutFile.toLocal8Bit().constData(), new_width, new_height, nBandCount, dataType, NULL);
if (pDDst == NULL) {
qDebug() << "not create init GDALDataset ";
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
return ;
}
pDDst->SetProjection(pdstSrcWKT);
pDDst->SetGeoTransform(gt);
qDebug() << "GDALCreateGenImgProjTransformer " << Qt::endl;
void* hTransformArg;
hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT,
FALSE, 0.0, 1);
qDebug() << "no proj ";
//(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7)
if (hTransformArg == NULL) {
qDebug() << "hTransformArg create failure";
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
return;
}
qDebug() << "has proj ";
double dGeoTrans[6] = { 0 };
int nNewWidth = 0, nNewHeight = 0;
if (GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg,
dGeoTrans, &nNewWidth, &nNewHeight)
!= CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
return ;
}
GDALWarpOptions* psWo = GDALCreateWarpOptions();
CPLSetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); // 使用所有可用的CPU核心
CPLSetConfigOption("GDAL_CACHEMAX", "16000"); // 设置缓存大小为500MB
// psWo->papszWarpOptions = CSLDuplicate(NULL);
psWo->eWorkingDataType = dataType;
psWo->eResampleAlg = eResample;
psWo->hSrcDS = (GDALDatasetH)pDSrc;
psWo->hDstDS = (GDALDatasetH)pDDst;
qDebug() << "GDALCreateGenImgProjTransformer" ;
psWo->pfnTransformer = GDALGenImgProjTransform;
psWo->pTransformerArg = GDALCreateGenImgProjTransformer(
(GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);
;
qDebug() << "GDALCreateGenImgProjTransformer has created" << Qt::endl;
psWo->nBandCount = nBandCount;
psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));
for (int i = 0; i < nBandCount; i++) {
psWo->panSrcBands[i] = i + 1;
psWo->panDstBands[i] = i + 1;
}
GDALWarpOperation oWo;
if (oWo.Initialize(psWo) != CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
return ;
}
qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << Qt::endl;
oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
GDALFlushCache(pDDst);
qDebug() << "ChunkAndWarpImage over" << Qt::endl;
// GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
// GDALDestroyWarpOptions(psWo);
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return ;
}
void ResampleByReferenceRasterB(QString InrasterAPath, QString RefrasterBPath, QString OutrasterCPath) {
// 注册所有GDAL驱动
GDALAllRegister();
// 打开参考栅格B
GDALDataset* refDS = (GDALDataset*)GDALOpen(RefrasterBPath.toUtf8().constData(), GA_ReadOnly);
if (!refDS) {
qDebug() << "无法打开参考栅格B：" << RefrasterBPath;
return;
}
// 获取参考栅格的地理变换、投影和尺寸
double geotransform[6];
if (refDS->GetGeoTransform(geotransform) != CE_None) {
qDebug() << "获取参考栅格的地理变换失败。";
GDALClose(refDS);
return;
}
const char* proj = refDS->GetProjectionRef();
int cols = refDS->GetRasterXSize();
int rows = refDS->GetRasterYSize();
GDALClose(refDS); // 获取信息后关闭参考栅格
// 打开输入栅格A
GDALDataset* srcDS = (GDALDataset*)GDALOpen(InrasterAPath.toUtf8().constData(), GA_ReadOnly);
if (!srcDS) {
qDebug() << "无法打开输入栅格A：" << InrasterAPath;
return;
}
// 获取输入栅格的波段数和数据类型
int nBands = srcDS->GetRasterCount();
if (nBands == 0) {
qDebug() << "输入栅格没有波段数据。";
GDALClose(srcDS);
return;
}
GDALDataType dataType = srcDS->GetRasterBand(1)->GetRasterDataType();
// 创建输出栅格C
GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (!driver) {
qDebug() << "无法获取GeoTIFF驱动。";
GDALClose(srcDS);
return;
}
GDALDataset* dstDS = driver->Create(
OutrasterCPath.toUtf8().constData(),
cols,
rows,
nBands,
dataType,
nullptr
);
if (!dstDS) {
qDebug() << "无法创建输出栅格：" << OutrasterCPath;
GDALClose(srcDS);
return;
}
// 设置输出栅格的地理变换和投影
dstDS->SetGeoTransform(geotransform);
dstDS->SetProjection(proj);
// 配置GDAL Warp选项
GDALWarpOptions* psWO = GDALCreateWarpOptions();
psWO->hSrcDS = srcDS;
psWO->hDstDS = dstDS;
psWO->nBandCount = nBands;
psWO->panSrcBands = (int*)CPLMalloc(nBands * sizeof(int));
psWO->panDstBands = (int*)CPLMalloc(nBands * sizeof(int));
for (int i = 0; i < nBands; ++i) {
psWO->panSrcBands[i] = i + 1;
psWO->panDstBands[i] = i + 1;
}
psWO->eResampleAlg = GRA_NearestNeighbour; // 使用最近邻重采样
// 初始化坐标转换器
psWO->pfnTransformer = GDALGenImgProjTransform;
psWO->pTransformerArg = GDALCreateGenImgProjTransformer(
srcDS, GDALGetProjectionRef(srcDS),
dstDS, GDALGetProjectionRef(dstDS),
FALSE, 0.0, 1
);
if (!psWO->pTransformerArg) {
qDebug() << "创建坐标转换器失败。";
GDALDestroyWarpOptions(psWO);
GDALClose(srcDS);
GDALClose(dstDS);
return;
}
// 执行Warp操作
GDALWarpOperation oWarp;
if (oWarp.Initialize(psWO) != CE_None) {
qDebug() << "初始化Warp操作失败。";
GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);
GDALDestroyWarpOptions(psWO);
GDALClose(srcDS);
GDALClose(dstDS);
return;
}
CPLErr eErr = oWarp.ChunkAndWarpImage(0, 0, cols, rows);
if (eErr != CE_None) {
qDebug() << "执行Warp操作失败。";
}
// 清理资源
GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);
GDALDestroyWarpOptions(psWO);
GDALClose(srcDS);
GDALClose(dstDS);
qDebug() << "重采样完成，结果已保存至：" << OutrasterCPath;
}
void CreateSARIntensityByLookTable(QString IntensityRasterPath,
QString LookTableRasterPath,
QString SARIntensityPath,
long min_rid, long max_rid,
long min_cid, long max_cid,
std::function<void(long, long)> processBarShow
)
{
gdalImage looktableds(LookTableRasterPath);
gdalImage geoIntensity(IntensityRasterPath);
gdalImage SARIntensity= CreategdalImageDouble(SARIntensityPath,max_rid-min_rid,max_cid-min_cid,1);
long blockYSize = Memory1GB / looktableds.width / 8 * 2;
Eigen::MatrixXd SARData = SARIntensity.getData(0, 0, SARIntensity.height, SARIntensity.width, 1);
SARData = SARData.array() * 0;
// 分块处理
for (int yOff = 0; yOff < looktableds.height; yOff += blockYSize)
{
processBarShow(yOff, looktableds.height);
qDebug() << "Process : [" << yOff * 100.0 / looktableds.height << " % ]";
Eigen::MatrixXd rowData = looktableds.getData(yOff, 0, blockYSize, looktableds.width, 1);
Eigen::MatrixXd colData = looktableds.getData(yOff, 0, blockYSize, looktableds.width, 2);
Eigen::MatrixXd geoData = geoIntensity.getData(yOff, 0, blockYSize, looktableds.width, 1);
for (long i = 0; i < rowData.rows(); i++) {
for (long j = 0; j < rowData.cols(); j++) {
long r =round( rowData(i,j))-min_rid;
long c = round(colData(i, j))-min_cid;
if (r >= 0 && r < SARIntensity.height && c >= 0 && c < SARIntensity.width) {
SARData(r, c) = SARData(r, c) + geoData(i, j);
}
}
}
}
SARIntensity.saveImage(SARData, 0, 0, 1);
qDebug() << "Process : [ 100 % ]";
processBarShow(1000,1000);
}
bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath)
{
Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
gdalImage img = CreategdalImageDouble(dataStrPath, data.rows(), data.cols(), 1, gt, "", false, true, true);
img.saveImage(data, 0, 0, 1);
return true;
}
View Git Blame Copy Permalink