【转】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 指定输入文件，可以采取两种方式

1、1.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天，h27v05是MODIS产品在全球的轨道行列号，中国一般是h27v05、h27v06、h28v05、h28v06，005表示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函数