1. 课后作业内容
2. 知识储备
2.1 相关函数的使用
3. 编程
3.1 编程内容
关于一些代码的注意事项，以及我的一些见解、看法都一一写在了代码内（以注释形式）
; 通过文件名称以及文件内的数据集名称获取数据集数据 function get_ds_data, file_path, ds_name ; 获取文件id file_id = hdf_sd_start(file_path, /read) ; 传入文件路径(file_path)，传入读写方式(/read) ; 获取数据集index ds_index = hdf_sd_nametoindex(file_id, ds_name) ; 传入数据集所在文件的id， 传入数据集的名称 ; 获取数据集的id ds_id = hdf_sd_select(file_id, ds_index) ; 传入数据集所在文件的id， 传入数据集的index ; 获取数据集的数据 hdf_sd_getdata, ds_id, ds_data ; 第一个位置传入数据集的id， 第二个传入变量ds_data用来接收获取的数据 ; 关闭已经打开的数据集以及文件————》这是专业素养！ok？ hdf_sd_endaccess, ds_id ; 传入数据集的id可以准确找到该数据集再内存中的位置然后关闭它 hdf_sd_end, file_id ; 类似地， 传入文件的id可以准确的找到该文件在内存中的位置然后精准关闭 ; 返回数据集的数据 return, ds_data ; 注意return后面有一个半角的逗号 end ; 通过文件路径、数据集名称、属性的名称获属性内容 function get_att_data, 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) ; 获取属性的index att_index = hdf_sd_attrfind(ds_id, att_name) ; 传入该属性所在数据集的id， 传入该属性的名称 ; 获取属性的内容 hdf_sd_attrinfo, ds_id, att_index, data=att_data ; 传入该属性所在数据集的id， 传入该属性的index， 第三个位置data=ds_att_data表示函数将获取的data传给我们的变量ds_att_data ; 关闭已经打开的数据集和文件 hdf_sd_endaccess, ds_id ; 传入数据集的id可以准确找到该数据集再内存中的位置然后关闭它 hdf_sd_end, file_id ; 类似地， 传入文件的id可以准确的找到该文件在内存中的位置然后精准关闭 ; 返回属性内容 return, att_data end function get_date, file_name ; 上面给出了文件的名称, 这里需要从名称里面获取日期 ; 观察其中一个名称可以发现规律_MYD04_3K.A2018121.0545.061.2018121172155.hdf ; 索引10开始往后7个字符(这里是2018121)表示日期，即2018年的第121天 ; 获取原始日期 file_date = strmid(file_name, 10, 7) ; 实际就是切片 ; 传入字符串形式的字符串，再传入开始提取的字符的开始处索引(从0开始), 再传入往后的提取字符的个数 ; 获取原始日期的年、天数(当然，其实我们上面不需要输出file_date,直接再file_name基础上获取即可，你说得对) file_days = strmid(file_date, 4, 3) file_year = strmid(file_date, 0, 4) ; fix()函数将字符串转化为int型(当然，你也可以在上面strmid外边直接套一个fix()) file_days = fix(file_days) ; 这里将其转化为int是为了后面计算（字符串无法计算） file_year = fix(file_year) ; 这里要获取年-月-日的样式，自己写一个程序计算太麻烦(当然，你不介意也可以写) ; 所以可以通过获取2018121的儒略日，通过imsl_daystodate()函数得到年月日 ; 但是我们得到2018121的儒略日需要使用imsl_datetodays()函数，即需要该天的年月日 ; 闭环了不是？不是的，我们可以获取上一年最后一天(这里即是2017-12-31)儒略日啊，再加上121即可 file_lastyear_julian = imsl_datetodays(31, 12, file_year - 1) ; 获取2017-12-31的儒略日， julian表示儒略日 file_year_julian = file_lastyear_julian + file_days ; 获取2018121的儒略日 imsl_daystodate, file_year_julian, file_day, file_month, file_year ; 第一个位置传入儒略日， 后面传入变量用来接收获取的三个值（年、月、日） ; 返回值 return, [file_year, file_month, file_day] end pro week_two_homework ; 所有modis文件所在文件夹的路径 folder_path = 'D:\IDL编程设计\实验数据\chapter_1\MODIS_2018_mod04_3k' ; 获取所有modis文件的文件路径(以array形式输出)————》使用file_search()函数获取 file_path_array = file_search(folder_path, '*.hdf') ; 第一个参数传入你需要检索的路径, 第二个参数传入检索的标准(*.hdf表示只搜索hdf文件) ; 我还是想再强调以下，检索不只是当前文件夹，还有当前文件夹的子文件夹等等 ; 获取所有文件的文件名称(以array形式输出)————》使用file_basename()函数获取 file_name_array = file_basename(file_path_array) ; 这里你既可以传入单个的路径得到该路径指向的文件的名称， 也可以传入存储有多个文件路径的数组(array)得到相同类型的存储有多个文件名称的数组 ; 三个待提取点的相关信息(依次是北京、上海、成都的经纬度) extract_lon_array = [116.40, 121.47, 104.06] extract_lat_array = [39.90, 31.23, 30.67] extract_area_array = ['北京', '上海', '成都'] ; 计算需要提取aod的点的个数(假设我们有若干点， 而且我们不知道是几个点————提升一点难度) extract_area_count = n_elements(extract_area_array) ; 传入数组可以返回该数组内所有元素的总个数 ; 有area_count个待提取点，所以需要循环每一个点 for extract_area_i = 0, extract_area_count - 1 do begin ; 注意，这里索引是从0开始，所以需要减1.如果你不习惯，可以设置extract_area_i=1避免减一(后面数组取值也要相应变化)，但是专业点我建议你还是不要这么去做，会被笑话的(我会笑话) ; 得到目前循环中待提取点的相关信息 extract_lon = extract_lon_array[extract_area_i] extract_lat = extract_lat_array[extract_area_i] ; 创建txt文件，等会儿得到的数据将会放在这个文件里面(由于每一个点的数据单独一个txt文件，所以在外边这一层创建txt文件) openw, 4, './week_two_test/date_lon_lat_aod_' + extract_area_array[extract_area_i] + '.txt' ; 这里说一下，我运行的时候openw1， 1是不可以的，他会说文件已经打开(这里我不是很理解？？？) ; 所以我修改了一下id编号，改成了4，没有报错（我也试过了3，也可以，不报错） ; 创建文件成功给点提示吧 print, '>>>' + extract_area_array[extract_area_i] + '位置aod提取的txt文件创建成功' ; 每一个点在n多个modis文件里面都有最精确的aod(气溶胶厚度), 所以每一个modis文件都需要求离待提取点的最近的点的aod ; 所以需要循环这些modis文件 ; 得到modis文件的总数量 file_count = n_elements(file_name_array) ; 进入循环 for file_i = 0, file_count - 1 do begin ; 获取当前循环下的modis文件的数据集数据(包括lon、lat、aod)以及aod数据集的两个属性(scale_factor、_FillValue——用于计算真实的aod数据) ; 获取三个数据集(lon、lat、aod) file_lat = get_ds_data(file_path_array[file_i], 'Latitude') file_lon = get_ds_data(file_path_array[file_i], 'Longitude') file_aod = get_ds_data(file_path_array[file_i], 'Image_Optical_Depth_Land_And_Ocean') ; 存储aod的数据集 ; 以上均为数组，应该能理解吧？-因为一个文件里面它不可能只有一个点的经纬度、aod等等数据，一个数据集存储有很多个点的数据(以数组形式存在，当然是二维，你可以使用hdf explorer软件打开看一下) ; 获取aod数据集的两个属性 aod_sf = get_att_data(file_path_array[file_i], 'Image_Optical_Depth_Land_And_Ocean', 'scale_factor') aod_fv = get_att_data(file_path_array[file_i], 'Image_Optical_Depth_Land_And_Ocean', '_FillValue') ; 计算真实aod file_aod = (file_aod ne aod_fv) * file_aod * aod_sf ; 获取所有点于待提取点的距离(两点之间的距离公式) distance = ((file_lon - extract_lon) ^ 2 - (file_lat - extract_lat) ^ 2) ^ 0.5 ; 找到最短距离 min_distance = min(distance) ; 传入一个数组，返回数组里面最小的元素 ; 由于最短距离小于0.1°的才会输出，所以这里加一个if判断 if min_distance lt 0.1 then begin ; 通过最小距离找到距离数组(即distance数组)里面的那个最小距离的下标(这里说清楚一点比较好， 这里下标是这样子的， 从第一行开始数(从0开始)， 数完第一行继续数第二行，一直到找到那个数的位置停下，譬如array(3, 3)中第一行第二列的下标是6) min_pos = where(min_distance eq distance) ; 传入一个数组(min_distance eq distance实际得到的是一个数组，大家想想看是不是。另外其实where会输出数组中所有不为0的数值的下标，记住0是假，其余数值是真，C语言应该学过) ; 通过找到的下标去找该点的经纬度以及aod(气溶胶厚度) min_lon = file_lon[min_pos] min_lat = file_lat[min_pos] min_aod = file_aod[min_pos] ; 判断min_aod是否是0(本来是_FiLLValue属性值，前面eq比较将其转化为0了)，如果是，那么该aod是无效的 if min_aod eq 0.0 then continue ; 跳出循环不输出 ; 获取当前文件的日期 file_date = get_date(file_name_array[file_i]) file_year = file_date[0] file_month = file_date[1] file_day = file_date[2] ; 输出(有日期、经纬度、AOD) printf, 4, file_year, file_month, file_day, min_lon, min_lat, min_aod, format='(i0, "-", i02, "-", i02, ",", 3(f0.3, :, ","))' endif endfor free_lun, 4 ; 记得关闭之前openw打开的文件 ; 成功将提取的aod存储也应该给点提示吧 print, '>>>' + extract_area_array[extract_area_i] + '位置成功提取aod并存储至文件中' endfor end
3.2 运行结果
不用在意为什么成都没有结果，因为我试过了，如果我没有限制aod=0,min_distance<0.1 它们三个城市的输出结果都比较多。所以只是加了比较多的限制
4.题外话
我在运行代码的时候遇到了几个问题
一：出现变量不存在的问题，这个很正常，因为变量实在太多了，难免会混淆。修改一下就好了
譬如：
二：出现错误——》文件已经打开，无法再次打开文件
我是将openw, 1, 路径名 修改为了————》openw, 3, 路径名
然后可以运行了，具体原因我也不清楚，已经做了有一会了，不想代码了
三： 出现浮点数非法，而且没有提示非法在哪个位置我好修改，
不过好在没有影响我的结果输出，我也没有搭理这个警告，因为太难找了
谢谢，刚刚找出来了， 是距离公式的((x - x1) ^ 2 - (y - y1) ^ 2) ^ 0.5错了，你看出来了吗？
修改了距离公式发现还是报错，结果发现更傻缺的错误（现在已经修改过来了）
