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

2155 lines
71 KiB
C++
Raw Permalink Normal View History

初始化模块
2025-05-27 15:05:39 +00:00
#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 操作
#include "gdal_priv.h"
#include "cpl_conv.h"
#include <iostream>
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 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;
}
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 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;
}
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;
}
void ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode) {
// 注册所有GDAL驱动
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
// 打开输入栅格文件
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");
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 cropRasterByLatLon(const char* inputFile, const char* outputFile, double minLon, double maxLon, double minLat, double maxLat) {
// 初始化 GDAL 库
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
// 打开栅格数据集
GDALDataset* poDataset = (GDALDataset*)GDALOpen(inputFile, GA_ReadOnly);
if (poDataset == nullptr) {
qDebug() << "Failed to open input raster." ;
return;
}
// 获取栅格数据的地理参考信息
double adfGeoTransform[6];
if (poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
qDebug() << "Failed to get geotransform." ;
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) {
qDebug() << "Failed to get GTiff driver." ;
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) {
qDebug() << "Failed to create output raster." ;
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((GDALDatasetH)(GDALDatasetH)poDataset);
GDALClose((GDALDatasetH)(GDALDatasetH)poOutDataset);
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
}
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 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.";
}
void resampleRaster(const char* inputRaster, const char* outputRaster, double targetPixelSizeX, double targetPixelSizeY) {
// 初始化GDAL
GDALAllRegister();
// 打开输入栅格文件
GDALDataset* poDataset = (GDALDataset*)GDALOpen(inputRaster, GA_ReadOnly);
if (poDataset == nullptr) {
qDebug() << "Failed to open raster file." ;
return;
}
// 获取原始栅格的空间参考
double adfGeoTransform[6];
if (poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
qDebug() << "Failed to get geotransform." ;
GDALClose(poDataset);
return;
}
// 获取原始栅格的尺寸
int nXSize = poDataset->GetRasterXSize();
int nYSize = poDataset->GetRasterYSize();
// 计算目标栅格的尺寸
double targetXSize = ceil((adfGeoTransform[1] * nXSize) / targetPixelSizeX);
double targetYSize = ceil(std::abs((adfGeoTransform[5] * nYSize) / targetPixelSizeY));
// 创建目标数据集（输出栅格）
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (poDriver == nullptr) {
qDebug() << "Failed to get GTiff driver." ;
GDALClose(poDataset);
return;
}
// 设置输出数据集的地理变换（坐标系）
double targetGeoTransform[6] = {
adfGeoTransform[0], targetPixelSizeX, 0, adfGeoTransform[3], 0, -targetPixelSizeY
};
GDALClose((GDALDatasetH)(GDALDatasetH)poDataset);
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
ResampleGDAL(inputRaster, outputRaster, targetGeoTransform, targetXSize, targetYSize,GRA_Bilinear);
qDebug() << "Resampling completed.";
}
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;
}
}
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;
}
}
}
long getProjectEPSGCodeByLon_Lat(double lon, double lat, ProjectStripDelta stripState)
{
long EPSGCode = 0;
switch (stripState) {
case ProjectStripDelta::Strip_3: {
EPSGCode = getProjectEPSGCodeByLon_Lat_inStrip3(lon, lat);
break;
};
case ProjectStripDelta::Strip_6: {
EPSGCode = getProjectEPSGCodeByLon_Lat_inStrip6(lon, lat);
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;
}
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;
}
}
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", "4000"); // 设置缓存大小为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;
}
/** 2025.03.25 下面的函数存在 Eigen使用****************************/
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;
}
Eigen::MatrixXd getGeoTranslationArray(QString in_path)
{
return Eigen::MatrixXd();
}
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";
}
ErrorCode DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QString demsloperPath)
{
omp_lock_t lock;
omp_init_lock(&lock);
QString DEMPath = dempath;
QString XYZPath = demxyzpath;
QString SLOPERPath = demsloperPath;
gdalImage demds(DEMPath);
gdalImage demxyz = CreategdalImageDouble(XYZPath, demds.height, demds.width, 3, demds.gt, demds.projection, true, true, false);// X,Y,Z
int64_t linecount = 0;
qDebug() << u8"start dem (lon,lat,ati) -> dem [X,Y,Z]";
int64_t line_invert= Memory1MB / 8.0 / demds.width * 200;
#pragma omp parallel for
for (int64_t max_rows_ids = 0; max_rows_ids < demds.height; max_rows_ids = max_rows_ids + line_invert) {
long temp_line_invert = max_rows_ids + line_invert < demds.height ? line_invert : demds.height - max_rows_ids;
double R = 0;
double dem_row = 0, dem_col = 0, dem_alt = 0;
Eigen::MatrixXd demdata = demds.getData(max_rows_ids, 0, temp_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++) {
Landpoint LandP{ 0,0,0 };
Point3 GERpoint{ 0,0,0 };
double rowidx = 0;
double colidx = 0;
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;
}
}
omp_set_lock(&lock);
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);
linecount = linecount + temp_line_invert;
qDebug() << "dem -> XYZ [" << linecount * 100.0 / demds.height << "] %";
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
}
// 计算坡向角
qDebug() << u8"start dem (lon,lat,ati) -> dem Sloper[X,Y,Z]";
gdalImage demsloperxyz = CreategdalImageDouble(SLOPERPath, demds.height, demds.width, 4, demds.gt, demds.projection, true, true, false);// X,Y,Z,cosangle
linecount = 0;
line_invert = Memory1MB / 8.0 / demds.width * 200;;
#pragma omp parallel for
for (int64_t start_ids = 1; start_ids < demds.height; start_ids = start_ids + line_invert) {
long temp_line_invert = start_ids + line_invert < demds.height ? line_invert : demds.height - start_ids;
long dem_rows = 0, dem_cols = 0;
//long startlineid = start_ids;
Eigen::MatrixXd demdata = demds.getData(start_ids - 1, 0, temp_line_invert + 2, demxyz.width, 1);
Eigen::MatrixXd demsloper_x = demsloperxyz.getData(start_ids , 0, temp_line_invert , demxyz.width, 1);
Eigen::MatrixXd demsloper_y = demsloperxyz.getData(start_ids , 0, temp_line_invert , demxyz.width, 2);
Eigen::MatrixXd demsloper_z = demsloperxyz.getData(start_ids , 0, temp_line_invert , demxyz.width, 3);
Eigen::MatrixXd demsloper_angle = demsloperxyz.getData(start_ids , 0, temp_line_invert , demxyz.width, 4);
dem_rows = demsloper_x.rows();
dem_cols = demsloper_x.cols();
for (long i = 0; i < dem_rows ; i++) {
Landpoint p0, p1, p2, p3, p4, pslopeVector, pp;
Vector3D slopeVector;
double sloperAngle = 0;
Vector3D Zaxis = { 0,0,1 };
double rowidx = 0, colidx = 0;
for (long j = 1; j < dem_cols - 1; j++) {
rowidx = i + 1;
colidx = j;
demds.getLandPoint(rowidx, colidx, demdata(i + 1, j + 1), p0);
demds.getLandPoint(rowidx - 1, colidx, demdata(i - 1 + 1, j), p1);
demds.getLandPoint(rowidx, colidx - 1, demdata(i + 1, j - 1), p2);
demds.getLandPoint(rowidx + 1, colidx, demdata(i + 1 + 1, j), p3);
demds.getLandPoint(rowidx, colidx + 1, demdata(i + 1, 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;
}
}
omp_set_lock(&lock);
demsloperxyz.saveImage(demsloper_x, start_ids , 0, 1);
demsloperxyz.saveImage(demsloper_y, start_ids , 0, 2);
demsloperxyz.saveImage(demsloper_z, start_ids , 0, 3);
demsloperxyz.saveImage(demsloper_angle, start_ids , 0, 4);
linecount = linecount + temp_line_invert;
qDebug() << "dem -> Sloper [" << linecount * 100.0 / demds.height << "] %";
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
}
omp_destroy_lock(&lock); // 劫伙拷斤拷
return ErrorCode::SUCCESS;
}
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);
}
/**​
* @brief 将VRT文件转换为ENVI格式（生成.dat和.hdr文件）
* @param inputVrtPath 输入VRT文件路径
* @param outputEnviPath 输出ENVI数据文件路径（不含扩展名，自动生成.dat和.hdr）
* @param papszOptions GDAL创建选项（如压缩参数、存储顺序等，默认NULL）
* @return 成功返回true，失败返回false
*/
bool ConvertVrtToEnvi(QString vrtPath, QString outPath) {
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
const char* inputVrtPath = vrtPath.toUtf8().constData();
const char* outputEnviPath = outPath.toUtf8().constData();
// 注册所有GDAL驱动[7](@ref)
GDALAllRegister();
// 打开输入VRT文件[1,3](@ref)
GDALDataset* poVRTDataset = (GDALDataset*)GDALOpen(inputVrtPath, GA_ReadOnly);
if (!poVRTDataset) {
qDebug() << "错误：无法打开VRT文件 " << inputVrtPath ;
return false;
}
// 获取ENVI驱动[4,7](@ref)
GDALDriver* poENVIDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
if (!poENVIDriver) {
qDebug() << "错误：ENVI驱动未加载" ;
GDALClose(poVRTDataset);
return false;
}
// 构造输出文件名（强制添加.dat后缀）[4](@ref)
std::string outputPath = std::string(outputEnviPath) + ".dat";
// 执行格式转换（自动生成.hdr头文件）[4,7](@ref)
GDALDataset* poENVIDataset = poENVIDriver->CreateCopy(
outputPath.c_str(), // 输出路径
poVRTDataset, // 输入数据集
FALSE, // 非严格模式（允许格式调整）
nullptr, // 用户自定义选项（如压缩参数）
nullptr, nullptr // 进度条和参数（暂不启用）
);
// 检查转换结果
bool success = (poENVIDataset != nullptr);
if (!success) {
qDebug() << "错误：ENVI文件创建失败" ;
}
else {
GDALClose(poENVIDataset);
}
// 关闭输入数据集
GDALClose(poVRTDataset);
return success;
}
// 遥感类
ErrorCode ResampleDEM(QString indemPath, QString outdemPath, double gridx, double gridy)
{
double gridlat = gridy / 110000.0;
double gridlon = gridx / 100000.0;
long espgcode = GetEPSGFromRasterFile(indemPath.toUtf8().constData());
if (espgcode == 4326) {
resampleRaster(indemPath.toUtf8().constData(), outdemPath.toUtf8().constData(), gridlon, gridlat);
return ErrorCode::SUCCESS;
}
else {
QMessageBox::warning(nullptr, u8"警告", u8"请输入WGS84坐标的DEM影像");
return ErrorCode::FAIL;
}
}
void BASECONSTVARIABLEAPI CloseAllGDALRaster()
{
GDALDestroyDriverManager();
return ;
}
ErrorCode Complex2AmpRaster(QString inComplexPath, QString outRasterPath)
{
gdalImageComplex inimg(inComplexPath);
gdalImage ampimg = CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);
long blocklines = Memory1GB / 8 / inimg.width;
blocklines = blocklines < 100 ? 100 : blocklines;
omp_lock_t lock;
omp_init_lock(&lock);
#pragma omp parallel for
for (int64_t startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {
Eigen::MatrixXd imgArrb1 = ampimg.getData(startrow, 0, blocklines, inimg.width, 1);
Eigen::MatrixXcd imgArr = inimg.getDataComplex(startrow, 0, blocklines, inimg.width, 1);
imgArrb1 = imgArr.array().abs();
omp_set_lock(&lock);
ampimg.saveImage(imgArrb1, startrow, 0, 1);
omp_unset_lock(&lock); //
}
omp_destroy_lock(&lock); //
qDebug() << u8"影像写入到：" << outRasterPath;
return ErrorCode::SUCCESS;
}
ErrorCode Complex2PhaseRaster(QString inComplexPath, QString outRasterPath)
{
gdalImageComplex inimg(inComplexPath);
gdalImage ampimg = CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);
long blocklines = Memory1MB / 8 / inimg.width*200;
blocklines = blocklines < 100 ? 100 : blocklines;
;
omp_lock_t lock;
omp_init_lock(&lock);
#pragma omp parallel for
for (int64_t startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {
Eigen::MatrixXd imgArrb1 = ampimg.getData(startrow, 0, blocklines, inimg.width, 1);
Eigen::MatrixXcd imgArr = inimg.getData(startrow, 0, blocklines, inimg.width, 1);
imgArrb1 = imgArr.array().arg();
omp_set_lock(&lock);
ampimg.saveImage(imgArrb1, startrow, 0, 1);
omp_unset_lock(&lock); //
}
omp_destroy_lock(&lock); //
qDebug() << "影像写入到：" << outRasterPath;
return ErrorCode::SUCCESS;
}
ErrorCode Complex2dBRaster(QString inComplexPath, QString outRasterPath)
{
gdalImageComplex inimg(inComplexPath);
gdalImage ampimg = CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);
long blocklines = Memory1MB / 8 / inimg.width * 200;
blocklines = blocklines < 100 ? 100 : blocklines;
omp_lock_t lock;
omp_init_lock(&lock);
#pragma omp parallel for
for (int64_t startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {
Eigen::MatrixXd imgArrb1 = ampimg.getData(startrow, 0, blocklines, inimg.width, 1);
Eigen::MatrixXcd imgArr = inimg.getData(startrow, 0, blocklines, inimg.width, 1);
imgArrb1 = imgArr.array().abs().log10() * 20.0;
omp_set_lock(&lock);
ampimg.saveImage(imgArrb1, startrow, 0, 1);
omp_unset_lock(&lock); //
}
omp_destroy_lock(&lock); //
qDebug() << "影像写入到：" << outRasterPath;
return ErrorCode::SUCCESS;
}
ErrorCode BASECONSTVARIABLEAPI amp2dBRaster(QString inPath, QString outRasterPath)
{
gdalImage inimg(inPath);
gdalImage dBimg = CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);
long blocklines = Memory1MB / 8 / inimg.width * 200;
blocklines = blocklines < 100 ? 100 : blocklines;
omp_lock_t lock;
omp_init_lock(&lock);
#pragma omp parallel for
for (int64_t startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {
Eigen::MatrixXd imgArrdB = dBimg.getData(startrow, 0, blocklines, inimg.width, 1);
Eigen::MatrixXd imgArr = inimg.getData(startrow, 0, blocklines, inimg.width, 1);
imgArrdB = imgArr.array().log10() * 20.0;
omp_set_lock(&lock);
dBimg.saveImage(imgArrdB, startrow, 0, 1);
omp_unset_lock(&lock); //
}
omp_destroy_lock(&lock); //
qDebug() << "影像写入到：" << outRasterPath;
return ErrorCode::SUCCESS;
}
void MultiLookRaster(QString inRasterPath, QString outRasterPath, long looklineNumrow, long looklineNumCol)
{
gdalImage inRaster(inRasterPath);
int outRows = inRaster.height / looklineNumrow;
int outCols = inRaster.width / looklineNumCol;
int bandnum = 1;
QString project = inRaster.projection;
Eigen::MatrixXd inRasterGt = inRaster.gt;
GDALDataType datetype = inRaster.getDataType();
gdalImage outRaster = CreategdalImage(outRasterPath, outRows, outCols, 1, inRasterGt, project, true, true, true, datetype);
Eigen::MatrixXd inRasterArr = inRaster.getData(0, 0, inRaster.height, inRaster.width, 1);
Eigen::MatrixXd outRasterArr = outRaster.getData(0, 0, outRows, outCols, 1);
// 多视处理
#pragma omp parallel for collapse(2)
for (int row = 0; row < outRows; ++row) {
for (int col = 0; col < outCols; ++col) {
// 计算输入矩阵的窗口位置[2,3](@ref)
const int inputRow = row * looklineNumrow;
const int inputCol = col * looklineNumCol;
// 动态计算实际窗口大小（处理边界情况）[3](@ref)
const int actualRows = (row == outRows - 1) ?
(inRaster.height - inputRow) : looklineNumrow;
const int actualCols = (col == outCols - 1) ?
(inRaster.width - inputCol) : looklineNumCol;
// 提取当前窗口数据块[7](@ref)
Eigen::MatrixXd window = inRasterArr.block(inputRow, inputCol, actualRows, actualCols);
// 计算多视平均值（自动处理数据类型转换）[2,7](@ref)
double sum = 0.0;
if constexpr (std::is_integral_v<typename Eigen::MatrixXd::Scalar>) {
sum = window.cast<double>().sum(); // 整型数据转换为双精度计算
}
else {
sum = window.sum();
}
outRasterArr(row, col) = sum / (actualRows * actualCols);
}
}
// 可选：处理残留边界（当输入尺寸不是整数倍时）
if ((inRaster.height % looklineNumrow != 0) || (inRaster.width % looklineNumCol != 0)) {
// 底部边界处理[3](@ref)
const int residualRowStart = outRows * looklineNumrow;
const int residualRowSize = inRaster.height - residualRowStart;
if (residualRowSize > 0) {
#pragma omp parallel for
for (int col = 0; col < outCols; ++col) {
const int inputCol = col * looklineNumCol;
const int actualCols = (col == outCols - 1) ?
(inRaster.width - inputCol) : looklineNumCol;
Eigen::MatrixXd window = inRasterArr.block(residualRowStart, inputCol,
residualRowSize, actualCols);
outRasterArr(outRows - 1, col) = window.mean();
}
}
// 右侧边界处理[3](@ref)
const int residualColStart = outCols * looklineNumCol;
const int residualColSize = inRaster.width - residualColStart;
if (residualColSize > 0) {
#pragma omp parallel for
for (int row = 0; row < outRows; ++row) {
const int inputRow = row * looklineNumrow;
const int actualRows = (row == outRows - 1) ?
(inRaster.height - inputRow) : looklineNumrow;
Eigen::MatrixXd window = inRasterArr.block(inputRow, residualColStart,
actualRows, residualColSize);
outRasterArr(row, outCols - 1) = window.mean();
}
}
}
// 保存结果
outRaster.saveImage(outRasterArr, 0, 0, 1);
}
/* 启动下面的函数，会导致编译器错误，可能与Eigen 库 本身的bug相关，不建议排查，太费时间，而且不一定能排查出来
*
***
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;
}
****/