• 【转】Envi调用MODIS Reprojection Tool(MRT)对MODIS产品进行批处理拼接、重投影、裁切


    1 熟悉MRT

     

    MODIS产品的类型不同,一景HDF格式的影像包含的波段也各不相同。MRT处理时需要选择处理波段,0表示不作处理,1表示处理,首先要确定影像的波段数。

     

     

    1、拼接

    Mrtmosaic.exe程序用来拼接影像。调用方式为:

    mrtmosaic -i "g:\n%1.txt" -s " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 " -o "g:\temp.hdf"

    命令行帮助如下

    Usage: mrtmosaic -i input_filenames_file -t -h -o output_filename

                     -s spectral_subset "b1 b2 ... bN"

                     -g filename for the log file

    需要指定三个参数:

    -I 指定输入文件,可以采取两种方式

    11.hdf 2.hdf 3.hdf

    2、把影像的完整路径保存到txt文件中,作为input参数。如-I MOD092008001.txt

    -s 指定需要处理的波段,同样可以采用两种方式

    1、直接给出,如 –s “0 0 0 0 0 0 0 0 0 0 0 0”,注意英文引号;

    2、指定一个txt路径,让程序读取;

    -o 指定输出路径,一般直接给出

    1、直接给出,如 –o g:\tmp.hdf。注意直接存为HDF格式,便于后续处理。

    2、重投影、裁切

    Resample.exe用来重投影、裁切是MRT程序的核心。调用方式为:

    命令行帮助如下:

    Usage: resample -p parameter_file [options]

     

    Options that override parameter file specifications:

       -i input_file_name

       -o output_file_name

       -r resampling_type [NN BI CC NONE]

       -t projection_type [AEA ER GEO HAM IGH ISIN LA LCC MERCAT MOL PS SIN TM UTM]

       -j projection_parameter_list "p1 p2 ... p15"

       -s spectral_subset "b1 b2 ... bN"

       If using the -s switch, the SDSs should be represented as an

       array of 0s and 1s. A '1' specifies to process that SDS;

       '0' specifies to skip that SDS. Unspecified SDSs will not be processed.

       If the -s switch is not specified, then all SDSs will be processed.

       -a spatial_subset_type [INPUT_LAT_LONG INPUT_LINE_SAMPLE OUTPUT_PROJ_COORDS]

       -l spatial_subset "ULlat ULlong LRlat LRlong"

                   -or- "ULline ULsample LRline LRsample (0-based)"

                   -or- "ULprojx ULprojy LRprojx LRprojy"

          NOTE: line/sample must be specified for the             highest resolution

     of all SDSs specified             to be processed in the product.

       -u UTM_zone

       -x pixel_size

       -g filename for the log file

    可以只指定1个参数:

    –p 读入prm参数文件,进行处理。例如resample -p "g:\prrmMOD092008001.prm"prm文件如下:

    INPUT_FILENAME = g:\tmp_%1.hdf输入文件

    SPECTRAL_SUBSET = ( 1 1 1 1) 前面mrtmosaic拼接结果有4个波段

     

    SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG 经纬度裁切方式

     

    SPATIAL_SUBSET_UL_CORNER = ( 33.0 108.0 ) 左上纬经度

    SPATIAL_SUBSET_LR_CORNER = ( 28.0 117.0 ) 右下纬经度

     

    OUTPUT_FILENAME = F:\MRT_out\myd%1.tif 输出路径,不同波段自动区分

     

    RESAMPLING_TYPE = NEAREST_NEIGHBOR 最近邻采样方法

     

    OUTPUT_PROJECTION_TYPE = UTM 输出文件投影方式utm

     

    OUTPUT_PROJECTION_PARAMETERS = (投影参数

     0.0 0.0 0.0

     0.0 0.0 0.0

     0.0 0.0 0.0

     0.0 0.0 0.0

     0.0 0.0 0.0 )

     

    DATUM = WGS84 大地水准面

     

    UTM_ZONE = 50投影分带带号

    如果不会设置,可以先在图形界面里设置一次,把设置的结果保存下来。

    打开刚刚保存的参数文件

    去掉以#号打头的注释,文件显示如下:

    INPUT_FILENAME = G:\TmpMosaic.hdf

     

    SPECTRAL_SUBSET = ( 1 1 1 1 1 1 1 1 1 1 1 1 )

     

     

    SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG

     

    SPATIAL_SUBSET_UL_CORNER = ( 32.0 107.0 )

    SPATIAL_SUBSET_LR_CORNER = ( 29.0 116.0 )

     

    OUTPUT_FILENAME = G:\mod092001001.tif

     

    RESAMPLING_TYPE = NEAREST_NEIGHBOR

     

    OUTPUT_PROJECTION_TYPE = UTM

     

    OUTPUT_PROJECTION_PARAMETERS = (

     0.0 0.0 0.0

     0.0 0.0 0.0

     0.0 0.0 0.0

     0.0 0.0 0.0

     0.0 0.0 0.0 )

     

    DATUM = WGS84

     

    UTM_ZONE = 50

    如果研究区的经纬度范围不同,可以修改

    SPATIAL_SUBSET_UL_CORNER = ( 32.0 107.0 )

    SPATIAL_SUBSET_LR_CORNER = ( 29.0 116.0 )两行

    如果需要不同的投影,如经纬度投影,可以修改

    OUTPUT_PROJECTION_TYPE = UTM

     

    OUTPUT_PROJECTION_PARAMETERS = (

     0.0 0.0 0.0

     0.0 0.0 0.0

     0.0 0.0 0.0

     0.0 0.0 0.0

     0.0 0.0 0.0 )

     

    DATUM = WGS84

     

    UTM_ZONE = 50

    如果拼接时只选择了1个波段,可以修改

    SPECTRAL_SUBSET = ( 1 )

    当然这句

    INPUT_FILENAME = G:\TmpMosaic.hdf批处理是必须改掉的了

    关于像素分辨率,一般留空,也就是默认不改变影像的像素大小。

    2关于DOS批处理裁切MOD产品

    单次处理

    Dir g:\mod11*A%1*.hdf /s /b >> g:\n%1.txt

    f:\mrt\modis\bin\mrtmosaic -i "g:\n%1.txt" -s " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 " -o "g:\temp_%1.hdf"

    f:\mrt\modis\bin\resample -p "g:\prm%1.prm"

    del g:\prm%1.prm /q /f

    del g:\n%1.txt /q /f

    del g:\temp_%1.hdf /q /f

    批量处理

    for /f %%i in (e:\lis.txt) do @单次处理文件名 %%i

    如果双击运行,用%%i号表示变量。在命令行中%i表示变量。

    其中e:\lis.txt是年份及天数

    2008001

    2008009

    2008017

    2008025

    2008033

    2008041

    2008049

    G:\n%1.txt在运行时为g:\n2008001.txt,其中包括2008年第1天的多景MODIS产品影像名。

    Dir g:\mod11*A%1*.hdf /s /b >> g:\n%1.txt

    用来列出g盘(数据存储盘)以mod11打头的指定日期的的hdf格式文件,即MODIS影像。

    这里要说明一下,MODIS产品命名的规则

    MYD11A2.A2002185.h27v05.005.2007221162948.hdf

    从左至右,前7位表示产品类型MYD11指下午星(MYD)11号产品(land surface temperature),A2表示处理级别。A2002185表示过境时间是2002年第185天,h27v05MODIS产品在全球的轨道行列号,中国一般是h27v05h27v06h28v05h28v06005表示HDF5文件格式,早期MODIS采用HDF4格式存储,所以有时也可见到004的产品。2007等一串数字表示影像处理的年月日时分秒,hdf表示存储格式为EOS-HDF

    了解了MODIS产品命名规则之后,可以用DOS dir命令来挑选符合条件的MODIS产品路径了。

    >> g:\n%1.txt 表示把结果转存为g:\n%1.txt

    %1表示接收的第一个参数,即2008001

    3Envi生成时间列表

    2008001

    2008009

    2008017

    2008025

    这样的列表文件可以用在Envi中定义一个pro来生成:

    pro daylist,startdate,enddate,result_txt_path

            ;example daylist,2008001,2008013,'example.txt'

            ;

            ;

            a=(startdate mod 1000) /8 *8+(startdate /1000)*1000+1

           

            b=((enddate mod 1000) /8 +1)*8+(enddate /1000)*1000+1

            help,a,b       

            step=8

            openw,lun,result_txt_path,/get_lun

            res=a;     

            for i=a,b,step do begin

               

                printf,lun,string(res,format='(%"%7d")')

               

                res=res+step

            endfor

            free_lun,lun

            print,'done';

            close,/all

           

     

     

    end

    在命令行窗口敲入

    daylist,2009001,2009033,'f:\g4.txt'

    打开f:\g4.txt如下:

    2009001

    2009009

    2009017

    2009025

    2009033

    2009041

    4Envi平台下整合前面工作

    留给读者思考

     提示:spawn函数

  • 相关阅读:
    SQL Server的Linked Servers
    Pycharm新建Python项目
    C#-使用Newtonsoft.Json实现json字符串与object对象互转
    C#-使用HttpListener创建http服务
    Win10系统使用Gitblit搭建局域网Git服务器
    C# HttpClient类库
    C#读写自定义的多字段配置文件
    Postman软件-请求HTTP接口
    SoapUI软件-测试Web Service接口
    Python开发 必备
  • 原文地址:https://www.cnblogs.com/HomeGIS/p/1628353.html
Copyright © 2020-2023  润新知