GDAL:影像RPC和GCP校正
程序员文章站
2022-07-14 14:51:19
...
李明录大哥的代码
这几个代码：
链接：https://pan.baidu.com/s/1kLZw5KGIpIi4DWLNl-vwxQ
提取码：xj06
//ImageCorrectionByRPCorGCP.h
#include "gdal_priv.h"
#include"gdalwarper.h"
#include "gdal_alg_priv.h"
#include "ogr_spatialref.h"
#include "ogrsf_frmts.h"
#include "process.h"
#include<iostream>
using namespace std;
/*!成功执行*/
const int RE_SUCCESS = 0;
/*!文件不存在*/
const int RE_FILENOTEXIST = 1;
/*!文件格式不被支持*/
const int RE_FILENOTSUPPORT = 2;
/*!图像数据类型不正确*/
const int RE_FILETYPEERROR = 3;
/*!创建图像失败*/
const int RE_CREATEFAILED = 4;
/*!输入参数错误*/
const int RE_PARAMERROR = 5;
/*!其他错误*/
const int RE_FAILED = 6;
/*!图像不存在公共区域*/
const int RE_NOSAMEEXTENT = 7;
/*!用户取消操作*/
const int RE_USERCANCEL = 8;
/*!文件已经被使用*/
const int RE_FILEISUESED = 9;
/*!不支持的像素深度*/
const int RE_DEPTHNOTSUPPORT = 10;
/*!波段数量不符合要求*/
const int RE_BANDCOUNTERROR = 11;
/*!文件不存在投影*/
const int RE_NOPROJECTION = 12;
/*!投影不一致*/
const int RE_PROJECTIONDIFF = 13;
class ImageCorrectbyRPCorGCP{
public:
/**
* @brief RPC校正
* @param pszSrcFile 输入文件路径
* @param pszDstFile 输出文件路径
* @param psRPC RPC信息
* @param dfPixErrThreshold RPC坐标反算迭代阈值
* @param papszOptions RPC校正高程设置参数
* @param eResampleMethod 重采样方式，具体参考GDALResampleAlg的定义
* @param pszFormat 输出文件格式，详细参考GDAL支持数据类型
* @param pProcess 进度条指针
* @return 返回值，表示计算过程中出现的各种错误信息
*/
int ImageWarpByRPC(const char * pszSrcFile, const char * pszDstFile, GDALRPCInfo *psRPC,
double dfPixErrThreshold = 0.0, char **papszOptions = NULL,
GDALResampleAlg eResampleMethod = GRA_NearestNeighbour, const char * pszFormat = "GTiff", CProcessBase* pProcess = NULL);
/**
* @brief TPS校正
* @param pszSrcFile 输入文件路径
* @param pszDstFile 输出文件路径
* @param nGCPCount GCP点对个数
* @param pasGCPList GCP点对列表
* @param pszDstWKT 输出文件投影，WKT格式
* @param dResX 输出图像X方向分辨率，默认为0，将输出分辨率设置为源图像分辨率
* @param dResY 输出图像Y方向分辨率，默认为0，将输出分辨率设置为源图像分辨率
* @param eResampleMethod 重采样方式，具体参考GDALResampleAlg的定义
* @param pszFormat 输出文件格式，详细参考GDAL支持数据类型
* @param pProcess 进度条指针
* @return 返回值，表示计算过程中出现的各种错误信息
*/
int ImageWarpByTPS(const char * pszSrcFile, const char * pszDstFile, int nGCPCount,
const GDAL_GCP *pasGCPList, const char * pszDstWKT, double dResX=0.0, double dResY=0.0,
GDALResampleAlg eResampleMethod = GRA_NearestNeighbour, const char * pszFormat = "GTiff", CProcessBase* pProcess = NULL);
private:
protected:
};
//ImageCorrectionByRPCorGCP.cpp
#include "ImageCorrectionByRPCorGCP.h"
/**
* @brief 图像校正坐标转换结构体
*/
struct TransformChain
{
/*! GDAL坐标转换函数指针 */
GDALTransformerFunc GDALTransformer;
/*! GDAL坐标转换参数 */
void * GDALTransformerArg;
/*! 输出图像6参数 */
double adfGeotransform[6];
/*! 输出图像逆6参数 */
double adfInvGeotransform[6];
};
/**
* @brief 创建一个转换函数
*
* 该函数用于创建一个将输出图像投影坐标转到输出图像的行列号坐标
* @sa DestroyGeoToPixelTransform
* @returns 返回GDAL坐标转换回调函数的参数 \ref GeoToPixelTransform
*/
void *CreateGeoToPixelTransform(GDALTransformerFunc GDALTransformer, void *GDALTransformerArg, double *padfGeotransform)
{
TransformChain *pChain = new TransformChain;
pChain->GDALTransformer = GDALTransformer;
pChain->GDALTransformerArg = GDALTransformerArg;
memcpy(pChain->adfGeotransform, padfGeotransform, sizeof(double) * 6);
if (!GDALInvGeoTransform(pChain->adfGeotransform, pChain->adfInvGeotransform))
{
// 如果输出图像六参数不能计算逆六参数，则不能进行变换处理
delete pChain;
return NULL;
}
return (void*)pChain;
}
/**
* @brief 析构转换参数
*/
void DestroyGeoToPixelTransform(void *GeoToPixelTransfomArg)
{
delete static_cast<TransformChain *>(GeoToPixelTransfomArg);
}
/**
* @brief 坐标转换函数
* @sa CreateGeoToPixelTransform
*/
int GeoToPixelTransform(void *pTransformerArg, int bDstToSrc, int nPointCount,
double *x, double *y, double *z, int *panSuccess)
{
TransformChain *pChain = static_cast<TransformChain*>(pTransformerArg);
if (pChain == NULL)
return FALSE;
if (!bDstToSrc) //坐标正变换
{
// 先调用GDAL库中的坐标转换函数，将坐标转为输出图像的投影坐标
if (!pChain->GDALTransformer(pChain->GDALTransformerArg, bDstToSrc, nPointCount, x, y, z, panSuccess))
return FALSE;
// 再从输出图像投影坐标系转到输出图像行列号
#pragma omp parallel for
for (int i = 0; i < nPointCount; ++i)
{
if (!panSuccess[i])
continue;
double xP = x[i];
double yP = y[i];
GDALApplyGeoTransform(pChain->adfInvGeotransform, xP, yP, &x[i], &y[i]);
}
}
else //坐标逆变换
{
// 先从输出图像行列号转到输出图像投影坐标系
#pragma omp parallel for
for (int i = 0; i < nPointCount; ++i)
{
double P = x[i];
double L = y[i];
GDALApplyGeoTransform(pChain->adfGeotransform, P, L, &x[i], &y[i]);
}
// 再调用GDAL库中坐标转换函数从输出图像投影坐标转换到原始坐标
if (!pChain->GDALTransformer(pChain->GDALTransformerArg, bDstToSrc, nPointCount, x, y, z, panSuccess))
return FALSE;
}
return TRUE;
}
/**
* @brief RPC校正
* @param pszSrcFile 输入文件路径
* @param pszDstFile 输出文件路径
* @param psRPC RPC信息
* @param dfPixErrThreshold RPC坐标反算迭代阈值
* @param papszOptions RPC校正高程设置参数
* @param eResampleMethod 重采样方式，具体参考GDALResampleAlg的定义
* @param pszFormat 输出文件格式，详细参考GDAL支持数据类型
* @param pProcess 进度条指针
* @return 返回值，表示计算过程中出现的各种错误信息
*/
int ImageCorrectbyRPCorGCP::ImageWarpByRPC(const char * pszSrcFile, const char * pszDstFile, GDALRPCInfo *psRPC,
double dfPixErrThreshold, char **papszOptions,
GDALResampleAlg eResampleMethod, const char * pszFormat, CProcessBase* pProcess){
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
// 打开原始图像并计算图像信息------------------------------------------------
GDALDataset *hSrcDS = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly);
if (hSrcDS == NULL){
return RE_FILENOTEXIST;
}
if (pProcess != NULL){
pProcess->ReSetProcess();
pProcess->SetMessage("Start RPC Correction....");
}
GDALDataType eDT = hSrcDS->GetRasterBand(1)->GetRasterDataType();//获取数据类型
int iBandCount = hSrcDS->GetRasterCount();//获取波段数
// 创建RPC坐标转换关系
void *hTransformArg = GDALCreateRPCTransformer(psRPC, FALSE, dfPixErrThreshold, papszOptions);//false 坐标反算
if (hTransformArg==NULL)
{
GDALClose(hSrcDS);
return RE_PARAMERROR;
}
// 使用SuggestedWarpOutput2函数计算输出图像四至范围、大小、六参数等信息
double adfGeoTransform[6] = { 0 };
double adfExtent[4] = { 0 };
int nPixels, nLines;
CPLErr cplerror = GDALSuggestedWarpOutput2(hSrcDS, GDALRPCTransform, hTransformArg, adfGeoTransform, &nPixels, &nLines, adfExtent, 0);
if (cplerror!=CE_None)
{
GDALClose(hSrcDS);
return RE_PARAMERROR;
}
// 创建输出图像
GDALDriver *hDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
GDALDataset *hDstDS = hDriver->Create(pszDstFile, nPixels, nLines, iBandCount, eDT, NULL);
hDstDS->SetGeoTransform(adfGeoTransform);
//设置投影方式1
// OGRSpatialReference osrs;
// osrs.SetWellKnownGeogCS("WGS84");
// char *pszWKT = NULL;
// osrs.exportToWkt(&pszWKT);
// hSrcDS->SetProjection(pszWKT);
//设置投影方式2
const char *WKT_WGS84 = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
hDstDS->SetProjection(WKT_WGS84);
//设置投影方式3
//hDstDS->SetProjection(SRS_WKT_WGS84);
// 构造GDALWarp的变换选项
GDALWarpOptions *psWO = GDALCreateWarpOptions();
psWO->papszWarpOptions = CSLDuplicate(NULL);
psWO->hSrcDS = hSrcDS;
psWO->hDstDS = hDstDS;
psWO->eWorkingDataType = eDT;
psWO->eResampleAlg = eResampleMethod;
psWO->pfnTransformer = GeoToPixelTransform;
psWO->pTransformerArg = CreateGeoToPixelTransform(GDALRPCTransform, hTransformArg, adfGeoTransform);
psWO->pfnProgress = GDALProgressFunc(ALGTermProgress);
psWO->pProgressArg = pProcess;
psWO->nBandCount = iBandCount;
psWO->panSrcBands = (int *)CPLMalloc(psWO->nBandCount*sizeof(int));
psWO->panDstBands = (int *)CPLMalloc(psWO->nBandCount*sizeof(int));
//psWO->panSrcBands = new int[iBandCount];
//psWO->panDstBands = new int[iBandCount];
for (int i = 0; i < iBandCount; i++)
{
psWO->panSrcBands[i] = i + 1;
psWO->panDstBands[i] = i + 1;
}
// 创建GDALWarp执行对象，并使用GDALWarpOptions来进行初始化
GDALWarpOperation oWO;
oWO.Initialize(psWO);
// 执行处理
//oWO.ChunkAndWarpImage(0, 0, nPixels, nLines);//1.X偏移量2.Y偏移量3.目标影像的宽度4.目标影像的高度
oWO.ChunkAndWarpMulti(0, 0, nPixels, nLines);//多线程处理
// 释放资源和关闭文件
DestroyGeoToPixelTransform(psWO->pTransformerArg);
GDALDestroyWarpOptions(psWO);
GDALDestroyRPCTransformer(hTransformArg);
GDALClose(hSrcDS);
GDALClose(hDstDS);
if (pProcess != NULL)
pProcess->SetMessage("RPC校正完成！");
return RE_SUCCESS;
}
/**
* @brief TPS校正
* @param pszSrcFile 输入文件路径
* @param pszDstFile 输出文件路径
* @param nGCPCount GCP点对个数
* @param pasGCPList GCP点对列表
* @param pszDstWKT 输出文件投影，WKT格式
* @param dResX 输出图像X方向分辨率，默认为0，将输出分辨率设置为源图像分辨率
* @param dResY 输出图像Y方向分辨率，默认为0，将输出分辨率设置为源图像分辨率
* @param eResampleMethod 重采样方式，具体参考GDALResampleAlg的定义
* @param pszFormat 输出文件格式，详细参考GDAL支持数据类型
* @param pProcess 进度条指针
* @return 返回值，表示计算过程中出现的各种错误信息
*/
int ImageCorrectbyRPCorGCP::ImageWarpByTPS(const char * pszSrcFile, const char * pszDstFile, int nGCPCount,
const GDAL_GCP *pasGCPList, const char * pszDstWKT, double dResX, double dResY,
GDALResampleAlg eResampleMethod, const char * pszFormat, CProcessBase* pProcess)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
// 打开原始图像并计算图像信息------------------------------------------------
GDALDataset *hSrcDS = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly);
if (hSrcDS == NULL){
return RE_FILENOTEXIST;
}
if (pProcess != NULL){
pProcess->ReSetProcess();
pProcess->SetMessage("Start GCP Correction....");
}
GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1)); //获取数据类型
int iBandCount = GDALGetRasterCount(hSrcDS);
// 创建TPS坐标转换关系
void *hTransformArg = GDALCreateTPSTransformer(nGCPCount, pasGCPList, FALSE);
if (hTransformArg == NULL)
{
GDALClose(hSrcDS);
return RE_PARAMERROR;
}
// 使用SuggestedWarpOutput函数计算输出图像四至范围、大小、六参数等信息
double adfGeoTransform[6];
double adfExtent[4];
int nPixels, nLines;
if (GDALSuggestedWarpOutput2(hSrcDS, GDALTPSTransform, hTransformArg,
adfGeoTransform, &nPixels, &nLines, adfExtent, 0) != CE_None)
{
GDALClose(hSrcDS);
return RE_PARAMERROR;
}
// 下面开始根据用户指定的分辨率来反算输出图像的大小和六参数等信息
double dResXSize = dResX;
double dResYSize = dResY;
//如果结果投影为投影坐标系统，则将其分辨率与原始影像一致
// if (dResXSize == 0.0 && dResYSize == 0.0)
// {
// double dGeoTrans[6] = { 0 };
// GDALGetGeoTransform(hSrcDS, dGeoTrans);
// dResXSize = ABS(dGeoTrans[1]);
// dResYSize = ABS(dGeoTrans[5]);
// }
// 如果用户指定了输出图像的分辨率
if (dResXSize != 0.0 || dResYSize != 0.0)
{
// 如果只指定了一个，使用自动计算的结果
if (dResXSize == 0.0) dResXSize = adfGeoTransform[1];
if (dResYSize == 0.0) dResYSize = adfGeoTransform[5];
// 确保分辨率符号正确
if (dResXSize < 0.0) dResXSize = -dResXSize;
if (dResYSize > 0.0) dResYSize = -dResYSize;
// 计算输出图像的范围
double minX = adfGeoTransform[0];
double maxX = adfGeoTransform[0] + adfGeoTransform[1] * nPixels;
double maxY = adfGeoTransform[3];
double minY = adfGeoTransform[3] + adfGeoTransform[5] * nLines;
// 按照用户指定的分辨率来计算图像的输出大小以及范围
nPixels = (int)(((maxX - minX) / dResXSize) + 0.5);
nLines = (int)(((minY - maxY) / dResYSize) + 0.5);
adfGeoTransform[0] = minX;
adfGeoTransform[3] = maxY;
adfGeoTransform[1] = dResXSize;
adfGeoTransform[5] = dResYSize;
}
// 创建输出图像
GDALDriver* hDiver = GetGDALDriverManager()->GetDriverByName(pszFormat);
GDALDataset *hDstDS = hDiver->Create(pszDstFile, nPixels, nLines, iBandCount, eDT, NULL);
hDstDS->SetProjection(pszDstWKT);
hDstDS->SetGeoTransform(adfGeoTransform);
// 构造GDALWarp的变换选项
GDALWarpOptions *psWO = GDALCreateWarpOptions();
psWO->papszWarpOptions = CSLDuplicate(NULL);
psWO->eWorkingDataType = eDT;
psWO->eResampleAlg = eResampleMethod;
psWO->hSrcDS = hSrcDS;
psWO->hDstDS = hDstDS;
psWO->pfnTransformer = GeoToPixelTransform;
psWO->pTransformerArg = CreateGeoToPixelTransform(GDALTPSTransform, hTransformArg, adfGeoTransform);
psWO->pfnProgress = GDALProgressFunc(ALGTermProgress);
psWO->pProgressArg = pProcess;
psWO->nBandCount = iBandCount;
psWO->panSrcBands = (int *)CPLMalloc(psWO->nBandCount*sizeof(int));
psWO->panDstBands = (int *)CPLMalloc(psWO->nBandCount*sizeof(int));
for (int i = 0; i < iBandCount; i++)
{
psWO->panSrcBands[i] = i + 1;
psWO->panDstBands[i] = i + 1;
}
// 创建GDALWarp执行对象，并使用GDALWarpOptions来进行初始化
GDALWarpOperation oWO;
oWO.Initialize(psWO);
// 执行处理
oWO.ChunkAndWarpImage(0, 0, nPixels, nLines);//单线程处理
//oWO.ChunkAndWarpMulti(0, 0, nPixels, nLines);//多线程处理
// 释放资源和关闭文件
DestroyGeoToPixelTransform(psWO->pTransformerArg);
GDALDestroyWarpOptions(psWO);
GDALDestroyTPSTransformer(hTransformArg);
GDALClose(hSrcDS);
GDALClose(hDstDS);
if (pProcess != NULL)
pProcess->SetMessage("TPS校正完成！");
return RE_SUCCESS;
}
//process.h
#include <limits>
#include <string>
using namespace std;
#ifndef MAX
#define MAX(a, b) (((a) > (b)) ? (a) : (b))
#endif // !1
#ifndef MIN
#define MIN(a, b) (((a) > (b)) ? (b) : (a))
#endif // !1
class CProcessBase
{
public:
/**
*@brief构造函数
*/
CProcessBase()
{
m_dPosition = 0.0;
m_iStepCount = 100;
m_iCurStep = 0;
m_bIsContinue = true;
}
/**
*@brief析构函数
*/
virtual~CProcessBase(){}
/**
*@brief设置进度信息
*@parampszMsg 进度信息
*/
virtual void SetMessage(const char *pszMsg) = 0;
/**
*@brief设置进度值
*@paramdPosition 进度值
*@return 返回是否取消的状态，true为不取消，false为取消
*/
virtual bool SetPosition(double dPosition) = 0;
/**
*@brief进度条前进一步，返回true表示继续，false表示取消
*@return 返回是否取消的状态，true为不取消，false为取消
*/
virtual bool StepIt() = 0;
/**
*@brief设置进度个数
*@paramiStepCount 进度个数
*/
virtual void SetStepCount(int iStepCount)
{
ReSetProcess();
m_iStepCount = iStepCount;
}
/**
*@brief获取进度信息
*@return 返回当前进度信息
*/
string GetMessage()
{
return m_strMessage;
}
/**
*@brief获取进度值
*@return 返回当前进度值
*/
double GetPosition()
{
return m_dPosition;
}
/**
*@brief重置进度条
*/
void ReSetProcess()
{
m_dPosition = 0.0;
m_iStepCount = 100;
m_iCurStep = 0;
m_bIsContinue = true;
}
/*!进度信息*/
string m_strMessage;
/*!进度值*/
double m_dPosition;
/*!进度个数*/
int m_iStepCount;
/*!进度当前个数*/
int m_iCurStep;
/*!是否取消，值为false时表示计算取消*/
bool m_bIsContinue;
};
/**
*@brief控制台进度条类
*
*提供控制台程序的进度条类接口，来反映当前算法的进度值
*/
class CConsoleProcess : public CProcessBase
{
public:
/**
*@brief构造函数
*/
CConsoleProcess()
{
m_dPosition = 0.0;
m_iStepCount = 100;
m_iCurStep = 0;
};
/**
*@brief析构函数
*/
~CConsoleProcess()
{
//remove(m_pszFile);
};
/**
*@brief设置进度信息
*@parampszMsg 进度信息
*/
void SetMessage(const char* pszMsg)
{
m_strMessage = pszMsg;
printf("%s\n", pszMsg);
}
/**
*@brief设置进度值
*@paramdPosition 进度值
*@return返回是否取消的状态，true为不取消，false为取消
*/
bool SetPosition(double dPosition)
{
m_dPosition = dPosition;
TermProgress(m_dPosition);
m_bIsContinue = true;
return true;
}
/**
*@brief进度条前进一步
*@return返回是否取消的状态，true为不取消，false为取消
*/
bool StepIt()
{
m_iCurStep++;
m_dPosition = m_iCurStep*1.0 / m_iStepCount;
TermProgress(m_dPosition);
m_bIsContinue = true;
return true;
}
private:
void TermProgress(double dfComplete)
{
static int nLastTick = -1;
int nThisTick = (int)(dfComplete*40.0);
nThisTick = MIN(40, MAX(0, nThisTick));
if (nThisTick < nLastTick && nLastTick >= 39)
nLastTick = -1;
if (nThisTick <= nLastTick)
return;
while (nThisTick > nLastTick)
{
nLastTick++;
if (nLastTick % 4 == 0)
fprintf(stdout, "%d", (nLastTick / 4) * 10);
else
fprintf(stdout, ".");
}
if (nThisTick == 40)
fprintf(stdout, "-done.\n");
else
fflush(stdout);
}
};
/**
*\brief调用GDAL进度条接口
*
*该函数用于将GDAL算法中的进度信息导出到CProcessBase基类中，供给界面显示
*
*@paramdfComplete 完成进度值，其取值为0.0到1.0之间
*@parampszMessage 进度信息
*@parampProgressArgCProcessBase的指针
*
*@return返回TRUE表示继续计算，否则为取消
*/
inline int ALGTermProgress(double dfComplete, const char* pszMessage, void* pProgressArg)
{
if (pProgressArg != NULL)
{
CProcessBase*pProcess = (CProcessBase*)pProgressArg;
pProcess->m_bIsContinue = pProcess->SetPosition(dfComplete);
if (pProcess->m_bIsContinue)
return TRUE;
else
return FALSE;
}
else
return TRUE;
}
//调用
#include "ImageCorrectionByRPCorGCP.h"
void CreateGCPList(GDAL_GCP **papasGCPs, int nGCPCount);
int main(){
/********************************************************************************************************/
//------------------------------------控制点影像校正-----------------------------------------------------//
/********************************************************************************************************/
const char* pszSrcFile1 = "C:/Users/Administrator/Documents/Tencent Files/328302992/FileRecv/GDAL书籍代码及数据/GDAL书籍数据/7-Data/orb_view3_beijing.tif";
const char* pszDstFile1 = "C:/Users/Administrator/Documents/Tencent Files/328302992/FileRecv/GDAL书籍代码及数据/GDAL书籍数据/7-Data/orb_view3_beijing_校正后.tif";
// 构造GCP列表
int nGCPCount = 14;
GDAL_GCP *pasGCPs = NULL;
CreateGCPList(&pasGCPs, nGCPCount);
// 构造输出文件投影，通过OGRSpatialReference类将EPSG代码专为WKT字符串
OGRSpatialReference OSR;
OSR.importFromEPSG(900913);
char *pszSRSWKT = NULL;
OSR.exportToWkt(&pszSRSWKT);
//开始进行控制点校正
CConsoleProcess *pProGCP = new CConsoleProcess;
ImageCorrectbyRPCorGCP ImageCorrectbyTPS;
ImageCorrectbyTPS.ImageWarpByTPS(pszSrcFile1, pszDstFile1, nGCPCount, pasGCPs, pszSRSWKT, 0.0, 0.0, GRA_Cubic, "GTiff", pProGCP);
delete pProGCP;
/********************************************************************************************************/
//------------------------------------RPC影像校正-------------------------------------------------------//
/********************************************************************************************************/
const char *pszSrcFile2 = "F:/测试数据/GF2/GF2_PMS2_E106.4_N38.3_20171210_L1A0002836848-MSS2.tiff";
const char *pszDstFile2 = "F:/测试数据/GF2/GF2_PMS2_E106.4_N38.3_20171210_L1A0002836848-MSS2_校正.tiff";
// 首先打开图像获取RPC信息，并构造为一个GDALRPCInfo结构体
// 如果RPC信息存储在文件中，可以使用GDALLoadRPCFile()和GDALLoadRPBFile()着两个函数来读取
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
GDALDataset *pSrcDS = (GDALDataset*)GDALOpen(pszSrcFile2, GA_ReadOnly);
char **papszRPCInfo = pSrcDS->GetMetadata("RPC");
GDALRPCInfo sRPC;
GDALExtractRPCInfo(papszRPCInfo, &sRPC);
GDALClose(pSrcDS);
// 设置RPC校正中用到的高程信息，如果不设置，高程值全部为0
char **papszRPCOptions = NULL;
// papszRPCOptions = CSLSetNameValue(papszRPCOptions, "RPC_HEIGHT", "10.5");
// papszRPCOptions = CSLSetNameValue(papszRPCOptions, "RPC_HEIGHT_SCALE", "1.0");
// papszRPCOptions = CSLSetNameValue(papszRPCOptions, "RPC_DEM", "F:\\DEM.tif");
// papszRPCOptions = CSLSetNameValue(papszRPCOptions, "RPC_DEMINTERPOLATION", "cubic");
CConsoleProcess *pProRPC = new CConsoleProcess;
ImageCorrectbyRPCorGCP ImageCorrectbyRPC;
ImageCorrectbyRPC.ImageWarpByRPC(pszSrcFile2, pszDstFile2, &sRPC, 0.0, papszRPCOptions, GRA_Cubic, "GTiff", pProRPC);
delete pProRPC;
return 0;
}
void CreateGCPList(GDAL_GCP **papasGCPs, int nGCPCount)
{
GDAL_GCP *pasGCPs = (GDAL_GCP *)CPLMalloc(sizeof(GDAL_GCP) * nGCPCount);
GDALInitGCPs(nGCPCount, pasGCPs);
pasGCPs[0].dfGCPPixel = 490.686;
pasGCPs[0].dfGCPLine = 549.4;
pasGCPs[0].dfGCPX = 1.2956e+07;
pasGCPs[0].dfGCPY = 4.85454e+06;
pasGCPs[1].dfGCPPixel = 1356.07;
pasGCPs[1].dfGCPLine = 325.889;
pasGCPs[1].dfGCPX = 1.29571e+07;
pasGCPs[1].dfGCPY = 4.85459e+06;
pasGCPs[2].dfGCPPixel = 1562.03;
pasGCPs[2].dfGCPLine = 1353.52;
pasGCPs[2].dfGCPX = 1.29571e+07;
pasGCPs[2].dfGCPY = 4.85321e+06;
pasGCPs[3].dfGCPPixel = 792.372;
pasGCPs[3].dfGCPLine = 1555.38;
pasGCPs[3].dfGCPX = 1.29561e+07;
pasGCPs[3].dfGCPY = 4.85316e+06;
pasGCPs[4].dfGCPPixel = 1698.31;
pasGCPs[4].dfGCPLine = 125.354;
pasGCPs[4].dfGCPX = 1.29576e+07;
pasGCPs[4].dfGCPY = 4.85476e+06;
pasGCPs[5].dfGCPPixel = 690.134;
pasGCPs[5].dfGCPLine = 72.8692;
pasGCPs[5].dfGCPX = 1.29563e+07;
pasGCPs[5].dfGCPY = 4.8551e+06;
pasGCPs[6].dfGCPPixel = 56.0791;
pasGCPs[6].dfGCPLine = 663.708;
pasGCPs[6].dfGCPX = 1.29554e+07;
pasGCPs[6].dfGCPY = 4.85451e+06;
pasGCPs[7].dfGCPPixel = 295.116;
pasGCPs[7].dfGCPLine = 1833.23;
pasGCPs[7].dfGCPX = 1.29554e+07;
pasGCPs[7].dfGCPY = 4.85294e+06;
pasGCPs[8].dfGCPPixel = 1948.79;
pasGCPs[8].dfGCPLine = 1716.42;
pasGCPs[8].dfGCPX = 1.29575e+07;
pasGCPs[8].dfGCPY = 4.85264e+06;
pasGCPs[9].dfGCPPixel = 1565.59;
pasGCPs[9].dfGCPLine = 1973.35;
pasGCPs[9].dfGCPX = 1.29569e+07;
pasGCPs[9].dfGCPY = 4.85241e+06;
pasGCPs[10].dfGCPPixel = 452.785;
pasGCPs[10].dfGCPLine = 1045.22;
pasGCPs[10].dfGCPX = 1.29558e+07;
pasGCPs[10].dfGCPY = 4.85391e+06;
pasGCPs[11].dfGCPPixel = 1496.64;
pasGCPs[11].dfGCPLine = 876.708;
pasGCPs[11].dfGCPX = 1.29571e+07;
pasGCPs[11].dfGCPY = 4.85384e+06;
pasGCPs[12].dfGCPPixel = 149.061;
pasGCPs[12].dfGCPLine = 1393.46;
pasGCPs[12].dfGCPX = 1.29553e+07;
pasGCPs[12].dfGCPY = 4.85354e+06;
pasGCPs[13].dfGCPPixel = 276.186;
pasGCPs[13].dfGCPLine = 450.768;
pasGCPs[13].dfGCPX = 1.29557e+07;
pasGCPs[13].dfGCPY = 4.85473e+06;
*papasGCPs = pasGCPs;
}
RPC校正
前：
后：
GCP校正
前：
后：