我将以MODIS气溶胶产品进行一些简单的分析.如果你觉得稍有繁琐, 你可以直接阅读下方关于HDF4的文件读取的一些常见函数说明.
file_id = hdf_sd_start(file_path, /read)
功能: 打开HDF4文件,返回该文件的id(内存地址,以后访问该文件可以通过其)
解释: 第一个参数传入hdf4文件的绝对路径(如果该hdf4文件在目前pro文件所在文件夹中, 那么相对路径也是可以的),第二个参数传入打开hdf4文件的方式,除了目前传入的/read(只读),你还可以传入/rdwr(可读可写),/create(创建新的hdf4文件)。函数返回该文件的内存地址
hdf_sd_info, file_id, dataset_num(接收数据的变量), attr_num(接收数据的变量)
功能:获取hdf4文件内的数据集个数,全局属性个数
解释:第一个参数传入hdf4文件的id;第二个参数传入一个变量名,用于接收该文件的数据集个数;第三个参数传入一个变量名,用于接收该文件的全局属性个数。
dataset_index = hdf_sd_nametoindex(file_id, dataset_name)
功能:获取hdf4文件中数据集名称为dataset_name在所有数据集列表(不需要了解)中的索引
解释:第一个参数传入hdf4文件的id;第二个参数传入数据集的名称(str类型)。函数返回该数据集的索引(index)
dataset_id = hdf_sd_select(file_id,dataset_index)
功能:获取文件id为file_id的hdf4文件下的数据集索引为dataset_index的数据集id
解释:第一个参数传入文件id;第二个参数传入数据集的索引(index)。函数返回数据集的id(内存地址)
get_sd_getinfo,dataset_id,name=变量1,natts=变量2
功能:获取数据集的名称,该数据集下的属性个数
解释:第一个参数传入数据集的id,后面关键字传变量,显然变量1接收函数返回的该数据集的名称,变量2接收函数返回的该数据集下的属性个数
get_sd_getdata,dataset_id,dataset_data(接收数据集内容的变量)
功能:获取数据集的内容
解释:第一个参数传入数据集的id。后面传入一个变量名dataset_data用于接收函数返回的数据集内容。
get_sd_attrinfo, sd_id, attr_index, name=变量1, data=变量2
功能:获取全局属性的内容和名称 或者 获取数据集的属性和名称
解释:第一个参数传入文件的id,那么第二个参数传入需要获取的全局属性的index;第一个参数传入数据集的id,那么第二个参数传入需要获取的数据集的属性的index。后面传入两个变量分别用于接收函数返回的属性的名称和内容。
hdf_sd_endaccess, dataset_id
功能:关闭数据集,如果打开了数据集一般是需要关闭的。
解释:传入数据集的id
hdf_sd_end, file_id
功能 :关闭hdf4文件,这个是必须要关闭的而不是一般。
解释:传入文件的id
其实前面的数据集关闭一般我是不会去使用的,因为我一般直接关闭整个文件,自然文件里面的数据集也就会被关闭。
接下来是我写的一个案例。
案例需求:提取MODIS气溶胶产品中的Image_Optical_Depth_Land_And_Ocean数据集内分别与北京经纬度(39.90°N,116.40°E)上海经纬度(121.47, 104.06) 成都(104.06, 30.67)最接近的点的AOD有效结果,结果分别输出到不同的txt文件中去。
当然,这个案例可能取的并不恰当,因为我还是用比较多的其它函数。但是,你可以尝试舍去这些函数,精简功能,其实是不影响阅读的。
; 读取数据集内容并返回 function get_dataset, file_path, dataset_name ; 函数名称定义为get_dataset, 传入参数: 文件路径, 数据集名称 ; 获取文件id file_id = hdf_sd_start(file_path, /read) ; 传入文件路径, 打开方式为read ; 获取数据集index dataset_index = HDF_SD_NAMETOINDEX(file_id, dataset_name) ; 传入参数: 文件id, 数据集名称 ; 获取数据集id dataset_id = HDF_SD_SELECT(file_id, dataset_index) ; 传入参数: 文件id, 数据集索引indedx ; 获取数据集内容 HDF_SD_GETDATA, dataset_id, dataset_data ; 传入数据集id, 返回数据集内容并赋值给dataset_data ; 关闭文件 hdf_sd_end, file_id ; 传入文件id ; 返回数据集内容 return, dataset_data end ; 获取数据集属性内容 function get_attr, file_path, dataset_name, attr_name ; 获取文件id file_id = hdf_sd_start(file_path, /read) ; 传入文件路径, 打开方式为read ; 获取数据集index dataset_index = HDF_SD_NAMETOINDEX(file_id, dataset_name) ; 传入参数: 文件id, 数据集名称 ; 获取数据集id dataset_id = HDF_SD_SELECT(file_id, dataset_index) ; 传入参数: 文件id, 数据集索引indedx ; 获取属性index attr_index = HDF_SD_ATTRFIND(dataset_id, attr_name) ; 获取属性内容 HDF_SD_ATTRINFO, dataset_id, attr_index, data=attr_data ; 关闭文件 hdf_sd_end, file_id ; 传入文件id ; 返回属性内容 return, attr_data end pro hdf4_dataset_read ; 获取下方路径下所有的hdf文件的路径 files_path = 'D:\IDL_program\experiment_data\chapter_1\MODIS_2018_mod04_3k' file_path_list = FILE_SEARCH(files_path, '*.hdf') ; 获取file_path_list中元素的个数 file_path_num = N_ELEMENTS(file_path_list) ; aod数据集的全称 aod_name = 'Image_Optical_Depth_Land_And_Ocean' ; 三个待提取点 extract_lon=[116.40,121.47,104.06] extract_lat=[39.90,31.23,30.67] extract_name = ['北京', '上海', '成都'] ; 循环获取每一个待提取点的最近aod数据 for extract_i = 0, 2 do begin ; 当前循环下待提取点的经纬度 longitude = extract_lon[extract_i] latitude = extract_lat[extract_i] ; 当前待提取点得到的aod数据存储路径 aod_path = files_path + '\min_distance_aod_' + extract_name[extract_i] + '.txt' ; 创建该txt文件 openw, 3, aod_path ; 循环获取每一个文件的数据集内容并进行处理 for file_path_i = 0, file_path_num - 1 do begin ; 获取文件路径 file_path = file_path_list[file_path_i] ; 从文件路径中获取日期 ; 获取文件名称 file_name = FILE_BASENAME(file_path, '.hdf') ; 获取日期(str类型) file_date = strmid(file_name, 10, 7) file_year = strmid(file_name, 10, 4) file_days = strmid(file_name, 14, 3) ; 获取日期(int类型) file_year = fix(file_year) file_days = fix(file_days) ; 日期换算 ; 计算上一年最后一天的儒略日 file_date_last_julian = imsl_datetodays(31, 12, file_year - 1) ; 得到file_date的儒略日 file_date_julian = file_date_last_julian + file_days ; 由儒略日得到年月日 IMSL_DAYSTODATE, file_date_julian, file_day, file_month, file_year ; 获取aod数据集内容 aod_data = get_dataset(file_path, aod_name) ; 获取经纬度数据集内容 lat_data = get_dataset(file_path, 'Latitude') long_data = get_dataset(file_path, 'Longitude') ; 对aod数据集进行预处理 ; 获取aod数据集的属性_FillValue aod_fv = get_attr(file_path, aod_name, '_FillValue') aod_fv = aod_fv[0] ; 获取aod数据集的属性scale_factor aod_sf = get_attr(file_path, aod_name, 'scale_factor') aod_sf = aod_sf[0] ; aod数据集处理 aod_data = (aod_data ne aod_fv) * aod_data * aod_sf + (aod_data eq aod_fv) * aod_fv ; 获取距离(39.90°N, 116.40°E)最近的点 ; 求解所有点距离北京的距离 distance_data = ((long_data - longitude) ^ 2 + (lat_data - latitude) ^ 2) ^ 0.5 ; 求解最小距离 min_distance = min(distance_data) ; 求解最小距离所在点的坐标 min_index = where(distance_data eq min_distance) ; 获取最小距离所在点的aod数据 min_distance_aod = aod_data[min_index] min_distance_aod = min_distance_aod[0] ; 获取最小距离所在点的经纬度数据 min_distance_longitude = long_data[min_index] min_distance_longitude = min_distance_longitude[0] min_distance_latitude = lat_data[min_index] min_distance_latitude = min_distance_latitude[0] ; 判断最小距离所在点的aod数据是否有效 ; 距离需要小于0.1°且不为0.0 if ((min_distance gt 0.0) and (min_distance lt 0.1)) then continue ; 判断aod是否为无效的-9999填充值 if min_distance_aod eq aod_fv then continue printf, 3, file_year, file_month, file_day, min_distance_longitude, min_distance_latitude, min_distance_aod, format='(%"%d-%d-%d,%0.4f,%0.4f,%0.4f")' endfor ; 之前给txt文件分配的内存地址应该收回 FREE_LUN, 3 endfor print, '>>> 程序结束!' end