1. 课堂内容
批量重投影
这里是重投影的整体思路：
1. 获取Modis Swath数据(这里只获取Lat、Lon、Aod(气溶胶厚度)三个数据集以及aod数据集的两个属性)，并对aod数据进行简单的处理
2. 调用二次开发接口以Lat、Lon数据作为原材料创建Glt文件
3. 调用二次开发接口以Glt文件作为参照，对Aod数据进行重投影并输出
这里主要是二次开发接口比较陌生，其实就是调用包我感觉。
另外对于二次接口的使用比较陌生可以使用ENVI先进行重投影的操作然后将ENVI的重投影操作和IDL的代码进行对比的话，你会发现代码和ENVI的重投影操作有很多相似之处，这将会有助于你理解代码的操作和运行
2. ENVI实操与IDL编程
function get_hdf4_ds, file_path, ds_name ; 构建函数, 用于得到hdf4文件中数据集的数据 ; 获取文件id file_id = hdf_sd_start(file_path, /read) ; 获取数据集index ds_index = hdf_sd_nametoindex(file_id, ds_name) ; 获取数据集的id ds_id = hdf_sd_select(file_id, ds_index) ; 传入数据集的id以及数据集的名称 ; 获取数据集的数据 hdf_sd_getdata, ds_id, ds_data ; 传入数据集id，在传入一个变量用于接函数获取的该数据集的数据 ; 关闭已经打开数据集以及文件 hdf_sd_endaccess, ds_id hdf_sd_end, file_id ; 返回值 return, ds_data end function get_hdf4_att_info, file_path, ds_path, att_name ; 构建函数: 用于获取hdf4文件的数据集的相关属性 ; 获取文件id file_id = hdf_sd_start(file_path, /read) ; 获取数据集的index ds_index = hdf_sd_nametoindex(file_id, ds_path) ; 获取数据集的id ds_id = hdf_sd_select(file_id, ds_index) ; 获取该数据集下的属性的index att_index = hdf_sd_attrfind(ds_id, att_name) ; 获取该数据集的某一属性内容 hdf_sd_attrinfo, ds_id, att_index, data=att_data ; 关闭数据集以及文件 hdf_sd_endaccess, ds_id hdf_sd_end, file_id ; 返回值 return, att_data end pro week_four_study1 ; 对modis产品(HDF4)进行重投影(两步走, 创建glt文件，借glt文件投影) ; 记录一下开始时间 start = systime(1) ; 由于这里需要使用envi的api接口，所以需要声明 compile_opt idl2 envi,/restore_base_save_files envi_batch_init ; 首先获取modis文件中的lon、lat(前面两个数据是用于创建glt即地理信息查找表)、aod数据,并存储在硬盘中，方便后续envi的api去调用文件 ; 路径 in_path = 'D:/IDL_program/experiment_data/chapter_3/modis_swath' out_path = 'D:/IDL_program/experiment_data/chapter_3/modis_swath/geo_out' ; 假如out_path目录不存在，那么创建 if file_test(out_path, /directory) eq 0 then begin file_mkdir, out_path endif ; 获取in_path中的所有modis数据(即获取hdf4文件) ; 获取所有hdf4文件的路径(以array形式返回) file_path_array = file_search(in_path, '*.hdf') ; 传入目录，传入需要查找的hdf4文件的后缀即.hdf ; 获取该目录下的所有hdf4文件的个数 file_count = n_elements(file_path_array) ; 获取lon、lat、aod数据集以及aod数据集的几个属性 for file_i = 0, file_count - 1 do begin ; 该循环下的文件的路径 file_path = file_path_array[file_i] ; 获取数据集数据 lon_data = get_hdf4_ds(file_path, 'Longitude') lat_data = get_hdf4_ds(file_path, 'Latitude') aod_data = get_hdf4_ds(file_path, 'Image_Optical_Depth_Land_And_Ocean') ; 获取aod数据集的属性(_FillValue, scale_factor) aod_fv = get_hdf4_att_info(file_path, 'Image_Optical_Depth_Land_And_Ocean', '_FillValue') aod_sf = get_hdf4_att_info(file_path, 'Image_Optical_Depth_Land_And_Ocean', 'scale_factor') ; 需要注意，这里返回的属性值是数组形式array(即便它只有一个float数),所以下面使用需要注意 ; 对aod数据进行处理 aod_data = (aod_data ne aod_fv[0]) * aod_fv[0] * aod_sf[0] ; 注意aod_fv和aod_sf均是数组 ; 输出路径(也许你会发现，每次循环会把之前的目录覆盖，但是不要担心，我们在本循环中就会使用到现在存储的数据) lon_path = out_path + '/lon_out.tiff' lat_path = out_path + '/lat_out.tiff' aod_path = out_path + '/aod_out.tiff' ; 存储三个数据集的数据 write_tiff, lon_path, lon_data, /float write_tiff, lat_path, lat_data, /float write_tiff, aod_path, aod_data, /float ; 将存储的文件给envi的api调用(或者说我们告诉api他需要的文件在哪里) envi_open_file, lon_path, r_fid=lon_id envi_open_file, lat_path, r_fid=lat_id envi_open_file, aod_path, r_fid=aod_id ; 传入存储数据集的路径，由关键字传参(r_fid=)得到该数据集(或者说是存储的tiff文件)的id(用变量lon_id等接收) ; 输出glt文件的路径 out_path_glt = out_path + '/glt.img' out_path_glt_hdr = out_path + '/glt.hdr' ; 没有这个行不行，等会儿试一下 ; 输入输出投影信息 in_prj = envi_proj_create(/GEOGRAPHIC) ; 输入是经纬度 out_prj = envi_proj_create(/GEOGRAPHIC) ; 输出也是经纬度 ; 创建glt文件 envi_glt_doit, $ ; $ 表示换行 i_proj=in_prj, x_fid=lon_id, y_fid=lat_id, x_pos=0, y_pos=0, $ ; 创建glt需要传入的相关信息 o_proj=out_prj, pixel_size=pixel_size, rotation=0.0, $ ; 输出glt需要传入的相关信息 out_name=out_path_glt, r_fid=glt_id ; 这个也是，输出的路径以及返回的创建好的glt文件的id ; 输出重投影结果的路径 out_georef_path = out_path + '/' + file_basename(file_path, '.hdf') + '_georef.img' out_georef_hdr_path = out_path + '/' + file_basename(file_path, '.hdf') + '_georef.hdr' ; 使用glt文件进行重投影 envi_georef_from_glt_doit, $ glt_fid=glt_id, pos=0, $ ; 传入glt文件的id以及图层数(由于只有一层填0即可) fid=aod_id, $ ; 传入待投影的aod文件的id out_name=out_georef_path, r_fid=georef_id ; 传入已经投影好的文件所存放的路径，以及返回该文件的id ; 关闭已经打开的文件 envi_file_mng, id=lon_id, /remove envi_file_mng, id=lat_id, /remove envi_file_mng, id=aod_id, /remove envi_file_mng, id=glt_id, /remove envi_file_mng, id=georef_id, /remove ; 重投影中间产生了一些文件需要删除,因为我们只需要重投影文件，至于经纬度文件、glt文件如果你想保留那么请不要删除 file_delete, [lon_path, lat_path, aod_path, out_path_glt, out_path_glt_hdr] ; 记录一下结束时间 stop = systime(1) endfor print, stop - start end
编译运行：
