IDL学习笔记（四）MODIS数据处理。MODIS数据介绍，以及Swath数据处理

MODIS传感器介绍

MODlS (Moderate Resolution lmaging Spectroradiometer 中分辨率成像光谱仪 ) :
卫星平台 : Terra（上午11点左右过境）、Aqua（下午1点半左右过境
空间分辨率 : 250 m、500 m、1 km(at nadir，星下点分辨率)
波段数 : 36 个
过境频次 : 2 - 4 次 / 天
扫描方式 : 穿轨扫描(cross-track-scanning
幅宽 : 2230 km

modis成像方式：穿轨扫描！一行一行拼起来的。由于是星下点分辨率，星下点分辨率最高，两侧可能分辨率会低一些，会有形变以及分辨率的降低。即：影像边缘空间分辨率远大于星下点。

MODIS 数据产品

Swath产品的范围和 level 1 级别数据范围完全类似
Grid 产品，范围规划到固定的范围上。

注意：

两大类产品：
1.使用SWATH产品必须进行重投影(地理参考)
2.GRID产品自带正弦投影，但一般也需要重投影成其他格式

Swath 数据

扫描后结果是 ：南北颠倒，以及有变形。越往边缘，像元越粗，空间越被压缩。这样的数据，需要进行重投影，这样才能使用：

Grid 数据

Grid正弦投影产品全球格网划分(等面积分块）。


类似于“h28v05”是指的是地理位置模块位置是哪一个。grid数据进行了一些升尺度降尺度的操作，已经使得每个像元分辨率都是一个分辨率了（例如：500m)。即 Grid产品自带投影，且每个像元分辨率一致。但为了方便显示、操作、画图，往往需要重投影。

MODIS Swath 数据重投影

MODIS Swath 数据重投影方法：

构建地理查找表，实现重投影。使用envi接口进行实现。
需要：

1~6：构建GLT
7~8：重投影

1.x方向坐标（经度）数据集
2.x方向（纬度）数据集
3.方向数据集坐标类型。一般都是选下图

4.目标输出格式选择
5.输出的像元分辨率
6.图像转置角度，一般设置为 0 

虽然输出的也是经纬度坐标，但是像元排列已经改变了。中纬度地区，一度大概110公里左右。

7.选择GLT文件（想对哪个文件做重投影，就可以选哪一个了）
8.待投影文件

如果有负号，就说明原始影像上没有这个经纬度，但是又最相近的一个，检索后赋值过来的（最近邻查找）。这样的重投影的方法是一个近似的处理，逼近于真实。（包含误差）

对应ENVI接口

1~6：构建GLT：

ENVI_GLT_DOIT

7~8：重投影：

ENVI_GROREF_FROM_GLT_DOIT

具体步骤:
1读取经纬度数据集、目标数据集
2 输出以上数据集
3 调用ENVI接口读取数据集(enviopen_file)
4调用ENVI接口制作GLT(envi_proj_create,envi_glt_doit)
5 调用ENVI接口从GLT对目标数据集作重投影(envi_georef from_glt_doit)
6 删除中间数据，保留目标结果

function hdf4_data_get,file_name,sds_name
  sd_id=hdf_sd_start(file_name,/read)
  sds_index=hdf_sd_nametoindex(sd_id,sds_name)
  sds_id=hdf_sd_select(sd_id,sds_index)
  hdf_sd_getdata,sds_id,data
  hdf_sd_endaccess,sds_id
  hdf_sd_end,sd_id
  return,data
end

function hdf4_attdata_get,file_name,sds_name,att_name
  sd_id=hdf_sd_start(file_name,/read)
  sds_index=hdf_sd_nametoindex(sd_id,sds_name)
  sds_id=hdf_sd_select(sd_id,sds_index)
  att_index=hdf_sd_attrfind(sds_id,att_name)
  hdf_sd_attrinfo,sds_id,att_index,data=att_data
  hdf_sd_endaccess,sds_id
  hdf_sd_end,sd_id
  return,att_data
end

pro modis_swath_glt
  compile_opt idl2
  envi,/restore_base_save_files
  envi_batch_init
  
  input_directory='E:\Data\IDL\chapter_3\modis_swath\'
  output_directory='E:\Data\IDL\chapter_3\modis_swath\geo_out\'
  directory_exist=file_test(output_directory,/directory)
  if (directory_exist eq 0) then begin
    file_mkdir,output_directory
  endif

  file_list=file_search(input_directory,'*.hdf')
  file_n=n_elements(file_list)

  for file_i=0,file_n-1 do begin
    start_time=systime(1)
    modis_lon_data=hdf4_data_get(file_list[file_i],'Longitude')
    modis_lat_data=hdf4_data_get(file_list[file_i],'Latitude')
    modis_target_data=hdf4_data_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean')
    scale_factor=hdf4_attdata_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean','scale_factor')
    fill_value=hdf4_attdata_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean','_FillValue')

    modis_target_data=(modis_target_data ne fill_value[0])*modis_target_data*scale_factor[0]

    out_lon=output_directory+'lon_out.tiff'
    out_lat=output_directory+'lat_out.tiff'
    out_target=output_directory+'target.tiff'
    write_tiff,out_lon,modis_lon_data,/float
    write_tiff,out_lat,modis_lat_data,/float
    write_tiff,out_target,modis_target_data,/float

    envi_open_file,out_lon,r_fid=lon_fid;打开经度数据，获取经度文件id
    envi_open_file,out_lat,r_fid=lat_fid;打开纬度数据，获取纬度文件id
    envi_open_file,out_target,r_fid=target_fid;打开目标数据，获取目标文件id

    out_name_glt=output_directory+file_basename(file_list[file_i],'.hdf')+'_glt.img'
    out_name_glt_hdr=output_directory+file_basename(file_list[file_i],'.hdf')+'_glt.hdr'
    in_proj=envi_proj_create(/geographic)
    out_proj=envi_proj_create(/geographic)
    envi_glt_doit,$ ; 下面的语言，为了好看，才输入输出分行，$是分行的意思
      i_proj=in_proj,x_fid=lon_fid,y_fid=lat_fid,x_pos=0,y_pos=0,$;指定创建GLT所需输入数据信息
      o_proj=out_proj,pixel_size=0.03,rotation=0.0,out_name=out_name_glt,r_fid=glt_fid;指定输出GLT文件信息

    out_name_geo=output_directory+file_basename(file_list[file_i],'.hdf')+'_georef.img'
    out_name_geo_hdr=output_directory+file_basename(file_list[file_i],'.hdf')+'_georef.hdr'
    envi_georef_from_glt_doit,$
      glt_fid=glt_fid,$;指定重投影所需GLT文件信息。 fid 的意思是 file的id
      fid=target_fid,pos=0,$;指定待投影数据id。如果目标需要多个层进行重投影，可以pos = [0,1,2]
      out_name=out_name_geo,background = 0,r_fid=geo_fid;指定输出重投影文件信息,geo_fid是重投影后的id,background默认是0

    envi_file_mng,id=x_fid,/remove  ;envi_file_mng关闭文件（不是删除）
    envi_file_mng,id=y_fid,/remove
    envi_file_mng,id=target_fid,/remove
    envi_file_mng,id=glt_fid,/remove
    envi_file_mng,id=geo_fid,/remove
    file_delete,[out_lon,out_lat,out_target,out_name_glt,out_name_glt_hdr] ;删除中间数据

    end_time=systime(1)
    print,'The program is end, this '+ string(end_time-start_time) + 's.'
  endfor
  envi_batch_exit,/no_confirm
end

调用ENVI接口，以下三句话是套话：

  compile_opt idl2
  envi,/restore_base_save_files
  envi_batch_init

1. compile_opt idl2 ：告诉程序准备调用接口
2. envi,/restore_base_save_files ：调用语句
3. envi_batch_init &