1. 实验内容
1.1 提取日期、经纬度、AOD(气溶胶厚度)
提取/ coarse_data /chapter_1/MODIS_2018_mod04_3k/ 目录下所有 MODIS 气溶胶产品中 Image_Optical_Depth_Land_And_Ocean 数据集内与北京经纬度（ 39.90°N ， 116.40°E ）最接近的点的 AOD 有效结果，并按每行内容为日期号 经度 纬度 AOD 格式输出到 IDL 控制台，如：
• 2018143 106.39740 38.243156 0.21600001
• .
• .
• .
• .
• 2018147 113.02930 37.370262 0.21000001
1.2 输出形式逗号分隔
将题目1中的输出形式改为用逗号分隔，如：
• 2013278,116.49,39.91,2.767
1.3 输出日期时转化为年-月-日形式
将题目2中的日期号，转换为年-月-日的形式输出，同时数值保留4位小数，如：
• 2013-01-09,116.49,39.91,2.767
2. 知识储备
2.1 本程序会使用的一些函数的介绍
3. 编程
另外，这里因为实验内容中2和3的内容基本相似，而且实验3的内容相对更难，所以这里不演示实验2的操作，相信大家都可以自己写出，加油！（！！！！！！浩楠哥哥好帅！！！！！！！）
3.1 编程内容
; 构建函数——提取文件中的气溶胶厚度数据集(ROD)的数据 function get_ds_data, file_path, file_ds_name ; 获取文件file_name以及文件内部file_ds_name的id ; 获取文件id file_id = hdf_sd_start(file_path, /read) ; 获取数据集的index file_ds_index = hdf_sd_nametoindex(file_id, file_ds_name) ; 获取数据集的id file_ds_id = hdf_sd_select(file_id, file_ds_index) ; 获取数据集的data hdf_sd_getdata, file_ds_id, file_ds_data ; 关闭数据集以及文件 hdf_sd_endaccess, file_ds_id hdf_sd_end, file_id ; 返回值————数据集的数据 return, file_ds_data ; 这里注意return后面有一个半角的逗号 end ; 构建函数————提取文件中的数据集的内部属性的数据 function get_att_data, file_path, file_ds_name, file_ds_att_name ; 获取文件的id file_id = hdf_sd_start(file_path, /read) ; 获取数据集的index file_ds_index = hdf_sd_nametoindex(file_id, file_ds_name) ; 获取数据集的id file_ds_id = hdf_sd_select(file_id, file_ds_index) ; 获取数据集下的属性file_ds_att_name的index file_ds_att_index = hdf_sd_attrfind(file_ds_id, file_ds_att_name) ; 获取指定属性的数据 hdf_sd_attrinfo, file_ds_id, file_ds_att_index, data=file_ds_att_data ; 关闭数据集以及文件————这是习惯 hdf_sd_endaccess, file_ds_id hdf_sd_end, file_id ; 返回值————属性 return, file_ds_att_data ; 这里注意return后面有一个半角的逗号 end ; 前面是函数构建，现在是程序的开始 pro week_two_test ; 本程序旨在解决提取某点的AOD... ; 所有文件所在文件夹路径 folder_path = 'D:\IDL编程设计\实验数据\chapter_1\MODIS_2018_mod04_3k' ; 获取文件夹下所有文件的————》通过file_search()函数获取文件夹下所有hdf文件的名称 file_path_list = file_search(folder_path, '*.hdf', count=file_count) ; 返回str数组(另外注意这里可以返回获取的文件数量， 使用count=————等价于下面的n_elements()) ; 上面函数中第一个参数传入文件夹的路径, 第二个参数传入 需要指定获取文件的类型(譬如这里是*.hdf)。 注意，这里搜索的文件不只是本文件夹，还有本文件夹的子文件等等。 ; 待提取点的经纬度 extract_lat = 39.90 extract_lon = 116.40 ; 寻找距离上面待提取点最近的点的AOD——思路——使用距离公式求得最近的点以及位置（经纬度和行列数）, 输出该点的经纬度以及AOD ; 首先每个文件都有要获取的数据，所以需要循环 ; 获取循环所需的文件的个数————》通过n_elements()函数获取数组的元素个数(其实这里多余，只是为了展示既可以像上面count=得到文件数，也可以使用函数获取，当然，现在这种方法普遍性更好) file_count = n_elements(file_path_list) ; 传入一个数组，返回数组的元素个数 for file_i = 0, file_count - 1 do begin ; 文件数-1不要忘了 ; 通过文件名获取日期 ; 先通过文件的路径获取文件名称 file_name = file_basename(file_path_list[file_i]) ; 在处理文件名称获取日期————》通过strmid()函数对字符串进行切片处理 file_date = strmid(file_name, 10, 7) ; 获取该日期的儒略日 ; 思路：直接获取需要该日期的年月日，而我们发现得到的file_date是年加该年已经过去的天数，所以我们获取上一年的最后一天的儒略日加上年已经过去的天数就是儒略日 ; 获取该日期的的年、已经过的天数 ; file_date_year = fix(strcompress(strmid(file_name, 10, 4))) ; fix将str型的转化为int型 file_date_year = strmid(file_name, 10, 4) ; fix将str型的转化为int型 ; file_date_days = fix(strcompress(strmid(file_name, 14, 3))) file_date_days = strmid(file_name, 14, 3) ; 获取上一年最后一天的儒略日 julian_lastyear = imsl_datetodays(31, 12, file_date_year - 1) ; 上一年需要减一 ; 获取该年的儒略日 file_date_julian = julian_lastyear + file_date_days ; 获取该年的年月日 imsl_daystodate, file_date_julian, file_day, file_month, file_year ; 获取索引为file_i的文件中AOD数据集的AOD数据 file_ds_aod = get_ds_data(file_path_list[file_i], 'Image_Optical_Depth_Land_And_Ocean') ; 获取索引为file_i的文件中经纬度数据集的经纬度数据 file_ds_lon = get_ds_data(file_path_list[file_i], 'Longitude') file_ds_lat = get_ds_data(file_path_list[file_i], 'Latitude') ; 获取属性值——》这里计算真实的aod数据需要两个属性值（一个scale_factor, 另一个是_filevalue） file_aod_sf = get_att_data(file_path_list[file_i], 'Image_Optical_Depth_Land_And_Ocean', 'scale_factor') file_aod_fv = get_att_data(file_path_list[file_i], 'Image_Optical_Depth_Land_And_Ocean', '_FillValue') ; 计算真实的aod数据 file_ds_aod = (file_ds_aod ne file_aod_fv[0]) * file_ds_aod * file_aod_sf[0] ; 由于获取的属性值是array，所以需要指定索引 ; 计算所有的点与待提取点的距离 distance = ((file_ds_lon - extract_lon) ^ 2 + (file_ds_lat - extract_lat) ^ 2) ^ 0.5 ; 找到最小的距离 min_distance = min(distance) ; 找到最小距离所在的点的行列数(因为是distance是数组， 所以是行列数) min_pos = where(min_distance eq distance) ; 如果aod是0， 那么说明aod是无效的， 因为前面我们对file_ds_aod ne file_aod_fv[0]， 所以，aod=0说明aod是无效值，那么我们不输出 if file_ds_aod[min_pos] eq 0.0 then begin continue endif ; 输出 print, [file_year, file_month, file_day, file_ds_lon[min_pos], file_ds_lat[min_pos], file_ds_aod[min_pos]], format='(i0, "-", i02, "-", i02, ",", 3(f0.3, :, ","))' ; 上面的fromat如果没有:， 那么输出时的末端正常来说会有一个逗号， 但是我们觉得逗号是多余的，所以加上:可以使得最后面的逗号不输出 ; 输出还需要注意f表示浮点型，i表示整型， 需要控制0的输出和C语言是一样的, 譬如f0.3中3表示小数位保留三位,采用四舍五入保留 ; 另外，我们print输出数字时通常数字前面会带有很多空格，那是因为该数字的有效位数有这么多，但是实际上我们输出时通常又不需要这些空格，可以使用f0， 那个0起这个作用 ; 再有， f012.4这些我觉得其实和C语言很相像， 所以这里讲也不一定会听的很明白，我的建议是自己去试试，你会有很对感悟的。 endfor end
3.2 输出结果：
IDL> week_two_test
% Compiled module: GET_DS_DATA.
% Compiled module: GET_ATT_DATA.
% Compiled module: WEEK_TWO_TEST.
2018-05-02,116.379,39.887,0.329
2018-05-03,116.395,39.899,0.408
2018-05-04,116.374,39.877,0.820
2018-05-04,116.354,39.889,0.802
2018-05-06,116.478,39.287,0.375
2018-05-18,116.386,39.903,0.744
2018-05-18,113.432,39.607,0.134
2018-05-22,116.382,39.886,0.594
2018-05-23,116.385,39.892,0.404
2018-05-23,116.300,40.478,0.196
2018-05-25,116.422,39.910,1.019
2018-05-27,116.770,38.303,0.406
2018-05-27,116.417,39.908,0.292
2018-05-28,116.397,39.892,0.435
2018-05-29,116.391,39.890,0.309
2018-05-30,116.405,39.912,0.359
2018-05-31,116.364,39.901,0.566
IDL>
