2120 lines
70 KiB
C++
2120 lines
70 KiB
C++
#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(poDataset);
|
||
GDALClose(poOutDataset);
|
||
}
|
||
|
||
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 = (adfGeoTransform[1] * nXSize) / targetPixelSizeX;
|
||
double targetYSize = std::abs((adfGeoTransform[5] * nYSize) / targetPixelSizeY);
|
||
|
||
// 创建目标数据集（输出栅格）
|
||
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
|
||
if (poDriver == nullptr) {
|
||
qDebug() << "Failed to get GTiff driver." ;
|
||
GDALClose(poDataset);
|
||
return;
|
||
}
|
||
|
||
// 创建输出数据集
|
||
GDALDataset* poOutDataset = poDriver->Create(outputRaster, (int)targetXSize, (int)targetYSize, poDataset->GetRasterCount(), GDT_Float32, nullptr);
|
||
if (poOutDataset == nullptr) {
|
||
qDebug() << "Failed to create output raster." ;
|
||
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.";
|
||
}
|
||
|
||
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)
|
||
{
|
||
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 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 * 2 / 8 / inimg.width;
|
||
blocklines = blocklines < 100 ? 100 : blocklines;
|
||
|
||
long startrow = 0;
|
||
for (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();
|
||
ampimg.saveImage(imgArrb1, startrow, 0, 1);
|
||
}
|
||
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 = Memory1GB * 2 / 8 / inimg.width;
|
||
blocklines = blocklines < 100 ? 100 : blocklines;
|
||
Eigen::MatrixXd imgArrb1 = Eigen::MatrixXd::Zero(blocklines, ampimg.width);
|
||
Eigen::MatrixXcd imgArr = Eigen::MatrixXcd::Zero(blocklines, inimg.width);
|
||
|
||
long startrow = 0;
|
||
for (startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {
|
||
|
||
imgArrb1 = ampimg.getData(startrow, 0, blocklines, inimg.width, 1);
|
||
imgArr = inimg.getData(startrow, 0, blocklines, inimg.width, 2);
|
||
imgArrb1 = imgArr.array().arg();
|
||
ampimg.saveImage(imgArrb1, startrow, 0, 1);
|
||
}
|
||
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 = Memory1GB * 2 / 8 / inimg.width;
|
||
blocklines = blocklines < 100 ? 100 : blocklines;
|
||
Eigen::MatrixXd imgArrb1 = Eigen::MatrixXd::Zero(blocklines, ampimg.width);
|
||
Eigen::MatrixXcd imgArr = Eigen::MatrixXcd::Zero(blocklines, inimg.width);
|
||
|
||
long startrow = 0;
|
||
for (startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {
|
||
|
||
imgArrb1 = ampimg.getData(startrow, 0, blocklines, inimg.width, 1);
|
||
imgArr = inimg.getData(startrow, 0, blocklines, inimg.width, 2);
|
||
imgArrb1 = imgArr.array().abs().log10() * 20.0;
|
||
ampimg.saveImage(imgArrb1, startrow, 0, 1);
|
||
}
|
||
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;
|
||
}
|
||
|
||
|
||
****/
|