modis数据预处理

MODIS数据预处理

1. 波段设置

Modis影像有三种打开方式,一般我们用打开外部文件的方式打开科学数据集,因为需要数

据集中的一些辅助信息(主要是太阳几何,卫星几何).但是这样打开之后显示的波段从1开始的,而数据集中对应的modis通道并不是这个顺序.通过菜单栏中的

basic_tools->preprocessing->data_specific utilities->view HDF dataset attributes 可以打开数据集里每个要素的属性表,在里面选中需要的HDF文件中的数据集，就会打开其属性表，波段内容如下

Attribute 3-5: \

\

对应打开的HDF文件里1KM辐亮度文件的波段数,一共16个波段.其中13/14波段比较特殊,都有hi和lo两组数据,它们是传感器高敏感度和低敏感度两种状态下获取到的DN值,分别对应于较暗地物和较亮地物,使用哪个文件根据需要而定.但是在太湖湖区,13/14波段大部分区域效果都不太好.值会很大,出现溢出.可能是由于太湖的高浑浊度.

2. 几何校正

几何校正有三种方法:

1) 用envi自带模块进行几何校正,通过菜单栏中的

Map->Georeferences MODIS

选中envi中已经打开的需要校正的数据集，输入研究区的地理位置，如下图左，投影用UTM，基准面用WGS-84，区域根据经纬度确定。输入完成，envi会自动校正，并执行去蝴蝶结效应算法，有点是能对我们需要的那些波段进行校正。缺点也很明显。如下图右，校正结束的图像会失去原始图像四个角的信息，这样就无法和GLT校正的图像很好的匹配起来，不利于一些后续的处理。

2) 用GLT，即是查找表法对图像进行几何校正

Map->Georeference from input Geometry->buid GLT 用来建立查找表。在弹出的对话框中选择查找表的XY信息，其中X对应图像经度信息，Y对应纬度信息。然后只需要规定投影、基准面和区位信息，就可以生成一个查找表文件。这个查找表文件的实质也是两幅图像，分别在每个像元上保存着经纬度值，并且像元位置是拉伸到我们规定的输出投影上面去了，而且是逐像元的拉伸。那么剩下的矫正工作就只是把想要矫正的信息和查找表一一匹配起来，因此速度也很快。

Map->Georeference from input Geometry->Georeference from GLT

就是上面所说的，把查找表运用到需要矫正的每一个像元上。这里只需要选择想要矫正的影响就好了。有一点需要注意：这里查找表的分辨率比辐亮度信息要低，因此如果直接用这个查找表对原始的辐亮度图像进行矫正，只能得到原始图像一小部分的信息。因为矫正时是按照像元行列号一一对应查询的。解决方案有两个：可以把原始图像分辨率降低，重采样成和查找表一样大小的图像，再用查找表进行几何校正，这样运行速度快但是损失原始图像信息；也可以将查找表重采样成和原始的辐亮度图像一样大小的图像，再用这个新的查找表对图像进行矫正。这样就不会损失原始图像信息，但是计算速度会大大降低。视应用选择。 3) IDL批处理

;forward_functionenvi_proj_create PRO Modis_gef_batch

envi, /restore_base_save_files ;恢复ENVI sav文件 envi_batch_init, log_file=’batch.txt’ ;开始批处理模式 inpath = DIALOG_PICKFILE(/DIRECTORY, $ TITLE=\

CD,inpath

filename = FILE_SEARCH('*.HDF') ;print,result

n = N_ELEMENTS(filename)

outpath = DIALOG_PICKFILE(/DIRECTORY, $ TITLE=\ FOR i=0,n-1 DO BEGIN in_name=inpath+filename[i] out_name ='ReGeo'+filename[i]

;设置校正方法

;0 = Radiance \\ Emissivity, 1 = Reflectance \\ Emissivity calib_method = 1

;设置输出方法

;0 = Standard, 1 = Projected, 2 = Standard and Projected out_method = 1

;设置输出投影

output_projection = envi_proj_create(/geographic) ;在输出时设置去除蝴蝶效应

convert_modis_data, in_file=in_name, out_path=outpath, $ out_root=out_name, /l1b, out_method=out_method, $

out_proj=output_projection, calib_method=calib_method, /bowtie, $ sd_pos=[1,3], /no_msg, background=0.0 ENDFOR

envi_batch_exit END

只能生成envi标准格式打开的一些文件，无法对具体的所有辐亮度信息进行矫正。

3. 辐射定标

辐射定标信息也通过打开HDF数据集属性信息查看

basic_tools->preprocessing->data_specific utilities->view HDF dataset attributes

Attribute 3-7: \

316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849, 316.97219849

Attribute 3-8: \

\

Attribute 3-9: \

0.00002352, 0.00001261, 0.00000720, 0.00000538, 0.00000462, 0.00000227, 0.00000168, 0.00000255, 0.00000141, 0.00000232, 0.00000292, 0.00002401, 0.00003249, 0.00002493, 0.00002598

Offsets和scales规定了辐射定标的系数。Modis辐射定标公式是 R = scales*(DN-offsets) 其中R代表辐射定标之后的辐亮度，scales和offsets分别代表列表中的定标系数，DN代表图像原始数码值。 4. 大气校正

6s大气校正，最关键的就是获取卫星、太阳几何信息。Envi操作稍显繁琐，大致分为以下步骤：

1) 用上面的批处理程序对modis进行几何校正（modis处理工具包也可以），关键是获

取到没有边角信息损失的几何校正影像。因为后面要和几何信息进行叠加；

2) 用GLT方法对太阳几何、卫星几何信息进行波段合成（要用另存为envi格式影像。

因为没有进行过几何校正，图像不包含经纬度信息，layer stacking时会报错，只能save as成一幅没有经纬度信息的多波段图像，这一步是为了减少几何校正步骤）、几何校正；

3) 对几何校正之后的几何信息（四个波段的图像）进行重采样，得到和1)矫正得到的

modis校准影像等行等列的图像；

4) 对3)、4)得到的结果进行波段叠加，就可以在湖心区，用查看z_profile的方式，查

看该点的太阳几何、卫星几何信息。

为了简化该步骤，可以用IDL辅助查找这些信息。大致思路和envi操作相似，不过在查找湖心像元几何信息的部分做了一些改进，直接用经纬度进行匹配，可以减少一些步骤。（envi也可以这么操作）

function GLT_regeo , filename ;, pos1,pos2,pos3,pos4 ;获取经纬度信息和需要进行几何校正的波段 ;输出文件6个波段信息 ;0:longtitude ;1:latitude

;2:sensor_zenith ;3:sensor_azimuth ;4:solar_zenith ;5:solar_azimuth

outname = strsplit(filename,'.',/extract)

ENVI_OPEN_DATA_FILE, filename, r_fid=fid_x, /hdf_sd, $ hdfsd_dataset=1, hdfsd_interleave=0

ENVI_OPEN_DATA_FILE, filename, r_fid=fid_y, /hdf_sd, $ hdfsd_dataset=0, hdfsd_interleave=0 ;

ENVI_OPEN_DATA_FILE, filename, r_fid=fid_sensor_zenith, /hdf_sd, $ hdfsd_dataset=13, hdfsd_interleave=0

ENVI_OPEN_DATA_FILE, filename, r_fid=fid_sensor_azimuth, /hdf_sd, $ hdfsd_dataset=14, hdfsd_interleave=0

ENVI_OPEN_DATA_FILE, filename, r_fid=fid_solar_zenith, /hdf_sd, $ hdfsd_dataset=16, hdfsd_interleave=0

ENVI_OPEN_DATA_FILE, filename, r_fid=fid_solar_azimuth, /hdf_sd, $ hdfsd_dataset=17, hdfsd_interleave=0

;查询图像信息

envi_file_query, fid_sensor_zenith, dims=dims

;将四个有用信息的图像合成一个图像，便于后续的几何校正工作 envi_doit, 'cf_doit', $

fid=[fid_sensor_zenith,fid_sensor_azimuth,fid_solar_zenith,fid_solar_azimuth], $ pos=[0,0,0,0], dims=dims, $ remove=0, $

; out_name=out_name, $ r_fid=info_fid ,/in_memory

;将四个单独的几何信息文件释放

envi_file_mng, id=fid_sensor_zenith, /remove envi_file_mng, id=fid_sensor_azimuth, /remove envi_file_mng, id=fid_solar_zenith, /remove envi_file_mng, id=fid_solar_azimuth, /remove ;建立查找表

i_proj = envi_proj_create(/geographic, datum = 'WGS-84') o_proj = envi_proj_create(/utm, datum = 'WGS-84', zone=50) envi_glt_doit, i_proj=i_proj, $ o_proj=o_proj, $

rotation = 0, $ ;保证图像北面朝上 r_fid=glt_fid, $

x_fid=fid_x, y_fid=fid_y, $

out_name = outname[0] + '_GLT',$ x_pos=0, y_pos=0; ,/in_memory

;根据查找表对数据进行几何校正

envi_doit, 'envi_georef_from_glt_doit', fid=info_fid, $

glt_fid=glt_fid, out_name=outname[0] + '_regeo', pos=[0,1,2,3], $ subset=dims, r_fid=regeo_fid envi_file_query,regeo_fid,dims = dims envi_file_mng, id=glt_fid, /remove

envi_file_mng, id=info_fid, /remove envi_file_mng, id=fid_x, /remove envi_file_mng, id=fid_y, /remove

;将最后生成数据的id号，和最后图像的维度信息，存成一个数组返回 print,regeo_fid

return,[regeo_fid,dims] end ;

pro get_modis_info

file_name = DIALOG_PICKFILE(TITLE='select a HDF file') if strlen(file_name) eq 0 then return info = GLT_regeo(file_name)

;info[0]:函数输出图像的id值

;info[1:5]:函数输出图像的纬度信息 ;

;用记录下的太湖图像中心部分的经纬度，转化成投影坐标，再转化成图像坐标，对返回的图像做行列号匹配，即是最佳位置。

;地理坐标转化为投影坐标 center_lat = 31.22723889 center_lon = 120.2057722 datum = 'WGS-84'

iproj = ENVI_PROJ_CREATE(/geographic)

oproj = ENVI_PROJ_CREATE(/utm,zone=50,datum = datum) ENVI_CONVERT_PROJECTION_COORDINATES, $ center_lon, center_lat, iproj, $ oxmap, oymap, oproj

;投影坐标转化为图片坐标

envi_convert_file_coordinates, info[0], xf, yf, oxmap, oymap print,xf,yf x = fix(xf) y = fix(yf) print,x,y

sensor_zenith = ENVI_GET_DATA(fid=info[0], dims=info[1:5], pos=0) print,info[0]

help,sensor_zenith

sensor_azimuth = ENVI_GET_DATA(fid=info[0], dims=info[1:5], pos=1) solar_zenith = ENVI_GET_DATA(fid=info[0], dims=info[1:5], pos=2) solar_azimuth = ENVI_GET_DATA(fid=info[0], dims=info[1:5], pos=3) print,sensor_zenith[x,y] print,sensor_azimuth[x,y] print,solar_zenith[x,y] print,solar_azimuth[x,y] end

之后就可以利用输出的几何信息输入6s进行大气校正，至此预处理结束。

envi_file_mng, id=info_fid, /remove envi_file_mng, id=fid_x, /remove envi_file_mng, id=fid_y, /remove

;将最后生成数据的id号，和最后图像的维度信息，存成一个数组返回 print,regeo_fid

return,[regeo_fid,dims] end ;

pro get_modis_info

file_name = DIALOG_PICKFILE(TITLE='select a HDF file') if strlen(file_name) eq 0 then return info = GLT_regeo(file_name)

;info[0]:函数输出图像的id值

;info[1:5]:函数输出图像的纬度信息 ;

;用记录下的太湖图像中心部分的经纬度，转化成投影坐标，再转化成图像坐标，对返回的图像做行列号匹配，即是最佳位置。

;地理坐标转化为投影坐标 center_lat = 31.22723889 center_lon = 120.2057722 datum = 'WGS-84'

iproj = ENVI_PROJ_CREATE(/geographic)

oproj = ENVI_PROJ_CREATE(/utm,zone=50,datum = datum) ENVI_CONVERT_PROJECTION_COORDINATES, $ center_lon, center_lat, iproj, $ oxmap, oymap, oproj

;投影坐标转化为图片坐标

envi_convert_file_coordinates, info[0], xf, yf, oxmap, oymap print,xf,yf x = fix(xf) y = fix(yf) print,x,y

sensor_zenith = ENVI_GET_DATA(fid=info[0], dims=info[1:5], pos=0) print,info[0]

help,sensor_zenith

sensor_azimuth = ENVI_GET_DATA(fid=info[0], dims=info[1:5], pos=1) solar_zenith = ENVI_GET_DATA(fid=info[0], dims=info[1:5], pos=2) solar_azimuth = ENVI_GET_DATA(fid=info[0], dims=info[1:5], pos=3) print,sensor_zenith[x,y] print,sensor_azimuth[x,y] print,solar_zenith[x,y] print,solar_azimuth[x,y] end

之后就可以利用输出的几何信息输入6s进行大气校正，至此预处理结束。

本文来源：https://www.bwwdw.com/article/d6et.html