我将以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