精华内容
下载资源
问答
  • ENVI MODIS多光谱几何校正操作流程,包括几何精校正及正射校正详细流程。MODIS初级产品数据都是未经过任何预处理的数据,图幅较大,空间分辨率底,无法目视手动选择控制点,本身的特点决定了比较适用于进行自动几何...
  • 用IDL开发,实现了读取MODIS头文件信息、几何校正,非常好用。
  • modis swath tool 几何校正的安装说明 可以从网站上直接查看 http://www.cnblogs.com/gisbingxin/archive/2009/04/22/1441550.html
  • MODIS图像几何校正程序,用IDL语言实现
  • MODIS数据几何校正(IDL)

    千次阅读 2020-01-06 15:09:17
    使用IDL对MODIS L1B数据进行几何校正


    使用IDL对MODIS L1B数据进行几何校正。


    1 下载mctk.sav文件

    mctk工具非常便利地将原始HDF文件数据集分条目在列表中显示,打开hdf文件后,选择相应的数据集和投影方式,即可生成具有坐标投影信息的MODIS数据文件。点击mctk.sav工具下载mctk.sav文件。


    2 配置mctk.sav文件

    将mctk.sav拷贝至…\Exelis\ENVI51\extensions文件夹下面,运行ENVI5.1,点击主界面的Extensions文件夹,如果显示mctk则配置成功。

    在这里插入图片描述


    3 使用IDL进行校正

    IDL代码如下:

    Level 1B example
    forward_function envi_proj_create
    forward_function mctk_create_bridges
    pro test_batch_modis_conversion_level_1b
      compile_opt IDL2
      
      t1=systime(1) ;获取程序开始运行的系统时间
      print, 'begin...'
      
      modis_l1b_file = 'E:\MOD03\MOD021KM.A2017244.0735.006.2017244214951.hdf'
      
      output_location = 'E:\MOD03-out\'
    
      output_rootname = 'level_1b'
      
      ;Calibration method schema is:
      ;0 = Radiance / Emissivity, 1 = Reflectance / Emissivity, 2 = Radiance / Brightness Temp
      calib_method = 1
    
      ;Output method schema is:
      ;0 = Standard, 1 = Projected, 2 = Standard and Projected
      out_method = 1
      
      output_projection = envi_proj_create(/geographic)
      ;output_projection = envi_proj_create(/geographic, datum='WGS-84')
      
      ;do not put the bridge creation/destruction code inside a loop
      bridges = mctk_create_bridges()
    
      ;Choosing linear interpolation
      interpolation_method = 0
    
      convert_modis_data, in_file = modis_l1b_file, out_path = output_location, $
        out_root = output_rootname, out_method = out_method, $
        interp_method = interpolation_method, $
        out_proj = output_projection, calib_method = calib_method, $
        sd_pos = [0, 1, 2, 3], /no_msg, background = 0.0, r_fid_array = r_fid_array, $
        r_fname_array = r_fname_array, bridges = bridges, msg = msg  
          
      mctk_destroy_bridges, bridges
      
      t2=systime(1)
      print,'耗时(min):',(t2-t1)/60.0
    
    end
    

    最后结果如下,左图为校正前的影像,右图为校正后的影像。

    在这里插入图片描述
    欢迎大家批评指正。

    展开全文
  • IDL编写的MODIS图像几何校正程序,很实用的的程序,大家试试看,不错的
  • MODIS数据几何校正算法设计及其IDL实现
  • MODIS 数据几何校正

    千次阅读 2018-07-11 19:27:37
    1、可以使用ENVI的Georeference MODIS 工具进行校正 2、使用ENVI 的Build GLT 进行校正 需要输入的X、Y分别为MODIS数据中的Longitude和Latitude数据 Build GLT 做好后,再使用Georeference from GLT工具 3、...

    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.com/appstore/

    4、使用MRT进行处理。

     

     

     

    展开全文
  • IDL编写的MODIS图像几何校正程序,方便大家对MODIS数据进行校正
  • MODIS探测器的大扫描角多探元并扫方式是造成扫描...结果表明:应用该种方法处理的MODIS数据几何校正精度高,方便进行亚像元级精确几何校正,而且可以避免2次或多次像素重定位导致信息丢失,不会出现“BowTie”现象。
  • 批量几何校正MODIS 3km AOD产品
  • modis数据的几何校正

    千次阅读 2017-02-14 17:07:33
    * Describe: 地理查找表校正,用于HDF数据的校正算法 * ****************************************************************************/ #define IMAGEGEOLOCWARP_H class ImageGeoLocWarp { public: ...

    头文件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

    展开全文
  • modis工具,对于处理021KM数据,整合mod03数据
  • 如何用python做modis的大气校正,有没有比较好的包,之前在网上看见有用py6s做的,但是安装和教程都很少,安装步骤里面的文件也找不到地方下载了
  • 目录前言MODIS Swath数据Python实现后记 前言 我是学遥感的,自己的主力编程工具其实一直是IDL。IDL这语言很小众,听过的人不多,但我确实也拿他做了不少事情,科研和教学都得用。目测未来很长一段时间也脱离不了它...

    前言

    我是学遥感的,自己的主力编程工具其实一直是IDL。IDL这语言很小众,听过的人不多,但我确实也拿他做了不少事情,科研和教学都得用。目测未来很长一段时间也脱离不了它,毕竟自己用了快10年,那么多代码,短时间内用其他语言全部重新实现一遍,也不太现实。这几年机器学习很火,自己其实也一直有兴趣,然而能抽出时间去学习是件很难的事情,毕竟不能光看理论而不去动手实践。要实践,恐怕最适合的语言就是Python了。Python我之前并没用过,我觉得想尽快熟悉这门语言的套路,可能就是重写一些之前写过的简单程序,这篇以及以后的博文,其实就是记录一下我学习的过程。顺便可能还能帮助别人解决一些遥感学科内的一些共性技术问题,以后教学上可能也能用,何乐不为。但我写的代码可能一点也不符合软件行业的规范,我只管实现功能,毕竟我也不是专门学计算机的人,大佬轻拍。第一篇文章,就写一写MODIS Swath数据的几何校正和批处理吧,这个可能是很多遥感专业的初学者常遇到的问题之一。

    MODIS Swath数据

    MODIS的产品数据,如果按存放的形式来讲,我认为可以分为两类,第一类是Swath数据,第二类是Grid数据。Swath数据实际上就是按传感器的实际扫描过程,一行一行地存放的,每一景的空间范围并不固定。MODIS的1级产品和大部分的2级大气产品,都属于Swath数据。Grid数据其实是在Swath数据的基础上做了后处理,在全球分块的基础上,把原有的数据全部归算到了固定空间范围的格网内,且每个像元代表的面积也是相等的。Grid数据常用的投影格式就是正弦投影,MODIS绝大部分的地表产品就是这一形式。由于MODIS传感器自身的特性,在扫描的过程中星下点像元的空间分辨率和边缘区域像元的空间分辨率有较大的差异,越往边缘,像元代表的空间范围就越大(如下图)。但在数据存放的过程中,是按像元挨像元的方式,最终出来的结果,和实际的地表范围相比,看起来就像是变形了。要还原出Swath数据代表的真实空间范围,就必须要对它做几何校正(Georeference)。

    在这里插入图片描述

    Python实现

    整体来说,MODIS Swath数据的几何校正,需要它数据里提供的经纬度数据集或者专门的MOD/MYD03数据来实现。简单地说,经纬度数据集里记录了每个像元的实际经纬度,根据这个信息再结合一些最近邻采样的办法,便可以把数据实际空间分辨率不统一的问题初步解决。这个过程其实我之前用IDL的时候手动实现过,具体过程和ENVI里的GLT工具有些类似,但我并没有调用ENVI的GLT接口。然而,既然是用Python,当然就有更简单的办法。
    这里用到的是GDAL库,然而这就是我遇到的第一个困难。GDAL库太庞大了,目前来说我只是大概知道一些我需要的东西,接口和关键字太多了,很难全搞清楚。我这里用到了其中的两个API:Translate和Warp。总体来说,Translate在这里面的作用是构建一个用于几何校正的文件,这个和ENVI里的GLT有点像,但又不完全一样。如果有兴趣,可以把我代码段里的os.remove(vrt_file)这句注释掉,这样在数据目录下生成的vrt文件就可以保留,用文本编辑器打开能看到里面都包含了些啥信息,这里就不详细讲了。Warp的作用则是根据vrt里的信息,对目标数据集进行几何校正,中间的一些关键字和输出结果有关,后面再讲。
    这个程序里只定义了一个函数,swath_georeference。输入参数是:
    in_file:输入swath文件名。
    geo_file:输入地理数据文件名。这部分功能其实我没有实现完,目前只是针对文件里本身就带了Longitude和Latitude的Swath文件,所以可以看到下面的in_file和geo_file我输入的同一个。后面再完善单独输入MOD/MYD03的部分。
    target_dataset_name:待几何校正的目标数据集名称。这个可以用HDF Explorer、Panoply之类的工具查看,也可以在for subdataset in subdatasets:这句代码下添加print(subdataset)看到每个数据集的信息。
    out_file:输出文件名。
    out_geo_range:输出的空间范围。如果不指定,就用None,代表输出的是数据实际包含的全部空间范围。
    x_res_out:输出的像元分辨率(x方向)。这里用的是经纬度单位,我用的是的MYD04_3K的气溶胶产品,所以设置了一个0.03°,和3 km大概等效的一个分辨率。
    y_res_out:输出的像元分辨率(y方向)。这里用的是经纬度单位,同上。
    这个代码实际跑的时候,其实只需要对应更改main函数下面的那几行就行了,即input_directory(输入路径)、output_directory(输出路径)、out_file(输出文件名,如果需要的话)、dataset_name(目标数据集名称)、geo_range(输出的经纬度范围,不指定则=None)、x_res(x方向的空间分辨率,°)、y_res(y方向的空间分辨率,°)。
    输入参数正确之后,便可开始对input_directory下所有文件的批量处理。

    import os
    import gdal
    
    
    def swath_georeference(in_file, geo_file, target_dataset_name, out_file, out_geo_range, x_res_out, y_res_out):
        dataset = gdal.Open(in_file)
        subdatasets = dataset.GetSubDatasets()
        dataset_pos = 0
        target_pos = -1
        for subdataset in subdatasets:
            subdataset_name = subdataset[0]
            if subdataset_name.endswith(target_dataset_name):
                target_pos = dataset_pos
            dataset_pos += 1
        if target_pos != -1:
            target_dataset = dataset.GetSubDatasets()[target_pos][0]
            vrt_file = os.path.splitext(in_file)[0] + '.vrt'
            vrt_out = gdal.Translate(vrt_file, target_dataset,
                                     format='vrt', unscale='true')
            georeferenced_data = gdal.Warp(out_file, vrt_out, multithread=True, outputBounds=out_geo_range,
                                           format='GTiff', geoloc=True, dstSRS="EPSG:4326",
                                           xRes=x_res_out, yRes=y_res_out, dstNodata=0.0, outputType=gdal.GDT_Float32)
            print('The georeference of {} is complete.'.format(in_file))
            os.remove(vrt_file)
        else:
            print('{} has no the target dataset.'.format(in_file))
    
    
    def main():
        file_prefix = '.hdf'
        input_directory = 'O:/coarse_data/chapter_1/MODIS_2018_mod04_3k/'
        output_directory = 'O:/coarse_data/chapter_1/MODIS_2018_mod04_3k/geo_out/'
        if not os.path.exists(output_directory):
            os.mkdir(output_directory)
        file_list = os.listdir(input_directory)
        for file_i in file_list:
            if file_i.endswith(file_prefix):
                file_name = input_directory + file_i
                out_file = output_directory + os.path.splitext(file_i)[0] + '_geo.tiff'
                dataset_name = 'Image_Optical_Depth_Land_And_Ocean'
                geo_range = None  # [90.0, 55.0, 120.0, 25.0]
                x_res = 0.03
                y_res = 0.03
                swath_georeference(file_name, file_name, dataset_name, out_file, geo_range, x_res, y_res)
    
    
    if __name__ == '__main__':
        main()
    
    

    后记

    有些地方我也没搞清楚,比如想输出投影格式如UTM这种怎么实现,好像并不是简单的修改dstSRS="EPSG:4326"这一项。目前这段代码只是大概能用,还需要进一步修改完善,毕竟这就是我入坑Python一天之后写出来的,有些东西得慢慢搞清楚,GDAL的API帮助也不太友好。目前只测试了MYD04_3K这个数据,没试过能否直接对1级数据正确处理,但理论上是可以的。Warp这个API里有一个multithread关键字,看说明是控制多线程处理的,但好像我理解有问题,开不开没啥变化,我的CPU使用率毫无波澜。
    关于geo_range的设置,我个人认为如果后续是对特定研究区的处理分析,那这里可以直接设定一个经纬度范围,如注释里的[90.0, 55.0, 120.0, 25.0],90.0和55.0代表左上角点的经纬度坐标,120.0和25.0代表右下角点的经纬度坐标,这样几何校正之后的结果都会固定在这一空间范围内,后续要编程做一些均值处理要方便些。如果后续是要对一天的结果做拼接,那可能设置成None更合适。
    相比ENVI里的GLT接口,GDAL做几何校正要快多了。如果对IDL调用GLT接口感兴趣,可以看:
    如何使用idl做GLT重投影(调用envi接口)
    后面我也看情况增加一些IDL做遥感数据处理的内容,我对IDL远比Python熟悉得多。毕竟能看到这篇文章的,多半都是圈内人,IDL还是有用武之地。

    展开全文
  • 有关遥感影像几何校正问题

    千次阅读 2020-03-05 15:37:19
    几何校正 正射校正 几何配准的区分 几何校正分为不同级别,正射校正可以说是几何校正的最高级别。我们一般所说的几何校正是消除因大气传输、传感器本身、地球曲率等因素造成的几何畸变,主要纠正或者赋予影像平面...
  • Python 利用GDAL对图像进行几何校正

    千次阅读 多人点赞 2018-07-16 20:09:35
    Python+GDAL 对遥感图像进行几何校正 Python以简单高效易上手著称,拥有强大的高级数据结构以及众多第三方包,真是快速编程的绝佳首选。 最近想要利用GDAL库对遥感图像进行几何校正,在网上搜了搜,大部分是...
  • ENVI遥感影像处理及几何校正

    万次阅读 2017-11-21 19:35:00
    一、图像几何校正的概述 1、几何校正方法: 1)利用卫星自带的地理定位文件进行几何校正。主菜单>>>Map>>Georeference传感器的名称,来启动这种矫正方法。 2)Image to Image几何校正。一幅图像没有经过...
  • NASA出品的modis几何校正软件 MODIS 1B数据需处理的过程(几何纠正,DN转反射率,拼接)
  • 2.基于自带定位信息的几何校正 2.1 内容介绍 图像的几何形变一般分为两大类:系统性和非系统性。系统性几何形变一般是由传感器本身引起的,有规律可循、具有可预测性,可以用传感器模型来校正,卫星地面接收站已经...
  • 【ENVI入门系列】03.基于自带定位信息的几何校正   (2014-09-22 17:33:59) 转载▼   分类: ENVI 版权声明:本教程涉及到的...2.1MODIS数据几何校正 2.2ASAR数据几何校正 2.3
  • (23条消息)Python 利用GDAL对图像进行几何校正_python_可乐Mentos的博客-CSDN博客 https://blog.csdn.net/qq_27045589/article/details/81062586 代码参考: python调用GDAL实现几何校正 - 楚彦 - 博客园 ...
  • [ENVI] 遥感图像几何校正及示例

    万次阅读 2017-11-06 17:20:15
    ENVI中几何校正:基于自带定位信息的几何校正、图像精校正(控制点选择)、自动影像配准。
  • ---恢复内容开始--- 一、图像几何校正的概述 1、几何校正方法: ...一幅图像没有经过几何校正的删个文件或者已经经过几何校正的栅格文件作为基准图,通过两幅图上选择同名点来配准另一幅栅格文件,使相...
  • ENVI4.7几何校正bug

    2010-05-08 10:49:00
    最近吧envi升级到了4.7,结果处理modis结合校正的时候总是提示出现错误,更可恶的是错误之后envi会把原有的modis的hdf原始 数据删除掉。经查证这是envi4.7的bug(好大的bug),现已有补丁,处理方法如下文: 原文:...
  • 我看网上都说mod03数据地理定位文件,但是没找到怎么使用mod03数据对modis其他数据进行几何校正,想请问有没有大佬做过的,求指导,万分感谢~~~

空空如也

空空如也

1 2 3 4 5 ... 12
收藏数 221
精华内容 88
关键字:

modis几何校正