1. 实验内容
2. 知识储备
这里就不会对之前学过的一些函数做出声明，只对其中会使用的二次开发接口做出介绍
3. ENVI实操对对应DL代码部分
这一部分实在没有精力，有需要再来吧，累死了。
简单说一下
1. 打开需要重投影的hdf4文件，选择好数据集(这里只需要选择lon、lat、aod三个数据)
2.打开toolbox工具箱，打开build glt应该是我记得——》这是进行创建glt地理查找表的入口
3. 选择输入输出的投影信息、x方向上的数据集、y方向的数据集。。。。。
等等等。。。。真的累了
4. 编程
function get_hdf4_ds, file_path, ds_name ; 函数功能: 获取数据集的数据 ; 获取文件的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) ; 获取数据集的数据 hdf_sd_getdata, ds_id, ds_data ; 关闭数据集以及文件 hdf_sd_endaccess, ds_id hdf_sd_end, file_id ; 返回值 return, ds_data end function get_hdf4_att, file_path, ds_name, att_name ; 函数功能: 获取数据集的某一属性内容 ; 获取文件的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) ; 获取该属性的indedx 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_test ; 批量重投影modis swath产品并输出为tiff格式(且限制范围为39-42, 115-118) ;总体思路: ; 1. 获取lon、lat、aod数据集以及aod数据集的两个属性并进行处理、最后存储在硬盘 ; 2. 上面的数据处理包括对aod的数据进行处理，以及对三个数据进行范围限制处理(因为我们只需要北京地区范围) ; 3. 通过envi二次开发接创建glt文件(需要使用上面存储在硬盘中的lon、lat文件) ; 4. 通过envi二次开发接口参照glt地理信息查找表对aod文件进行重投影 ; 5. 将上面得到的重投影文件转化为geotiff格式 ; 由于这里需要使用envi的api接口，所以需要声明(其实我也不知道这是什么) compile_opt idl2 envi,/restore_base_save_files envi_batch_init ; 路径 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) eq 0 then begin ; 不存在函数返回0，否则返回1 file_mkdir, out_path endif ; 获取in_path中所有的hdf文件的路径 file_path_array = file_search(in_path, '*.hdf') ; 传入目录，指定搜索的文件格式是'*.hdf'——》返回该目录下的所有hdf文件 ; 获取hdf文件的总数量 file_count = n_elements(file_path_array) ; 循环每一个hdf文件，获取其中的数据，并进行一系列的处理最后输出为北京地区的tiff格式 for file_i = 0, file_count - 1 do begin ; 获取该循环下的文件路径 file_path = file_path_array[file_i] ; 获取该循环下的文件名称 file_name = file_basename(file_path, '.hdf') ; 获取该循环下的文件的数据集以及属性数据 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_fv = get_hdf4_att(file_path, 'Image_Optical_Depth_Land_And_Ocean', '_FillValue') aod_sf = get_hdf4_att(file_path, 'Image_Optical_Depth_Land_And_Ocean', 'scale_factor') ; 对aod数据集进行处理 aod_data = (aod_data ne aod_fv[0]) * aod_data * aod_sf[0] ; 你为什么老是忘记这是数组所以需要[] ; 存储时，我们只需要39-42， 115-119的经纬度范围,所以需要对范围进行类似裁剪 ; 获取在目标范围的所有的像元的下标(以数组形式返回) pos = where((lon_data gt 115.0) and (lon_data lt 119.0) and (lat_data gt 39.0) and (lat_data lt 42.0), pos_pixel_count) ; 119.0 39.0这些最好弄一个浮点型而不是整型会更好 ; 像元的个数 ; pos_pixel_count = n_elements(pos) ; 注意，这里本来是既可以n_elements()函数获取pos的元素个数，也可以再where()函数运行时指定/count(譬如现在就是) ; 但是，这里需要注意，如果where里面全为0(也就是说该图像没有一个像元在北京地区，那么where会返回一个数值)，那么pos会返回-1 ; 使用where()函数接收的count参数会返回0，而使用n_elements()函数会返回1，这回导致下面这个if判断出现错误 ; 所以推荐使用where()函数去取count即数组元素个数 ; 如果像元个数小于0，那么也就是说这个文件一点也不包含北京地区，那么我们就不需要对其进行后续的处理了 if pos_pixel_count eq 0 then continue ; 那么结束循环 ; 获取lon_data数据的总的行列数(其它lat_data、aod_data的行列数肯定也是这个) lon_size = size(lon_data) ; 传入数据，返回该数据的的一些信息 ; aod_size是一个数组，第0个数字表示aod_data数组的维度，第1个数字表示aod_data数组的总列数，第2个数表示总行数，第3个数表示元素类型的代号，第4个数表示aod_data数组的所有元素的个数 ; 我们需要总列数 max_cols = lon_size[1] ; 找到最左最右最上最下的行列号 rows = pos / max_cols ; 这个自己列个简单的数组试一下就知道是怎么回事，直接讲比较费劲 cols = pos mod max_cols ; 这个也是 up_row = min(rows) down_row = max(rows) left_col = min(cols) right_col = max(cols) ; 对输出的lon_data、lat_data、aod_data的范围进行限制处理(下面的取值中[a:b, c:d]表示从第a列到第b列，从第c行到第d行) out_lon_data = lon_data[left_col:right_col, up_row:down_row] out_lat_data = lat_data[left_col:right_col, up_row:down_row] out_aod_data = aod_data[left_col:right_col, up_row:down_row] ; 将得到的数据存储在硬盘中，方便后envi二次开发接口api调用这些数据文件 ; 输出路径 out_lon_path = out_path + '/out_lon.tiff' out_lat_path = out_path + '/out_lat.tiff' out_aod_path = out_path + '/out_aod.tiff' ; 输出 write_tiff, out_lon_path, out_lon_data, /float write_tiff, out_lat_path, out_lat_data, /float write_tiff, out_aod_path, out_aod_data, /float ; 使用envi的二次开发接口api打开上面存储的三个数据文文件——》传入文件的路径，使用关键字传参返回该文件的id envi_open_file, out_lon_path, r_fid=lon_id ; 用关键字r_fid返回打开文件的id，这里用变量lon_id接收(下面类似) envi_open_file, out_lat_path, r_fid=lat_id envi_open_file, out_aod_path, r_fid=aod_id ; 创建glt文件(需要使用lon和lat数据文件) ; 创建的glt文件的输出路径 out_glt_path = out_path + '/out_glt.img' out_glt_hdr_path = out_path + '/out_glt.hdr' ; 输入和输出的投影信息 in_proj = envi_proj_create(/GEOGRAPHIC) out_proj = envi_proj_create(/GEOGRAPHIC) ; 调用envi_glt_doit创建glt文件 envi_glt_doit, $ i_proj=in_proj, x_fid=lon_id, y_fid=lat_id, x_pos=0, y_pos=0, $ ; 输入设置 o_proj=out_proj, pixel_size=pixel_size, rotation=0.0, $ ; 输出设置pixel_size=pixel_size表示选择默认 out_name=out_glt_path, r_fid=glt_id ; 输出 ; 使用glt文件对aod数据进行重投影 ; 重投影结果的输出路径 out_georef_path = out_path + '/georef.img' out_georef_hdr_path = out_path + '/georef.hdr' ; 开始进行重投影 ENVI_GEOREF_FROM_GLT_DOIT, $ glt_fid=glt_id, pos=0, $ ; 输入glt文件的id以及所在层数 fid=aod_id, $ ; 输入需要进行重投影的文件的id out_name=out_georef_path, $ ; 传入输出的路径 r_fid=georef_id ; 输出重投影文件的id ; 对上面的重投影转化为tiff格式 ; 获取重投影文件的map_info信息————》因为转化为geotiff需要地理信息 map_info = envi_get_map_info(fid=georef_id) ; 传入文件id，返回该文件的map_info信息(结构体形式) ; 但是map_info不是所有的信息都是有用的(对于我们创建geo_info结构体而言) ; 获取map_info的第一行(索引为0开始)的数据——》返回数组，包含四个数值，前两数是某一像元的行列数(一般是左上角的原点，行列数为0, 0)，后两个数是对应前面像元的经纬度数值 pixel_geo = map_info.(1) ; 取值方式和之前的不一样，需要注意 ; 获取map_info第2行的数据——》一个像元的宽和长分别表示的距离(单位是°)，实际上就是空间分辨率嘛 pixel_space = map_info.(2) ; 获取重投影文件的数据 ; 首先需要获取重投影文件的dims——》调用envi_file_quer接口并传入重投影文件的id，通过关键字传参获取该重投影文件的dims envi_file_query, georef_id, dims=georef_dims ; 再获取重投影文件的数据——》 georef_data = envi_get_data(fid=georef_id, pos=0, dims=georef_dims) ; 创建geo_info结构体 geo_info={$ MODELPIXELSCALETAG:[pixel_space[0],pixel_space[1],0.0],$ MODELTIEPOINTTAG:[0.0,0.0,0.0,pixel_geo[2],pixel_geo[3],0.0],$ ; 这里需要注意:我他妈傻子了，填的是pixel_geo[0],[1]——》错了，这两个数分别是左上角点的列号行号，后面两个数是左上角点的经度和纬度，所以是填pixel_geo[2], pexel_geo[3] GTMODELTYPEGEOKEY:2,$ GTRASTERTYPEGEOKEY:1,$ GEOGRAPHICTYPEGEOKEY:4326,$ GEOGCITATIONGEOKEY:'GCS_WGS_1984',$ GEOGANGULARUNITSGEOKEY:9102,$ GEOGSEMIMAJORAXISGEOKEY:6378137.0,$ GEOGINVFLATTENINGGEOKEY:298.25722} ; 转化为tiff文件的输出路径 out_georef_tiff_path = out_path + '/' + file_name + '_georef.tiff' ; 转化为geotiff文件 write_tiff, out_georef_tiff_path, georef_data, /float, geotiff=geo_info ; 关闭打开的文件 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 file_delete, [out_lon_path, out_lat_path, out_aod_path, out_glt_path, out_glt_hdr_path, $ out_georef_path, out_georef_hdr_path] endfor ; 调用接口结束声明 envi_batch_exit,/no_confirm end
运行结果展示：
产生的geotiff文件：
5. 题外话
5.1 n_element()与一些函数自带的count参数返回的区别
; 注意，这里本来是既可以n_elements()函数获取pos的元素个数，也可以再where()函数运行时指定/count(譬如现在就是)
; 但是，这里需要注意，如果where里面全为0(也就是说该图像没有一个像元在北京地区，那么where会返回一个数值)，那么pos会返回-1
; 使用where()函数接收的count参数会返回0，而使用n_elements()函数会返回1，这回导致下面这个if判断出现错误
; 所以推荐使用where()函数去取count即数组元素个数
5.2 发现一个难以发现的错误
补(代码已经订正>>>2022-7-14)：这里发现了重大错误，虽然原来的代码可以运行，也能出图，但是，出的每一幅图的左上角的经纬度都是0°，0°。这是由于我写Geotiff结构体时，误操作了，将每一幅图的左上角点的列号行号作为左上角点的经度纬度填上去了。
下面是错误代码：
下面是修改之后的代码：
