modis数据的几何校正

最新推荐文章于 2024-04-16 22:11:17 发布
最新推荐文章于 2024-04-16 22:11:17 发布
阅读量6.2k 收藏 2
点赞数 1
分类专栏： 开源GIS

头文件ImageGeoLocWarp.h

#pragma once

#include "gdalwarper.h"
#include<string>
using namespace std;


/***************************************************************************
*
* Time: 2017-02-14
* Project: 遥感处理算法 
* Purpose: 地理查找表校正
* Author:  wuxin


* Describe: 地理查找表校正，用于HDF数据的校正算法
*
****************************************************************************/
#define IMAGEGEOLOCWARP_H
class  ImageGeoLocWarp
{
public:
/**
* @brief 构造函数，地理查找表校正
* @param pszSrcFile 原始数据路径
* @param pszDstFile 结果数据路径
* @param pszFormat 结果数据格式
* @param pProcess 进度条指针
*/
ImageGeoLocWarp(const char* pszSrcFile, const char* pszDstFile,
const char* pszFormat = "HFA");


~ImageGeoLocWarp(void);
/**
* @brief 获取输出参数信息
* @param iHeight 输出图像高度
* @param iWidth 输出图像宽度
* @param padfGeoTransform 输出GeoTransform六参数
* @param padfExtent 输出图像四至范围
* @return 是否计算成功，以及错误代码
*/
int GetSuggestedWarpOutput(int &iHeight, int &iWidth, double *padfGeoTransform = NULL, double *padfExtent = NULL);


/**
* @brief 执行校正算法
* @param dResX 输出图像横向分辨率，默认为0，表示自动计算输出象元大小
* @param dResY 输出图像纵向分辨率，默认为0，表示自动计算输出象元大小
* @param resampling 重采样方式，默认为双线性内插
* @return 是否计算成功，以及错误代码
*/
int DoGeoLocWarp(double dResX = 0.0, double dResY = 0.0, GDALResampleAlg resampling = GRA_Bilinear);


private:
//输入参数
const char* m_pszSrcFile; /*! 原始数据路径 */
const char* m_pszDstFile; /*! 结果数据路径 */
const char* m_pszFormat; /*! 结果数据格式 */
GDALResampleAlg m_Resampling; /*! 重采样方式 */
string m_strWkt; /*! 输出文件投影串，这里只能是WGS84经纬度 */
};

头文件IMGALG_H.h


****************************************************************************/


#ifndef IMGALG_H
#define IMGALG_H


/**
* \file ImgCore.h
* @brief 核心类型定义
*
* 导出接口（使用C语言方式），核心类型定义
*/


/**
* 忽略在MS Widows平台上的警告 lml 2010-10-19
* warning C4100: “*”: 未引用的形参
* warning C4190: “identifier1”有指定的 C 链接，但返回了与 C 不兼容的 UDT“identifier2”
* warning C4251: 类“type”需要由类“type2”的客户端使用 dll 接口
* warning C4275: 非 DLL 接口类键“identifier”作为 DLL 接口类键“identifier”的基使用
* warning C4305: 从“type1”到“type2”截断
* warning C4309: 截断常数值
* warning C4819: 该文件包含不能在当前代码页(936)中表示的字符。请将该文件保存为 Unicode 格式以防止数据丢失
* warning C4996: 使用了非标准扩展: 限定名中使用了枚举
*/
#if defined(_MSC_VER) && (_MSC_VER >= 1400)
#pragma warning(disable: 4100 4190 4251 4275 4305 4309 4819 4996 )
#endif


/**
* @brief 是否使用Vld内存泄露监测工具
*/
#if _USEVLD
#if _DEBUG //在debug模式下检测内存泄露
#include "vld.h"
#endif
#endif


/**
* @brief 是否使用LOG工具进行写日志
*/
#if _USELOG
#define USE_LOG4CPP
#endif


#include <float.h>
#include <algorithm>
#include <deque>
#include <fstream>
#include <limits>
#include <map>
#include <stack>
#include <string>
#include <vector>
using namespace std;


/**
* @brief 导出符号定义
*/
#ifdef IMGALG_EXPORTS
#define IMGALG_API __declspec(dllexport)
#else
#define IMGALG_API __declspec(dllimport)
#endif


/**
* @brief 定义NULL
*/
#ifndef NULL
#  define NULL  0
#endif


/**
* @brief 定义FALSE
*/
#ifndef FALSE
#  define FALSE 0
#endif


/**
* @brief 定义TRUE
*/
#ifndef TRUE
#  define TRUE  1
#endif


#ifndef MAX
/*! 求最大值 */
#  define MIN(a, b)      ((a<b) ? a : b)
/*! 求最小值 */
#  define MAX(a, b)      ((a>b) ? a : b)
#endif


/**
* @brief 定义ABS，求绝对值
*/
#ifndef ABS
#  define ABS(x)        ((x<0) ? (-1*(x)) : x)
#endif


/**
* @brief 定义PI=3.141592653...以及度和弧度转换
*/
#ifndef M_PI
/*! 定义圆周率PI */
# define M_PI  3.1415926535897932384626433832795
/*! 弧度转度 */
# define DEG_PER_RAD      ((double)(180.0/M_PI))
/*! 度转弧度 */
# define RAD_PER_DEG      ((double)(M_PI/180.0))
#endif


/**
* @brief 定义平方
*/
#ifndef M_SQUARE
# define M_SQUARE(x)  (x)*(x)
#endif


/**
* @brief 定义立方
*/
#ifndef M_CUBE
# define M_CUBE(x)  (x)*(x)*(x)
#endif


/*! 判断浮点数是否NaN值 */
inline bool isnan(const float& v)  { return _isnan(v) ? true : false; }
/*! 判断double数是否NaN值 */
inline bool isnan(const double& v) { return _isnan(v) ? true : false; }
/*! 获取double的NaN值 */
inline double nan() { return numeric_limits<double>::quiet_NaN(); }


/**
* @brief float类型的极值
*/
#ifndef FLT_EQUALS
/*! 浮点数是否相等 */
#define FLT_EQUALS(x, y)  (fabs((double)x-y)<FLT_EPSILON)
/*! 浮点数是否相等(指定比较阈值) */
#define FLT_EQUALS_N(x, y, z)  (fabs((double)x-y)<z)
#endif


#ifndef FLT_ZERO
/*! 浮点数是否为0 */
#define FLT_ZERO(x)  (fabs(x)<FLT_EPSILON)
#endif


/**
* @brief 释放数组
*/
#define RELEASE(x) if(x!=NULL) {delete []x; x = NULL;}


/**
* @brief 将全局区域设为操作系统默认区域
*/
#define SET_LOCAL { locale::global(locale("")); setlocale(LC_ALL,"Chinese-simplified"); }
/**
* @brief 还原全局区域设定
*/
#define REVERT_LOCAL locale::global(locale("C"))




#ifndef EQUAL
#if defined(WIN32) || defined(WIN32CE)
/*! 比较字符串是否相等 */
#  define EQUALN(a, b, n)           (_strnicmp(a, b, n) == 0)
/*! 比较字符串是否相等 */
#  define EQUAL(a, b)              (_stricmp(a, b) == 0)
#else
/*! 比较字符串是否相等 */
#  define EQUALN(a, b, n)           (strncasecmp(a, b, n) == 0)
/*! 比较字符串是否相等 */
#  define EQUAL(a, b)              (strcasecmp(a, b) == 0)
#endif
#endif


/*! byte */
typedef unsigned char byte;
/*! 8U */
typedef unsigned char DT_8U;
/*! 16U */
typedef unsigned short DT_16U;
/*! 16S */
typedef short DT_16S;
/*! 32U */
typedef unsigned int DT_32U;
/*! 32S */
typedef int DT_32S;
/*! 32F */
typedef float DT_32F;
/*! 64F */
typedef double DT_64F;


/*! 成功执行 */
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;
/*! 创建文件 */
const int RE_CREATEFILE=14;




#endif

源文件 ImageGeoLocWarp.cpp

#include "stdafx.h"
#include "ImageGeoLocWarp.h"
#include "ImageGeoLocWarp.h"//上面的头文件
#include "gdal_priv.h"
#include "gdalwarper.h"
#include "ogr_srs_api.h"
#include "IMGALG_H.h"








int RasterDelete(const char* pszFile)  
{  
GDALAllRegister();  


//打开图像  
GDALDataset *pDS = (GDALDataset *)GDALOpen(pszFile, GA_ReadOnly);  
if (pDS == NULL)  
return remove(pszFile);  


GDALDriver *pDriver = pDS->GetDriver();  
if( pDriver == NULL )  
{  
GDALClose((GDALDatasetH) pDS);  
return remove(pszFile);  
}  


GDALClose((GDALDatasetH) pDS);  


if(pDriver->Delete(pszFile) == CE_None)  
return 0;  
else  
return remove(pszFile);  
}  


ImageGeoLocWarp::ImageGeoLocWarp(const char* pszSrcFile, const char* pszDstFile,
const char* pszFormat)
{
m_pszSrcFile = pszSrcFile; /*! 原始数据路径 */
m_pszDstFile = pszDstFile; /*! 结果数据路径 */
m_pszFormat = pszFormat; /*! 结果数据格式 */

m_Resampling= GRA_Bilinear; /*! 重采样方式，默认为双线性 */
m_strWkt = SRS_WKT_WGS84; //投影只能是WGS84
}




ImageGeoLocWarp::~ImageGeoLocWarp(void)
{
}
int ImageGeoLocWarp::GetSuggestedWarpOutput(int &iHeight, int &iWidth, double *padfGeoTransform, double *padfExtent)
{
if(m_pszSrcFile == NULL || m_pszDstFile == NULL)
return RE_PARAMERROR;


GDALAllRegister();


GDALDatasetH hSrcDS = GDALOpen( m_pszSrcFile, GA_ReadOnly ); //打开文件
if(hSrcDS == NULL )
{



return RE_FILENOTEXIST;
}


if(m_strWkt.empty()) //目标投影不存在
{
GDALClose( hSrcDS );
return RE_PARAMERROR;
}


// 创建坐标转换参数
char** papszWarpOptions = NULL;
papszWarpOptions = CSLSetNameValue( papszWarpOptions, "METHOD", "GEOLOC_ARRAY" );
papszWarpOptions = CSLSetNameValue( papszWarpOptions, "DST_SRS", m_strWkt.c_str() );


void *hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszWarpOptions );
if(hTransformArg == NULL )
{



CSLDestroy( papszWarpOptions );
GDALClose( hSrcDS );


return RE_PARAMERROR;
}


GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;


// 获取目标文件大小以及转换参数信息
double adfDstGeoTransform[6] = {0};
double adfExtent[4];
int nPixels = 0, nLines = 0;


CPLErr eErr = GDALSuggestedWarpOutput2( hSrcDS, psInfo->pfnTransform, hTransformArg, 
adfDstGeoTransform, &nPixels, &nLines, adfExtent, 0 );


iWidth = nPixels;
iHeight = nLines;
if(padfGeoTransform != NULL)
memcpy(padfGeoTransform, adfDstGeoTransform, sizeof(double)*6);
if(padfExtent != NULL)
memcpy(padfExtent, adfExtent, sizeof(double)*4);


GDALDestroyGenImgProjTransformer(hTransformArg);
CSLDestroy( papszWarpOptions );
GDALClose( hSrcDS );


if(eErr != CE_None)
return RE_PARAMERROR;
else
return RE_SUCCESS;
}


int ImageGeoLocWarp::DoGeoLocWarp(double dResX, double dResY, GDALResampleAlg resampling)
{
m_Resampling == resampling;


if(m_pszSrcFile == NULL || m_pszDstFile == NULL)
return RE_PARAMERROR;





GDALAllRegister();


GDALDatasetH hSrcDS = GDALOpen( m_pszSrcFile, GA_ReadOnly ); //打开文件
if(hSrcDS == NULL )
{



return RE_FILENOTEXIST;
}


if(m_strWkt.empty()) //目标投影不存在
{



GDALClose( hSrcDS );
return RE_PARAMERROR;
}


// 创建坐标转换参数
char** papszWarpOptions = NULL;
papszWarpOptions = CSLSetNameValue( papszWarpOptions, "METHOD", "GEOLOC_ARRAY" );
papszWarpOptions = CSLSetNameValue( papszWarpOptions, "DST_SRS", m_strWkt.c_str() );


void *hTransformArg = GDALCreateGenImgProjTransformer2( hSrcDS, NULL, papszWarpOptions );
if(hTransformArg == NULL )
{



CSLDestroy( papszWarpOptions );
GDALClose( hSrcDS );


return RE_PARAMERROR;
}


GDALTransformerInfo* psInfo = (GDALTransformerInfo*)hTransformArg;


// 获取目标文件大小以及转换参数信息
double adfDstGeoTransform[6] = {0};
double adfExtent[4];
int nPixels = 0, nLines = 0;


CPLErr eErr = GDALSuggestedWarpOutput2( hSrcDS, psInfo->pfnTransform, hTransformArg, 
adfDstGeoTransform, &nPixels, &nLines, adfExtent, 0 );


if(eErr != CE_None)
{



GDALDestroyGenImgProjTransformer(hTransformArg);
CSLDestroy( papszWarpOptions );
GDALClose( hSrcDS );


return RE_PARAMERROR;
}


// 如果用户指定了输出图像的分辨率
if ( dResX != 0.0 || dResY != 0.0 )
{
// 如果只指定了一个，使用自动计算的结果
if ( dResX == 0.0 ) dResX = adfDstGeoTransform[1];
if ( dResY == 0.0 ) dResY = adfDstGeoTransform[5];


// 确保分辨率符号正确
if ( dResX < 0.0 ) dResX = -dResX;
if ( dResY > 0.0 ) dResY = -dResY;


// 计算输出图像的范围
double minX = adfDstGeoTransform[0];
double maxX = adfDstGeoTransform[0] + adfDstGeoTransform[1] * nPixels;
double maxY = adfDstGeoTransform[3];
double minY = adfDstGeoTransform[3] + adfDstGeoTransform[5] * nLines;


// 按照用户指定的分辨率来计算图像的输出大小以及范围
nPixels = ( int )((( maxX - minX ) / dResX ) + 0.5 );
nLines  = ( int )((( minY - maxY ) / dResY ) + 0.5 );
adfDstGeoTransform[0] = minX;
adfDstGeoTransform[3] = maxY;
adfDstGeoTransform[1] = dResX;
adfDstGeoTransform[5] = dResY;
}


// 获得波段数目
int nBandCount = GDALGetRasterCount(hSrcDS);
// 创建输出文件数据类型与源文件相同
GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, 1));


// 创建输出结果文件
GDALDriverH hDriver = GDALGetDriverByName( m_pszFormat);
if(hDriver == NULL )
{



GDALDestroyGenImgProjTransformer(hTransformArg);
CSLDestroy( papszWarpOptions );
GDALClose( hSrcDS );


return RE_FILENOTSUPPORT;
}


GDALDatasetH  hDstDS = GDALCreate( hDriver, m_pszDstFile, nPixels, nLines, nBandCount, eDT, NULL );
if(hDstDS == NULL)
{



GDALDestroyGenImgProjTransformer(hTransformArg);
CSLDestroy( papszWarpOptions );
GDALClose( hSrcDS );


return RE_CREATEFAILED;
}


// 设置投影等信息
GDALSetProjection( hDstDS, m_strWkt.c_str() );
GDALSetGeoTransform( hDstDS, adfDstGeoTransform );


if (hTransformArg != NULL)
GDALSetGenImgProjTransformerDstGeoTransform( hTransformArg, adfDstGeoTransform);


// 逐个波段获取颜色表
GDALColorTableH hCT;


for (int i=1; i<=nBandCount; i++)//band从1开始算
{
hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS, i) );
if( hCT != NULL )
GDALSetRasterColorTable( GDALGetRasterBand(hDstDS, i), hCT );
}


// 创建转换选项
CSLDestroy( papszWarpOptions ); //先释放之前的
papszWarpOptions = NULL;
papszWarpOptions = CSLSetNameValue(papszWarpOptions, "INIT_DEST", "0");


GDALTransformerFunc pfnTransformer = NULL;
pfnTransformer = GDALGenImgProjTransform; //坐标转换函数


GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
psWarpOptions->papszWarpOptions = CSLDuplicate(papszWarpOptions);
psWarpOptions->hSrcDS = hSrcDS;
psWarpOptions->hDstDS = hDstDS;
psWarpOptions->nBandCount = nBandCount;


// 申请内存
psWarpOptions->panSrcBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );  
psWarpOptions->panDstBands = (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );


for (int j=0; j<nBandCount; j++)
{
psWarpOptions->panSrcBands[j] = j + 1;  
psWarpOptions->panDstBands[j] = j + 1;
}


psWarpOptions->eResampleAlg == (GDALResampleAlg)resampling;

psWarpOptions->papszWarpOptions = papszWarpOptions;


psWarpOptions->pfnTransformer = pfnTransformer;
psWarpOptions->pTransformerArg = hTransformArg;


CPLErr Err = CE_None;
GDALWarpOperation oOperation;
if( oOperation.Initialize( psWarpOptions ) == CE_None )
Err = oOperation.ChunkAndWarpImage( 0, 0, nPixels, nLines);
else
{
CSLDestroy( papszWarpOptions );
GDALDestroyGenImgProjTransformer( hTransformArg );


GDALClose( hDstDS );
GDALClose( hSrcDS );


RasterDelete(m_pszDstFile); //删除结果数据
return RE_FAILED;
}


CSLDestroy( papszWarpOptions );
GDALDestroyGenImgProjTransformer( hTransformArg );


GDALClose( hDstDS );
GDALClose( hSrcDS );





if (Err != CE_None)
return RE_FAILED;


return RE_SUCCESS;
}


///*! Warp Resampling Algorithm */
//typedef enum {
//  /*! Nearest neighbour (select on one input pixel) */ GRA_NearestNeighbour=0,
//  /*! Bilinear (2x2 kernel) */                         GRA_Bilinear=1,
//  /*! Cubic Convolution Approximation (4x4 kernel) */  GRA_Cubic=2,
//  /*! Cubic B-Spline Approximation (4x4 kernel) */     GRA_CubicSpline=3,
//  /*! Lanczos windowed sinc interpolation (6x6 kernel) */ GRA_Lanczos=4
//} GDALResampleAlg;


/**
* 重采样函数(GDAL)
* @param pszSrcFile        输入文件的路径
* @param pszOutFile        写入的结果图像的路径
* @param fResX             X转换采样比，默认大小为1.0，大于1图像变大，小于1表示图像缩小
* @param fResY             Y转换采样比，默认大小为1.0
* @param nResampleMode     采样模式，有五种，具体参见GDALResampleAlg定义，默认为双线性内插
* @param pExtent           采样范围，为NULL表示计算全图
* @param pBandIndex        指定的采样波段序号，为NULL表示采样全部波段
* @param pBandCount        采样的波段个数，同pBandIndex一同使用，表示采样波段的个数
* @param pszFormat         写入的结果图像的格式
* @param pProgress         进度条指针
* @return 成功返回0，否则为其他值
*/




//int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, float fResX , float fResY, GDALResampleAlg nResampleMode,
//    LT_Envelope* pExtent, int* pBandIndex, int *pBandCount, const char* pszFormat)
//{
//    
//
//    GDALAllRegister();
//    GDALDataset *pDSrc = (GDALDataset *)GDALOpen(pszSrcFile, GA_ReadOnly);
//    if (pDSrc == NULL)
//    {
//        
//
//        return RE_FILENOTEXIST;
//    }
//
//    GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
//    if (pDriver == NULL)
//    {
//      
//
//        GDALClose((GDALDatasetH) pDSrc);
//        return RE_CREATEFILE;
//    }
//
//    int iBandCount = pDSrc->GetRasterCount();
//    string strWkt = pDSrc->GetProjectionRef();
//    GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();
//
//    double dGeoTrans[6] = {0};
//    pDSrc->GetGeoTransform(dGeoTrans);
//
//    int iNewBandCount = iBandCount;
//    if (pBandIndex != NULL && pBandCount != NULL)
//    {
//        int iMaxBandIndex = pBandIndex[0];    //找?出?最?大?的?波?段?索?引?序?号?
//        for (int i=1; i<*pBandCount; i++)
//        {
//            if (iMaxBandIndex < pBandIndex[i])
//                iMaxBandIndex = pBandIndex[i];
//        }
//
//        if(iMaxBandIndex > iBandCount)
//        {
//          
//
//            GDALClose((GDALDatasetH) pDSrc);
//            return RE_PARAMERROR;
//        }
//        
//        iNewBandCount = *pBandCount;
//    }
//
//    LT_Envelope enExtent;
//    enExtent.setToNull();
//
//    if (pExtent == NULL)    //全?图?计?算?
//    {
//        double dPrj[4] = {0};    //x1,x2,y1,y2
//        ImageRowCol2Projection(dGeoTrans, 0, 0, dPrj[0], dPrj[2]);
//        ImageRowCol2Projection(dGeoTrans, pDSrc->GetRasterXSize(), pDSrc->GetRasterYSize(), dPrj[1], dPrj[3]);
//        enExtent.init(dPrj[0], dPrj[1], dPrj[2], dPrj[3]);
//
//        pExtent = &enExtent;
//    }
//    
//    dGeoTrans[0] = pExtent->getMinX();
//    dGeoTrans[3] = pExtent->getMaxY();
//    dGeoTrans[1] = dGeoTrans[1] / fResX;
//    dGeoTrans[5] = dGeoTrans[5] / fResY;
//
//    int iNewWidth  = static_cast<int>( (pExtent->getMaxX() - pExtent->getMinX() / ABS(dGeoTrans[1]) + 0.5) );
//    int iNewHeight = static_cast<int>( (pExtent->getMaxX() - pExtent->getMinX() / ABS(dGeoTrans[5]) + 0.5) );
//
//    GDALDataset *pDDst = pDriver->Create(pszOutFile, iNewWidth, iNewHeight, iNewBandCount, dataType, NULL);
//    if (pDDst == NULL)
//    {
//        if(pProgress != NULL)
//            pProgress->SetProgressTip("创?建?输?出?文?件?失?败?！?");
//
//        GDALClose((GDALDatasetH) pDSrc);
//        return RE_FILENOTEXIST;
//    }
//
//    pDDst->SetProjection(strWkt.c_str());
//    pDDst->SetGeoTransform(dGeoTrans);
//
//    GDALResampleAlg eResample = (GDALResampleAlg) nResampleMode;
//
//    if(pProgress != NULL)
//    {
//        pProgress->SetProgressTip("正?在?执?行?重?采?样?...");
//        pProgress->SetProgressTotalStep(iNewBandCount*iNewHeight);
//    }
//
//    int *pSrcBand = NULL;
//    int *pDstBand = NULL;
//    int iBandSize = 0;
//    if (pBandIndex != NULL && pBandCount != NULL)
//    {
//        iBandSize = *pBandCount;
//        pSrcBand = new int[iBandSize];
//        pDstBand = new int[iBandSize];
//
//        for (int i=0; i<iBandSize; i++)
//        {
//            pSrcBand[i] = pBandIndex[i];
//            pDstBand[i] = i+1;
//        }
//    }
//    else
//    {
//        iBandSize = iBandCount;
//        pSrcBand = new int[iBandSize];
//        pDstBand = new int[iBandSize];
//
//        for (int i=0; i<iBandSize; i++)
//        {
//            pSrcBand[i] = i+1;
//            pDstBand[i] = i+1;
//        }
//    }
//    
//    void *hTransformArg = NULL, *hGenImgPrjArg = NULL;
//    hTransformArg = hGenImgPrjArg = GDALCreateGenImgProjTransformer2((GDALDatasetH) pDSrc, (GDALDatasetH) pDDst, NULL);
//    if (hTransformArg == NULL)
//    {
//        if(pProgress != NULL)
//            pProgress->SetProgressTip("转?换?参?数?错?误?！?");
//
//        GDALClose((GDALDatasetH) pDSrc);
//        GDALClose((GDALDatasetH) pDDst);
//        return RE_PARAMERROR;
//    }
//    
//    GDALTransformerFunc pFnTransformer = GDALGenImgProjTransform;
//    GDALWarpOptions *psWo = GDALCreateWarpOptions();
//
//    psWo->papszWarpOptions = CSLDuplicate(NULL);
//    psWo->eWorkingDataType = dataType;
//    psWo->eResampleAlg = eResample;
//
//    psWo->hSrcDS = (GDALDatasetH) pDSrc;
//    psWo->hDstDS = (GDALDatasetH) pDDst;
//
//    psWo->pfnTransformer = pFnTransformer;
//    psWo->pTransformerArg = hTransformArg;
//
//    psWo->pfnProgress = GDALProgress;
//    psWo->pProgressArg = pProgress;
//
//    psWo->nBandCount = iNewBandCount;
//    psWo->panSrcBands = (int *) CPLMalloc(iNewBandCount*sizeof(int));
//    psWo->panDstBands = (int *) CPLMalloc(iNewBandCount*sizeof(int));
//    for (int i=0; i<iNewBandCount; i++)
//    {
//        psWo->panSrcBands[i] = pSrcBand[i];
//        psWo->panDstBands[i] = pDstBand[i];
//    }
//
//    RELEASE(pSrcBand);
//    RELEASE(pDstBand);
//
//    GDALWarpOperation oWo;
//    if (oWo.Initialize(psWo) != CE_None)
//    {
//        if(pProgress != NULL)
//            pProgress->SetProgressTip("转?换?参?数?错?误?！?");
//
//        GDALClose((GDALDatasetH) pDSrc);
//        GDALClose((GDALDatasetH) pDDst);
//
//        return RE_PARAMERROR;
//    }
//
//    oWo.ChunkAndWarpImage(0, 0, iNewWidth, iNewHeight);
//    
//    GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
//    GDALDestroyWarpOptions( psWo );
//    GDALClose((GDALDatasetH) pDSrc);
//    GDALClose((GDALDatasetH) pDDst);
//
//    if(pProgress != NULL)
//        pProgress->SetProgressTip("重?采?样?完?成?！?");
//
//    return RE_SUCCESS;
//}



主函数：jhjz200.cpp


// jhjz200.cpp : 定义控制台应用程序的入口点。
//


// jhjz.cpp : 定义控制台应用程序的入口点。
//


#include "stdafx.h"
#include "IMGALG_H.h"
#include "ImageGeoLocWarp.h"


int _tmain(int argc, _TCHAR* argv[])
{
return 0;
}


#include "gdal_priv.h"
#include "cpl_conv.h" //for CPLMalloc()


#include <string>
#include <vector>
using namespace std;
#include "gdal_priv.h"
#include"ImageGeoLocWarp.h"


 
int main()
{
            





//int iRev = RE_SUCCESS;
对HDF数据中的一个字数据进行处理，下面分别是子数据集的路径和输出文件路径
const char* pszSrc = "HDF4_EOS:EOS_SWATH:\"F:\\GDAL\\data1\\AQUA_20160912_140059_MOD021KM.hdf\":MODIS_SWATH_Type_L1B:EV_1KM_RefSB";
//const char* pszSrc = "F:\\GDAL\\data1\\AQUA_20160912_140059_MOD021KM.hdf";
//const char* pszDst = "D:\\Data\\hdf\\H1BCLR101106020418636_12.tif";


//ImageGeoLocWarp warp(pszSrc, pszDst, "HFA");
//iRev = warp.DoGeoLocWarp(0.0, 0.0, GRA_NearestNeighbour);


//if(iRev == RE_SUCCESS)
// printf("success!\n");
//else
// printf("failed!\n");


//
//::system("pause");
//return 0;
getchar();//不让程序一闪而过
}


参考文献：

http://blog.csdn.net/liminlu0314/article/details/8629298


http://blog.csdn.net/liminlu0314/article/details/7072007

ArcGIS林业制图技术
关注
专栏目录
MODIS几何校正（IDL源代码）
11-30
用IDL开发，实现了读取MODIS头文件信息、几何校正，非常好用。
MODIS图像几何校正程序
11-03
MODIS图像几何校正程序，用IDL语言实现
MRTSwath用于modis L1b数据几何校正
10-29
modis工具，对于处理021KM数据，整合mod03数据
Radon变换图像几何校正(Matlab)
07-02
基于雷登变换的图像几何校正代码(Matlab)
基于数据矢量特性的MODIS数据几何校正 (2013年)
05-17
MODIS探测器的大扫描角多探元并扫方式是造成扫描...结果表明：应用该种方法处理的MODIS数据几何校正精度高，方便进行亚像元级精确几何校正，而且可以避免2次或多次像素重定位导致信息丢失，不会出现“BowTie”现象。
MODIS数据几何校正算法设计及其IDL实现
07-18
### MODIS数据几何校正算法设计及其IDL实现 #### 概述 MODIS（Moderate Resolution Imaging Spectroradiometer）是搭载在美国EOS（Earth Observing System）系列卫星上的一个重要传感器，以其独特的中等空间分辨率...
MODIS Swath数据的几何校正-Python批处理
kymgoal的博客
01-03 3330
目录前言MODIS Swath数据Python实现后记 前言 我是学遥感的，自己的主力编程工具其实一直是IDL。IDL这语言很小众，听过的人不多，但我确实也拿他做了不少事情，科研和教学都得用。目测未来很长一段时间也脱离不了它，毕竟自己用了快10年，那么多代码，短时间内用其他语言全部重新实现一遍，也不太现实。这几年机器学习很火，自己其实也一直有兴趣，然而能抽出时间去学习是件很难的事情，毕竟不能光看理论而不去动手实践。要实践，恐怕最适合的语言就是Python了。Python我之前并没用过，我觉得想尽快熟悉这门语
MODIS 几何校正操作流程.docx
03-27
总结来说，MODIS几何校正和正射校正是遥感图像处理的重要步骤，涉及到数据的预处理、控制点的选取、校正参数的设定以及校正效果的验证。通过这些步骤，可以提高图像的空间精度，使其更好地服务于土地覆盖分类、生态...
IDL编写的MODIS图像几何校正程序
07-31
至于标签“idl modis几何校正程序”，它强调了程序使用的编程语言和主要功能，便于搜索和分类。用户可以通过这些标签找到相关的资源和讨论。 在压缩包中，只有一个文件名为“IDL编写的MODIS图像几何校正程序”，这...
使用GDAL对HDF数据进行校正
weixin_34218579的博客
02-28 503
在使用HDF的数据时，首先要对其进行校正处理。比如MODIS数据中，有个去除蝴蝶结现象。一般出现这种情况会在低分辨率卫星数据中（个人理解），比如气象卫星，海洋卫星等，这类数据一般的数据分辨率都很低，几百米到几公里，主要用于大范围的检测。 下面以MODIS数据为例，简单说明一下蝴蝶结现象出现的原因。 MODIS图像在其扫描线宽度方向由扫描条带组成，条带宽度为10（1000米分辨率）、20（50...
ENVI：如何对自带GLT表的图像进行几何校正？
炒茄子的博客
08-14 1462
1.1 打开ENVI软件1.2 打开MODIS数据1.3 打开工具箱寻找 Georeference By Sensor工具进行几何校正1.3.1 工具所在位置1.3.2 打开 MODIS 几何校正工具并操作校正的参数设置上面的图片操作中选择投影类型不能随便选择，需要查看当前img文件的元数据了解到图像的投影类型（这里是UTM横轴墨卡托投影），基准面是WGS-84椭球体，带号可以计算也可以直接看软件正中最下方的显示（这里显示是）。2.1 打开ASAR数据。...
ENVI实战—一文学会使用传感器自带信息配准工具进行几何校正
最新发布
Promising_GEO的博客
04-16 766
②影像校正：如图3所示，打开“几何校正”工具箱，打开“传感器几何校正”，选择“ASAR几何校正”，利用ENVI自带的ASAR信息配准工具展开校正，选择导入的影像文件，选择“Geographic Lat/Lon”，设置相关输出参数，完成校正。①读取影像：打开文件，点击“打开为”，选择“光学传感器”，选择“EOS”，点击“MODIS”，选择.hdf文件打开，此时MODIS影像可加载至窗口中。③设置输出参数：将投影方式修改为“UTM”，对输出文件进行命名，点击确定，完成MODIS数据的几何校正。
MODIS 数据几何校正
阳光在手指间
07-11 5181
1、可以使用ENVI的Georeference MODIS 工具进行校正 2、使用ENVI 的Build GLT 进行校正 需要输入的X、Y分别为MODIS数据中的Longitude和Latitude数据 Build GLT 做好后，再使用Georeference from GLT工具 3、使用MCTK，可从ENVI APP SOTRE下载：http://www.enviidl.c...
MODIS数据几何校正（IDL）
Lianj.Qin
01-06 4812
使用IDL对MODIS L1B数据进行几何校正。
Python GDAL批处理MODIS产品
江北20190411的博客
05-07 2926
批量几何校正MODIS产品 from osgeo import gdal, osr import os import time ''' 批量读取文件夹中hdf格式的modis产品数据并进行几何校正(WGS84) ''' def readHdfWithGeo(hdfFloder, saveFloder): #获取输入文件夹中所有的文件名 hdfNameList = os.listdi...
python gdal 重采样_用GDAL对影像重采样的一些要点
weixin_39709674的博客
12-17 434
1、RasterIO只能使用最临近插值法做影像金字塔时，开始用的GDALDataSet的RasterIO通过控制行列宽度来抽层，代码写起来很简洁，奇怪一直没见设置重采样参数的地方，急于实现功能，没有怎么深究，基本功能全都做好后发现分层质量很差。文档里边没找到资料介绍这些，查了很久在Maillist上发现RasterIO() level resampling should all be "neare...
【ArcGIS】GEE下载的MODIS数据存在坐标偏移的问题
博客
09-10 2520
使用Google Earth Engine 以及地区矢量数据进行批量下载、拼接、重投影MODIS数据后，用arcgis加载显示，发现发现从GEE上导出的MODIS数据在本地加载上矢量图后发现存在坐标偏移的问题 查阅资料了解到这是GEE的坐标加密保护，终于找到了方法。 问题出在从GEE导出的数据其投影信息的参数和原本MODIS官网中下载的数据给的投影信息存在差异。 解决步骤如下： 1.在nasa官网下载MODIS影像 hdf格式查不到投影信息 试了用MRT转成TIF也无法保留原来的投影信息 需要自己在env
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值