268 lines
8.7 KiB
C++
268 lines
8.7 KiB
C++
|
/**
|
|||
|
|
* @file RasterProjectBase.cpp
|
|||
|
|
* @brief None
|
|||
|
|
* @author 陈增辉 (3045316072@qq.com)
|
|||
|
|
* @version 2.5.0
|
|||
|
|
* @date 24-6-4
|
|||
|
|
* @copyright Copyright (c) Since 2024 中科卫星应用研究院 All rights reserved.
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
#include <QDebug>
|
|||
|
|
#include "RasterToolBase.h"
|
|||
|
|
#include "gdal_priv.h"
|
|||
|
|
#include "ogr_spatialref.h"
|
|||
|
|
#include "cpl_conv.h" // for CPLMalloc()
|
|||
|
|
#include <QTextCodec>
|
|||
|
|
#include <iostream>
|
|||
|
|
#include <QFile>
|
|||
|
|
|
|||
|
|
namespace RasterToolBase {
|
|||
|
|
long getProjectEPSGCodeByLon_Lat(double lon, double lat, ProjectStripDelta stripState)
|
|||
|
|
{
|
|||
|
|
long EPSGCode = 0;
|
|||
|
|
switch(stripState) {
|
|||
|
|
case ProjectStripDelta::Strip_3: {
|
|||
|
|
break;
|
|||
|
|
};
|
|||
|
|
case ProjectStripDelta::Strip_6: {
|
|||
|
|
break;
|
|||
|
|
}
|
|||
|
|
default: {
|
|||
|
|
EPSGCode = -1;
|
|||
|
|
break;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
qDebug() << QString(" EPSG code : %1").arg(EPSGCode);
|
|||
|
|
return EPSGCode;
|
|||
|
|
}
|
|||
|
|
long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat)
|
|||
|
|
{
|
|||
|
|
// EPSG 4534 ~ 4554 3 度带
|
|||
|
|
// 首先判断是否是在 中国带宽范围
|
|||
|
|
// 中心经度范围 ：75E ~ 135E 实际范围 73.5E ~ 136.5E,
|
|||
|
|
// 纬度范围 3N ~ 54N，放宽到 0N~ 60N
|
|||
|
|
if(73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内
|
|||
|
|
long code = trunc((lon - 73.5) / 3) + 4534;
|
|||
|
|
return code;
|
|||
|
|
} else { // 非中国境内 使用 高斯克吕格
|
|||
|
|
bool isSouth = lat < 0; // 简单判断南北半球，这里仅为示例，实际应用可能需要更细致的逻辑
|
|||
|
|
long prefix = isSouth ? 327000 : 326000;
|
|||
|
|
// std::string prefix = isSouth ? "327" : "326";
|
|||
|
|
lon = fmod(lon + 360.0, 360.0);
|
|||
|
|
long zone = std::floor((lon + 180.0) / 3.0);
|
|||
|
|
prefix = prefix + zone;
|
|||
|
|
return prefix;
|
|||
|
|
}
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat)
|
|||
|
|
{
|
|||
|
|
// EPSG 4502 ~ 4512 6度带
|
|||
|
|
// 首先判断是否是在 中国带宽范围
|
|||
|
|
// 中心经度范围 ：75E ~ 135E 实际范围 72.0E ~ 138E,
|
|||
|
|
// 纬度范围 3N ~ 54N，放宽到 0N~ 60N
|
|||
|
|
if(73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内
|
|||
|
|
long code = trunc((lon - 72.0) / 6) + 4502;
|
|||
|
|
return code;
|
|||
|
|
} else { // 非中国境内 使用 UTM// 确定带号，6度带从1开始到60，每6度一个带
|
|||
|
|
int zone = static_cast<int>((lon + 180.0) / 6.0) + 1;
|
|||
|
|
bool isSouth = lon < 0; // 判断是否在南半球
|
|||
|
|
long epsgCodeBase = isSouth ? 32700 : 32600; // 计算EPSG代码
|
|||
|
|
long prefix = static_cast<int>(epsgCodeBase + zone);
|
|||
|
|
return prefix;
|
|||
|
|
}
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
QString GetProjectionNameFromEPSG(long epsgCode)
|
|||
|
|
{
|
|||
|
|
qDebug() << "============= GetProjectionNameFromEPSG ======================";
|
|||
|
|
OGRSpatialReference oSRS;
|
|||
|
|
|
|||
|
|
// 设置EPSG代码
|
|||
|
|
if(oSRS.importFromEPSG(epsgCode) != OGRERR_NONE) {
|
|||
|
|
qDebug() << "epsgcode not recognition";
|
|||
|
|
return "";
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// 获取并输出坐标系的描述（名称）
|
|||
|
|
const char* pszName = oSRS.GetAttrValue("GEOGCS");
|
|||
|
|
if(pszName) {
|
|||
|
|
qDebug() << "Coordinate system name for EPSG " + QString::number(epsgCode)
|
|||
|
|
<< " is: " + QString::fromStdString(pszName);
|
|||
|
|
return QString::fromStdString(pszName);
|
|||
|
|
} else {
|
|||
|
|
qDebug() << "Unable to retrieve the name for EPSG " + QString::number(epsgCode);
|
|||
|
|
return "";
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// char* wkt = NULL;
|
|||
|
|
// // 转换为WKT格式
|
|||
|
|
// oSRS.exportToWkt(&wkt);
|
|||
|
|
//
|
|||
|
|
// qDebug() << wkt;
|
|||
|
|
//
|
|||
|
|
// // 从WKT中解析投影名称，这里简化处理，实际可能需要更复杂的逻辑来准确提取名称
|
|||
|
|
// std::string wktStr(wkt);
|
|||
|
|
// long start = wktStr.find("PROJCS[\"") + 8; // 找到"PROJCS["后的第一个双引号位置
|
|||
|
|
// // 从start位置开始找下一个双引号，这之间的内容即为投影名称
|
|||
|
|
// int end = wktStr.find('\"', start);
|
|||
|
|
// QString projName = QString::fromStdString(wktStr.substr(start, end - start));
|
|||
|
|
//
|
|||
|
|
// // 释放WKT字符串内存
|
|||
|
|
// CPLFree(wkt);
|
|||
|
|
|
|||
|
|
// return projName;
|
|||
|
|
}
|
|||
|
|
long GetEPSGFromRasterFile(QString filepath)
|
|||
|
|
{
|
|||
|
|
qDebug() << "============= GetEPSGFromRasterFile ======================";
|
|||
|
|
// QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器
|
|||
|
|
// QByteArray byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray
|
|||
|
|
//,这个应该会自动释放 const char* charArray = byteArray.constData(); //
|
|||
|
|
// 获取QByteArray的const char*指针
|
|||
|
|
|
|||
|
|
{
|
|||
|
|
if(QFile(filepath).exists()){
|
|||
|
|
qDebug() << "info: the image found.\n";
|
|||
|
|
}else{
|
|||
|
|
return -1;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
GDALAllRegister();
|
|||
|
|
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注册GDAL驱动
|
|||
|
|
// std::cout<<filepath.toLocal8Bit().constData()<<std::endl;
|
|||
|
|
// 打开影像文件
|
|||
|
|
GDALDataset* poDataset;
|
|||
|
|
poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
|
|||
|
|
if(NULL==poDataset) {
|
|||
|
|
qDebug() << "Error: Unable to open the image.\n";
|
|||
|
|
return -1;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// 获取影像的投影信息
|
|||
|
|
const char* pszProjection = poDataset->GetProjectionRef();
|
|||
|
|
|
|||
|
|
qDebug() << QString::fromUtf8(pszProjection);
|
|||
|
|
|
|||
|
|
// 创建SpatialReference对象
|
|||
|
|
OGRSpatialReference oSRS;
|
|||
|
|
if(oSRS.importFromWkt((char**)&pszProjection) != OGRERR_NONE) {
|
|||
|
|
qDebug() << ("Error: Unable to import projection information.\n");
|
|||
|
|
GDALClose(poDataset);
|
|||
|
|
return -1;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
long epsgCode = atoi(oSRS.GetAuthorityCode(nullptr)); // 获取EPSG代码
|
|||
|
|
|
|||
|
|
if(epsgCode != 0) {
|
|||
|
|
GDALClose(poDataset);
|
|||
|
|
qDebug() << QString("file %1 :epsg Code %2").arg(filepath).arg(epsgCode);
|
|||
|
|
return epsgCode;
|
|||
|
|
} else {
|
|||
|
|
qDebug() << "EPSG code could not be determined from the spatial reference.";
|
|||
|
|
GDALClose(poDataset);
|
|||
|
|
return -1;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
std::shared_ptr<PointRaster> GetCenterPointInRaster(QString filepath)
|
|||
|
|
{
|
|||
|
|
qDebug() << "============= GetCenterPointInRaster ======================";
|
|||
|
|
// QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器
|
|||
|
|
// QByteArray byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray
|
|||
|
|
//,这个应该会自动释放 const char* charArray = byteArray.constData(); //
|
|||
|
|
// 获取QByteArray的const char*指针
|
|||
|
|
GDALAllRegister();
|
|||
|
|
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
|
|||
|
|
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
|
|||
|
|
// std::cout<<filepath.toLocal8Bit().constData()<<std::endl;
|
|||
|
|
GDALDataset* poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
|
|||
|
|
if(nullptr==poDataset||NULL==poDataset) {
|
|||
|
|
qDebug() << "Could not open dataset";
|
|||
|
|
return nullptr; // 表示读取失败
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double x, y, z;
|
|||
|
|
|
|||
|
|
bool flag = false;
|
|||
|
|
|
|||
|
|
{
|
|||
|
|
double adfGeoTransform[6];
|
|||
|
|
if(poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
|
|||
|
|
qDebug() << "Failed to get GeoTransform";
|
|||
|
|
return nullptr;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double dfWidth = poDataset->GetRasterXSize();
|
|||
|
|
double dfHeight = poDataset->GetRasterYSize();
|
|||
|
|
|
|||
|
|
// 计算中心点坐标（像素坐标）
|
|||
|
|
double dfCenterX = adfGeoTransform[0] + dfWidth * adfGeoTransform[1] / 2.0
|
|||
|
|
+ dfHeight * adfGeoTransform[2] / 2.0;
|
|||
|
|
double dfCenterY = adfGeoTransform[3] + dfWidth * adfGeoTransform[4] / 2.0
|
|||
|
|
+ dfHeight * adfGeoTransform[5] / 2.0;
|
|||
|
|
|
|||
|
|
OGRSpatialReference oSRS;
|
|||
|
|
oSRS.importFromWkt(poDataset->GetProjectionRef());
|
|||
|
|
|
|||
|
|
if(oSRS.IsGeographic()) {
|
|||
|
|
qDebug() << "Center coords (already in geographic): (" + QString::number(dfCenterX)
|
|||
|
|
+ ", " + QString::number(dfCenterY) + ")";
|
|||
|
|
flag = true;
|
|||
|
|
x = dfCenterX;
|
|||
|
|
y = dfCenterY;
|
|||
|
|
} else {
|
|||
|
|
// 如果不是地理坐标系，转换到WGS84
|
|||
|
|
OGRSpatialReference oSRS_WGS84;
|
|||
|
|
oSRS_WGS84.SetWellKnownGeogCS("WGS84");
|
|||
|
|
|
|||
|
|
OGRCoordinateTransformation* poCT =
|
|||
|
|
OGRCreateCoordinateTransformation(&oSRS, &oSRS_WGS84);
|
|||
|
|
if(poCT == nullptr) {
|
|||
|
|
qDebug() << "Failed to create coordinate transformation";
|
|||
|
|
return nullptr;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// double dfLon, dfLat;
|
|||
|
|
if(poCT->Transform(1, &dfCenterX, &dfCenterY)) {
|
|||
|
|
qDebug() << "Center coords (transformed to WGS84): ("
|
|||
|
|
+ QString::number(dfCenterX) + ", " + QString::number(dfCenterY)
|
|||
|
|
<< ")";
|
|||
|
|
flag = true;
|
|||
|
|
x = dfCenterX;
|
|||
|
|
y = dfCenterY;
|
|||
|
|
} else {
|
|||
|
|
qDebug() << "Transformation failed.";
|
|||
|
|
}
|
|||
|
|
OGRCoordinateTransformation::DestroyCT(poCT);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
if(nullptr==poDataset||NULL==poDataset){}else{
|
|||
|
|
GDALClose(poDataset);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if(flag) {
|
|||
|
|
std::shared_ptr<PointRaster> RasterCenterPoint = std::make_shared<PointRaster>();
|
|||
|
|
RasterCenterPoint->x = x;
|
|||
|
|
RasterCenterPoint->y = y;
|
|||
|
|
RasterCenterPoint->z = 0;
|
|||
|
|
return RasterCenterPoint;
|
|||
|
|
} else {
|
|||
|
|
return nullptr;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long epsg_code)
|
|||
|
|
{
|
|||
|
|
OGRSpatialReference oSRS;
|
|||
|
|
if(oSRS.importFromEPSG(epsg_code) == OGRERR_NONE) {
|
|||
|
|
if(oSRS.IsGeographic()) {
|
|||
|
|
return CoordinateSystemType::GeoCoordinateSystem;
|
|||
|
|
} else if(oSRS.IsProjected()) {
|
|||
|
|
return CoordinateSystemType::ProjectCoordinateSystem;
|
|||
|
|
}
|
|||
|
|
} else {
|
|||
|
|
return CoordinateSystemType::UNKNOW;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
} // namespace RasterToolBase
|