IDL学习笔记（二）IDL处理卫星数据

HDF文件

数据集属性

数据集名称，是“：”前的一部分，不是long_name中的文字

通用属性

查看hdf的软件为HDF Explorer。也可以尝试使用Panoply进行简单的数据可视化（该软件也可以读nc文件）。

常用HDF4操作函数

HDF_SD_START函数是打开一个hdf文件。
HDF_SD_FILEINFO函数是查询hdf文件信息的。

HDF_SD_SELECT是选择数据集。
HDF_SD_NAMETOINDEX,已知数据集名称的情况下，获取数据集id号。

HDF_SD_GETDATA，获取指定数据集的id号所包含的数据。
HDF_SD_ATTRINFO，获取hdf属性信息。

使用完之后必须要关闭这个文件，不然的话相当于这个文件一直被占用，需要关闭、释放后掉。

常用的HDF5操作函数

读取HDF文件的一般步骤

1.指定文件名

2.打开文件

3.获取数据集/属性名称

4.获取数据集/属性ID

5.获取数据集/属性内容

6.下一步处理

HDF4文件读取-----数据信息查询

流程：
指定文件名–打开文件–获取文件数据集、属性的个数—循环获取数据集信息—（循环）—选中当前循环的数据集----获取当前数据集名称、属性个数、循环获取数据集属性信息----（是否完成循环）----否----获取属性名。-----是----判断数据集循环是否结束，若是则释放文件，若否循环继续。最后，释放文件

把光标放在函数上，按f1会直接跳转帮助

hdf_sd_fileinfo,sd_id,sds_num,att_num
文件名，文件id，（再任意定义两个文件名装返回结果），DataSets(数据集个数）, Attributes(属性个数）--------获取文件数据集/属性个数

循环获取数据集信息：

pro hdf4_file_info_extractmodis_data = 'E:\Data\modis\2019\05\MOD09GA.A2019130.h27v05.061.2020293213719.hdf'sd_id = hdf_sd_start(modis_data,/read)hdf_sd_fileinfo,sd_id,sds_num,att_num;print,sds_num;print,att_numfor sds_i = 0, sds_num - 1 do beginsds_id = hdf_sd_select(sd_id, sds_i)hdf_sd_getinfo,sds_id,name = sds_name,natts = sds_att_num ;此时这一行的左边赋值给右边，name是关键字，特例print, 'The sds name of '+strcompress(string(sds_i))+':'print,sds_att_numprint,'The include attributes name:'if sds_att_num gt 0 then begin ;循环获取属性信息for att_i=0,sds_att_num-1 do beginhdf_sd_attrinfo, sds_id, att_i, name=sds_att_name;通过name关键字传递属性名称给右边的sds_att_nameprint,sds_att_nameendforprint,string('****************************************')endifhdf_sd_endaccess, sds_id;关闭释放数据集endforend

1. sds_id = hdf_sd_select(sd_id, sds_i)选中数据集.在sd_id文件下，sds_i的地址id是多少传递给sds_id
2. hdf_sd_getinfo,sds_id,name = sds_name,natts = sds_att_num获取在sds_id位置上，文件的名字是什么，在这个文件下，又有多少数据集的属性att_num，赋值给natts
3. hdf_sd_getinfo,sds_id,name = sds_name,natts = sds_att_num ;此时这一行的左边赋值给右边，name是关键字，特例
4. hdf_sd_attrinfo, sds_id, att_i, name=sds_att_name;通过name关键字传递属性名称给右边的sds_att_name
5. hdf_sd_endaccess, sds_id;关闭释放数据集
6. hdf_sd_attrinfo, sds_id, att_i, name=sds_att_name只是显示了数据集名字，如果想显示内容，可以加一个关键词，data = sds_data_info,获取到的信息很多时候是一个数组

pro hdf4_file_info_extractmodis_data = 'E:\Data\modis\2019\05\MOD09GA.A2019130.h27v05.061.2020293213719.hdf'sd_id = hdf_sd_start(modis_data,/read)hdf_sd_fileinfo,sd_id,sds_num,att_num;print,sds_num;print,att_numfor sds_i = 0, sds_num - 1 do beginsds_id = hdf_sd_select(sd_id, sds_i)hdf_sd_getinfo,sds_id,name = sds_name,natts = sds_att_num ;此时这一行的左边赋值给右边，name是关键字，特例print, 'The sds name of '+strcompress(string(sds_i))+':'print,sds_att_numprint,'The include attributes name:'if sds_att_num gt 0 then begin ;循环获取属性信息for att_i=0,sds_att_num-1 do beginhdf_sd_attrinfo, sds_id, att_i, name=sds_att_name,data = sds_data_info;通过name关键字传递属性名称给右边的sds_att_nameprint,sds_att_namehelp,sds_attnamehelp,sds_data_infoendforprint,string('****************************************')endifhdf_sd_endaccess, sds_id;关闭释放数据集endforend

如果想把sds_att_name、sds_attname同时输出的时候呢？

    if sds_att_num gt 0 then begin ;循环获取属性信息for att_i=0,sds_att_num-1 do beginhdf_sd_attrinfo, sds_id, att_i, name=sds_att_name,data = sds_data_info;通过name关键字传递属性名称给右边的sds_att_nameprint,sds_att_name + sds_data_infoendforprint,string('****************************************')endif

此时有内部转换错误，因为字符串和数据集内变量相加，不能直接相加，需要强制类型转换为string

    if sds_att_num gt 0 then begin ;循环获取属性信息for att_i=0,sds_att_num-1 do beginhdf_sd_attrinfo, sds_id, att_i, name=sds_att_name,data = sds_data_info;通过name关键字传递属性名称给右边的sds_att_nameprint,sds_att_name + string(sds_data_info)endforprint,string('****************************************')endif

此时可以相加，但是结果如图所示

valid_range           0
valid_range  2147483647
_FillValue          -1

中间有很长的空格，此时是因为强制类型转换后，把转换前前面有多少个位也显示出来了，因此需要添加函数strcompress，使得空格不再那么长

    if sds_att_num gt 0 then begin ;循环获取属性信息for att_i=0,sds_att_num-1 do beginhdf_sd_attrinfo, sds_id, att_i, name=sds_att_name,data = sds_data_info;通过name关键字传递属性名称给右边的sds_att_nameprint,sds_att_name + strcompress(string(sds_data_info))endforprint,string('****************************************')endif

此时结果：

The include attributes name:
long_nameNumber of additional observations per row
unitsnone
valid_range 0 valid_range 2147483647
_FillValue -1

如果一个空格都不想要，就可以添加关键字“/remove_all”

    if sds_att_num gt 0 then begin ;循环获取属性信息for att_i=0,sds_att_num-1 do beginhdf_sd_attrinfo, sds_id, att_i, name=sds_att_name,data = sds_data_info;通过name关键字传递属性名称给右边的sds_att_nameprint,sds_att_name + strcompress(string(sds_data_info),/remove_al)endforprint,string('****************************************')endif

结果为：

The include attributes name:
long_nameNumberofadditionalobservationsperrow
unitsnone
valid_range0 valid_range2147483647
_FillValue-1

此时结果不好区分，可以认为在中间加一个冒号

 print,sds_att_name + ':'+strcompress(string(sds_data_info),/remove_al)

获取全局信息，以及释放文件：

pro hdf4_file_info_extractmodis_data = 'E:\Data\modis\2019\05\MOD09GA.A2019130.h27v05.061.2020293213719.hdf'sd_id = hdf_sd_start(modis_data,/read)hdf_sd_fileinfo,sd_id,sds_num,att_num;print,sds_num;print,att_numfor sds_i = 0, sds_num - 1 do beginsds_id = hdf_sd_select(sd_id, sds_i)hdf_sd_getinfo,sds_id,name = sds_name,natts = sds_att_num ;此时这一行的左边赋值给右边，name是关键字，特例print, 'The sds name of '+strcompress(string(sds_i))+':'print,sds_att_numprint,'The include attributes name:'if sds_att_num gt 0 then begin ;循环获取属性信息for att_i=0,sds_att_num-1 do beginhdf_sd_attrinfo, sds_id, att_i, name=sds_att_name,data = sds_data_info;通过name关键字传递属性名称给右边的sds_att_nameprint,sds_att_name + ':'+strcompress(string(sds_data_info),/remove_al)endforprint,string('****************************************')endifhdf_sd_endaccess, sds_id;关闭释放数据集endforfor att_i = 0, att_num-1 do beginhdf_sd_attrinfo,sd_id,att_i, name = att_nameprint,'global att name:'print,att_nameendforhdf_Sd_end,sd_idend

HDF4文件读取示例-----目标数据TIFF输出

数组和数字 sf[0] 相乘 VS 数组和数组 sf 相乘 注意区分！

pro hdf4_dataset_readfilename = 'E:\Data\modis\2019\05\MOD09GA.A2019130.h27v05.061.2020293213719.hdf'result_name = 'E:\Data\modis\2019\result\MOD09GA.A2019130.h27v05.061.2020293213719.tiff'modis_sd_id = hdf_sd_start(filename,/read)modis_sds = 'sur_refl_b02_1'modis_sds_index = hdf_sd_nametoindex(modis_sd_id,modis_sds);找到文件地址modis_sds_id = hdf_sd_select(modis_sd_id,modis_sds_index);转换为文件id号hdf_sd_getdata,modis_sds_id,modis_band2 ;由id找到文数据集，赋值给modis_band2modis_sds = 'sur_refl_b03_1'modis_sds_index = hdf_sd_nametoindex(modis_sd_id,modis_sds)modis_sds_id = hdf_sd_select(modis_sd_id,modis_sds_index)hdf_sd_getdata,modis_sds_id,modis_band3modis_att_index = hdf_sd_attrfind(modis_sds_id,'scale_factor')hdf_sd_attrinfo, modis_sds_id,modis_att_index,data=sfmodis_att_index = hdf_sd_attrfind(modis_sds_id,'_FillValue')hdf_sd_attrinfo,modis_sds_id,modis_att_index,data=fvhdf_sd_endaccess,modis_sds_idhdf_SD_end,modis_sd_idmodis_target_data = modis_band2modis_target_data = (modis_target_data ne fv[0])*modis_target_data*sf[0] ;注意这里乘以的是一个数，而不是一个恶搞数组write_tiff,result_name, modis_target_data,/float ;必须要加一个 /float 否则可能输出的是一个黑白图end

保存出去的图为：

保存成功！

提取modis产品中数据，与某一点经纬度最接近的点有效结果，并按每行内容为日期号、经度、纬度 AOD格式输出到IDL控制台

out_data = strmid(file_basename(file_list[file_i],10,7):

1.从file_list[file_i]中， 第10个位置， 提取7个字符出来，赋值给out_data
2.file_basename()意思是去除路径，只在文件名进行读取操作

 print，outdata,modis_lon_data,modis[pos],modis_lat_data[pos]

此时输出的是一个数字一行，是1行3列的模式输出，那么如何能让3个数字都在一行呢？

  print，[outdata,modis_lon_data,modis[pos],modis_lat_data[pos]]

这样就会以数组形式输出。但前面的日期会自动变成浮点数，该如何解决呢？

  print，outdata, [modis_lon_data,modis[pos],modis_lat_data[pos]]

或

  print，outdata, string([modis_lon_data,modis[pos],modis_lat_data[pos]])

就可以解决了

输出格式用逗号分隔：

  print，outdata, [modis_lon_data,modis[pos],modis_lat_data[pos]]，format = '(A,",",3(F0.3,:,","))'

format = ‘(A,“,”,3(F0,:,“,”))’ 格式控制， 第一个A代表是第一个以字符串形式输出，后接一个 ， 号，然后输出三个浮点数。F0.3意思是保留三位小数。此时其实不加 [modis_lon_data,modis[pos],modis_lat_data[pos]] 的方括号也是可以的。：的意思是，如果换行的话，就不要后面的那个逗号

如何把modis的年积日，换成日期号？

IDL有函数，可以把儒略日自动转换为年月日。也有把一个日期转换为儒略日的功能。
如果modis产品日期写的是2018121，这就是意味着，从2017年12月31号过了121天，那么就只需要先算出2017.12.31的儒略日，就可以转而计算2018121对应的儒略日了。
strmid(file_basename(file_list[file_i],10,7),是字符串
fix()把目标转换为整形数

  print，out_year,out_month,out,day, [modis_lon_data,modis[pos],modis_lat_data[pos]]，format = '(I0，“-”, I02，“-”, I02，“-”, I02,",",3(F0.3,:,","))'

1.I0 表示保留数据本身的样子，本身自己前面的空格全部去除
2.“-”意思是中间以 - 链接
3.I02保留本身的5，但需要补足2位，2表示有几个字符占位，若缺少，自动以0补全

function hdf4_data_get, file_name,sds_namesd_id = hdf_sd_start(file_name,/read)sds_index = hdf_sd_nametoindex(sd_id,sds_name)sds_id = hdf_sd_select(sd_id,sds_index)hdf_sd_getdata,sds_id,datahdf_sd_endaccess,sds_idhdf_sd_end,sd_idreturn,dataendfunction hdf4_attdata_get, file_name, sds_name, att_namesd_id = hdf_sd_start(file_name)sds_index = hdf_sd_nametoindex(sd_id,sds_name)sds_id = hdf_sd_select(sd_id, sds_index)att_index = hdf_sd_attrfind(sds_id,att_name)hdf_sd_attrinfo, sds_id,att_index,data=att_datahdf_sd_endaccess,sds_idhdf_sd_end,sd_idreturn, att_dataendpro mod04_nearest_point_value_extracting;程序功能：从modis-mod04气溶胶数据集中提取出特定经纬度点位的产品结果;待提取点位坐标extract_lon = 116.40extract_lat = 39.90point_name = 'BeiJing'data_path = 'E:/Data/IDL/chapter_1/MODIS_2018_mod04_3k/'file_list = file_search(data_path,'*.hdf',count = file_n) ;查找data_path路径下，所有以hdf结尾的文件（包括子路径下的）print,file_listprint,file_nfile_n = n_elements(file_list) ;计算文件数量print,file_n ;两种方式获取 file_nout_file = data_path+'point_value_'+point_name+'.txt'openw,1,out_file ;width = 80000,/appendfor file_i = 0, file_n-1 do beginout_data = strmid(file_basename(file_list[file_i]),10,7)date = fix(strmid(file_basename(file_list[file_i]),14,3))out_year = strmid(file_basename(file_list[file_i]),10,4)out_year_fix = fix(strmid(file_basename(file_list[file_i]),10,4))date_julian = imsl_datetodays(31,12,out_year_fix-1)imsl_daystodate,date_julian + date, day, month, yearout_month = monthout_day = dayprint,out_year,out_month,out_dayprint,file_list[file_i]help,file_list[file_i]modis_lon_data = hdf4_data_get(file_list[file_i],'Longitude')modis_lat_data = hdf4_data_get(file_list[file_i],'Latitude')modis_aod_data = hdf4_data_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean')scale_factor=hdf4_attdata_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean','scale_factor')fill_value = hdf4_attdata_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean','_FillValue')modis_aod_data = (modis_aod_data ne fill_value[0])*modis_aod_data*scale_factor[0]x = modis_lon_data - extract_lony = modis_lat_data - extract_latdistance = sqrt(x^2 + y^2)min_dis = min(distance)pos = where(distance eq min_dis)extract_aod = modis_aod_data[pos]if extract_aod gt 0.0 then beginprint,out_data,modis_lon_data[pos],modis_lat_data[pos],extract_aod,format = '(A,",",3(F0.3,:,","))'print,out_year,out_month,out_day,day, modis_lon_data[pos],modis_lat_data[pos],format ='(I0,"-",I02,"-",I02,",",3(F0.3,:,","))'endifendforfree_lun,1end

注意变量名不要写错，不要漏写“-”，注意英文标点不要写成中文。

其中将日期转换功能单独写为一个功能函数为：

  function hdf4_modis_data_read,file_base_nameout_date=strmid(file_basename(file_base_name),10,7)date=fix(strmid(file_basename(file_base_name),14,3))out_year=fix(strmid(file_basename(file_base_name),10,4))date_julian=imsl_datetodays(31,12,out_year-1) imsl_daystodate,date_julian+date,day,month,yearout_month=monthout_day=dayreturn,[out_year,out_month,out_day]end

调用函数为：

for file_i = 0, file_n-1 do begindata= hdf4_modis_data_read(file_list[file_i])out_year=data[0]out_month=data[1]out_day=data[2]print,out_year,out_month,out_dayprint,file_list[file_i]help,file_list[file_i]modis_lon_data = hdf4_data_get(file_list[file_i],'Longitude')modis_lat_data = hdf4_data_get(file_list[file_i],'Latitude')modis_aod_data = hdf4_data_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean')scale_factor=hdf4_attdata_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean','scale_factor')fill_value = hdf4_attdata_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean','_FillValue')modis_aod_data = (modis_aod_data ne fill_value[0])*modis_aod_data*scale_factor[0]x = modis_lon_data - extract_lony = modis_lat_data - extract_latdistance = sqrt(x^2 + y^2)min_dis = min(distance)pos = where(distance eq min_dis)extract_aod = modis_aod_data[pos]if extract_aod gt 0.0 then beginprint,out_year,out_month,out_day, modis_lon_data[pos],modis_lat_data[pos],extract_aod,format = '(A,",",3(F0.3,:,","))'print,out_year,out_month,out_day, modis_lon_data[pos],modis_lat_data[pos],format ='(I0,"-",I02,"-",I02,",",3(F0.3,:,","))'endifendforfree_lun,1end

(多数据）筛选数据

function hdf4_data_get, file_name,sds_namesd_id = hdf_sd_start(file_name,/read)sds_index = hdf_sd_nametoindex(sd_id,sds_name)sds_id = hdf_sd_select(sd_id,sds_index)hdf_sd_getdata,sds_id,datahdf_sd_endaccess,sds_idhdf_sd_end,sd_idreturn,dataendfunction hdf4_attdata_get, file_name, sds_name, att_namesd_id = hdf_sd_start(file_name)sds_index = hdf_sd_nametoindex(sd_id,sds_name)sds_id = hdf_sd_select(sd_id, sds_index)att_index = hdf_sd_attrfind(sds_id,att_name)hdf_sd_attrinfo, sds_id,att_index,data=att_datahdf_sd_endaccess,sds_idhdf_sd_end,sd_idreturn, att_dataendfunction hdf4_modis_data_read,file_base_nameout_date=strmid(file_basename(file_base_name),10,7)date=fix(strmid(file_basename(file_base_name),14,3))out_year=fix(strmid(file_basename(file_base_name),10,4))date_julian=imsl_datetodays(31,12,out_year-1) ;IMSL相关的函数说明在ENVI安装目录../Exelis/IDL85/help/pdf/advmathstats.pdf中imsl_daystodate,date_julian+date,day,month,yearout_month=monthout_day=dayreturn,[out_year,out_month,out_day]endpro mod04_nearest_point_value_extracting;程序功能：从modis-mod04气溶胶数据集中提取出特定经纬度点位的产品结果;待提取点位坐标extract_lon=[116.40,121.47,104.06]extract_lat=[39.90,31.23,30.67]point_name = ['BeiJing','ShangHai','ChengDu']data_path = 'E:/Data/IDL/chapter_1/MODIS_2018_mod04_3k/'file_list = file_search(data_path,'*.hdf',count = file_n) ;查找data_path路径下，所有以hdf结尾的文件（包括子路径下的）;print,file_list;print,file_nfile_n = n_elements(file_list) ;计算文件数量;print,file_n ;两种方式获取 file_nfor num_point_i = 0, n_elements(point_name)-1 do beginout_file = data_path+'point_value_'+point_name[num_point_i]+'.txt'print,out_fileopenw,1,out_file,/append ;width = 80000,/appendfor file_i = 0, file_n-1 do beginout_data= hdf4_modis_data_read(file_list[file_i])out_year=out_data[0]out_month=out_data[1]out_day=out_data[2]print,out_year,out_month,out_dayprint,file_list[file_i]help,file_list[file_i]modis_lon_data = hdf4_data_get(file_list[file_i],'Longitude')modis_lat_data = hdf4_data_get(file_list[file_i],'Latitude')modis_aod_data = hdf4_data_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean')scale_factor=hdf4_attdata_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean','scale_factor')fill_value = hdf4_attdata_get(file_list[file_i],'Image_Optical_Depth_Land_And_Ocean','_FillValue')modis_aod_data = (modis_aod_data ne fill_value[0])*modis_aod_data*scale_factor[0];最小距离像元位置下标计算lon_minus = abs(modis_lon_data - extract_lon[num_point_i])lat_minus = abs(modis_lat_data - extract_lat[num_point_i])distance = sqrt(lon_minus^2+lat_minus^2)min_pos = where(distance eq min(distance)) ;不是直接输出min(distance)if (modis_aod_data[min_pos] le 0.0) then continueprint,'The three output formats of the nearest point to' + point_name[num_point_i] + 'in file' + file_basename(file_list[file_i]) + ':'print,out_year,out_month,out_day, string([modis_lon_data[min_pos[0]],modis_lat_data[min_pos[0]],modis_aod_data[min_pos[0]]])print,out_year,out_month,out_day, modis_lon_data[min_pos],modis_lat_data[min_pos],format = '(A,",",3(F0.3,:,","))'print,out_year,out_month,out_day, modis_lon_data[min_pos],modis_lat_data[min_pos],format ='(I0,"-",I02,"-",I02,",",3(F0.3,:,","))'print,'*************************************************************************************************************************************************************************'if min(distance) gt 0.1 then continueprint,1,out_year,out_month,out_day, modis_lon_data[min_pos[0]],modis_lat_data[min_pos[0]],format ='(I0,"-",I02,"-",I02,",",3(F0.3,:,","))'; 注意写入的是数组里的数值，而不是数组！否则文件中大概率为空值！endforfree_lun,1endforend