IDL学习笔记（三）OMI数据处理。hdf5文件读取，图像反转，GeoTiff区别，月季年均值计算提取输出，单位转换，运行时间计算

modis Level 2 grid 数据是全球格网化数据。一天的数据全在其中。
modis Level 1 和 2 数据是一景一景的影像。

hdf5文件读取

hdf explor中，看起来像 “ 文件夹 ”一样，通俗可以叫做组。hdf5 与 hdf4不同，hdf5 不能只有数据集的名字，必须也有完整组名，才能被完整找到读出。

读取HDF5数据的一般步骤·
1、打开HDF5文件，获取文件id·
2、获取数据集名
3、将数据集名转换为数据集id
4、从数据集id处读取数据·
5、关闭数据集和文件

hdf5文件读取：

function hdf5_data_get,file_name,dataset_namefile_id = h5f_open(file_name)dataset_id = h5d_open(file_id,dataset_name) ;直接返回id 而不是indexdata = h5d_read(dataset_id)h5f_close,dataset_idh5d_close,file_idreturn,datadata = !null ;销毁data, 清空内存end

注意与hdf4文件读取，相区分：

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,data
end

单位转换，输出hdf5数据集的图像，并检查图像经纬度是否正确，若错误，则进行翻转等操作

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,data
endfunction hdf5_data_get,file_name,dataset_namefile_id = h5f_open(file_name)dataset_id = h5d_open(file_id,dataset_name) ;直接返回id 而不是indexdata = h5d_read(dataset_id)h5d_close,dataset_id ; h5d 关闭数据集h5f_close,file_id ; h5f 关闭文件return,datadata = !null ;销毁data, 清空内存endpro omi_no2_outputfile_name = 'E:/Data/IDL/chapter_2/NO2/2017/OMI-Aura_L3-OMNO2d_2017m0101_v003-2018m0627t042221.he5'dataset_name = '/HDFEOS/GRIDS/ColumnAmountNO2/Data Fields/ColumnAmountNO2TropCloudScreened'out_directory = 'E:/Data/IDL/chapter_2/result/';测试out_directory的目录（/directory）是否存在，若不存在，则创建dir_test = file_test(out_directory, /directory) if dir_test eq 0 then beginfile_mkdir,out_directoryendifout_name = out_directory + file_basename(file_name,'.h5')+'.tiff' ; 舍弃文件路径，只获取基础路径，但同时舍去".he5"data_temp = hdf5_data_get(file_name,dataset_name)data_temp = ((data_temp gt 0.0)*data_temp*10.0^10.0)/(!const.NA);(!const.NA)是阿伏伽德罗常数，该句是转换单位 最后单位是mol/km2;单位换完后，数据会变成“百”“千”量级，方便灰度拉伸;由于图像左右经度是正确的，但上下纬度是颠倒的，因此需要做一个上下反转（用到 rotate）函数data_temp = rotate(data_temp,7)write_tiff,out_name,data_temp,/float;,geotiff = geo_infoend

可以通过envi进行经纬度的查看。

GeoTiff 和 一般的 Tiff 区别

Geotiff和一般tiff的区别：
Geotiff有包含地理坐标、投影相关的描述信息

如何把上述 tiff 改为 GeoTiff 呢？

添加结构体、添加wirte_tiff后面的关键词。

结构体：

  ;投影参数geo_info = {$ ;结构体MODELPIXELSCALETAG:[0.25, 0.25, 0.0],$ ;X /Y /Z 方向的像元分辨率MODELTIEPOINTTAG:[0.0, 0.0, 0.0, -180.0, 90.0, 0.0],$ ;角点坐标。前三个指定位置，后三个指定数值。坐标转换信息，前三个0.0表示栅格图像第0，0，0个像素位置（z方向一般不存在），（若为0.5）就代表中心位置。后面-180.0代表x方向第 0 个位置对应的经度是 -180度，90.0代表y方向第0个位置对应的纬度是90.0度;MODELTIEPOINTTAG:[0.5, 0.0, 0.0, -179.875, 90.0, 0.0],$ ; 当起点中心偏移0.25，那么后面也要改GTMODELTYPEGEOKEY:2,$GTRASTERTYPEGEOKEY:1,$GEOGRAPHICTYPEGEOKEY:4326,$GEOGCITATIONGEOKEY:'GCS_WGS_1984',$GEOGANGULARUNITSGEOKEY:9102,$GEOGSEMIMAJORAXISGEOKEY:6378137.0,$GEOGINVFLATTENINGGEOKEY:298.25722}

一般只需要修改
1.MODELPIXELSCALETAG:[0.25, 0.25, 0.0],$ ;前两个0.25（分辨率，因为一般为2维）
2. MODELTIEPOINTTAG:[0.0, 0.0, 0.0, -180.0, 90.0, 0.0],$ ;一般需要修改-180.0，90.0这两个数字
3. 下面的数字不需要修改。涉及到GeoTiff的存储信息。如果只有上面两个信息，软件读取可能会乱赋值，所以需要写。
4. " $ "代表换行（另起一行）

目标影像下GeoTiff信息如何获取？

获取目标投影geotiff信息的最简方式
1.找一个已有目标投影格式的tiff文件
2.idl读取其中的geotiff结构体

tiff_file = 'D:/world_dem.tif'
data = read_tiff ( tiff_file , geotiff = geo_info )
print, geo_info
help, geo_info

3.修改该结构体中与起始坐标、分辨率的相关参数
4.将修改后的结构体写入目标文件

注意：可以对目标影像进行，重采样，裁剪等操作，看看哪些参数没有改变，无变化参数就是结构体需要的数字。

关于Geotiff的更多说明：

https://gis.stackexchange.com/questions/62837/understanding-geotifftagshttp://duff.ess.washington.edu/data/raster/drg/docs/geotiff.txt

计算不同年份数据，XX月均值计算。季度均值计算。年均值计算。（仅计算有效值，忽略无效值）

function h5_data_get,file_name,dataset_namefile_id=h5f_open(file_name)dataset_id=h5d_open(file_id,dataset_name)data=h5d_read(dataset_id)h5d_close,dataset_idh5f_close,file_idreturn,datadata=!null
endpro omi_no2_average_calculating;输入输出路径设置start_time = systime(1) ; 程序按秒计时，开始计时，与最后的end_time相呼应in_path = 'E:\Data\IDL\chapter_2\NO2\'out_path = 'E:\Data\IDL\chapter_2\NO2\average\'dir_test = file_test(out_path,/directory) ;检测目标文件夹是否存在if dir_test eq 0 then beginfile_mkdir, out_pathendiffilelist = file_search(in_path,'*NO2*.he5') ; 搜索NO2文件下面的，所以hdf5文件file_n = n_elements(filelist)group_name = '/HDFEOS/GRIDS/ColumnAmountNO2/Data Fields/' ; hdf5 需要组的路径，才能读取数据集target_dataset = 'ColumnAmountNO2CloudScreened'dataset_name = group_name + target_dataset; 月份存储数组初始化data_total_month = fltarr(1440,720,12) ; 1440列和720行可以在hdf explor里看到,12分别代表12个月份data_valid_month = fltarr(1440,720,12)data_avr_month = fltarr(1440,720,12); 季节存储数组初始化data_total_season = fltarr(1440,720,4) data_valid_season = fltarr(1440,720,4)data_avr_season = fltarr(1440,720,4); 处理年份设置、年份储存数组初始化year_start = 2017year_n = 2data_total_year = fltarr(1440,720,year_n) data_valid_year = fltarr(1440,720,year_n)data_avr_year = fltarr(1440,720,year_n)for file_i = 0, file_n - 1 do beginprint, filelist[file_i]data_temp = hdf5_data_get(filelist[file_i],dataset_name);help,data_temp ;Array[1440, 720]data_temp = ((data_temp gt 0.0)*data_temp/!const.NA)*(10.0^10.0);转换单位为 mol/km2data_temp = rotate(data_temp,7); 图像反转layer_i = fix(strmid(file_basename(filelist[file_i]),24,2)) - 1 ;从名字中去除月份的信息，从名字中，第24个字母开始取，取两个，但是由于1月要放在0层，2月要放在1层，因此需要-1;月均值data_total_month[*, *, layer_i] = data_temp + data_total_month[*, *, layer_i]data_valid_month[*, *, layer_i] = (data_temp gt 0) + data_valid_month[*, *, layer_i] ;有效是1，无效是0。所有大于0的，记录为1，加入进去，这样，就可以算出每一个位置上出现过几次有效值了;季月均值if(layer_i ge 2) and (layer_i le 4) then begindata_total_season[*, *,  0] = data_tempdata_valid_season[*, *,  0] = (data_temp gt 0.0)endifif(layer_i ge 5) and (layer_i le 7) then begindata_total_season[*, *,  1] = data_tempdata_valid_season[*, *,  1] = (data_temp gt 0.0)endifif(layer_i ge 8) and (layer_i le 10) then begindata_total_season[*, *,  2] = data_tempdata_valid_season[*, *,  2] = (data_temp gt 0.0)endifif(layer_i ge 11) or (layer_i le 1) then begindata_total_season[*, *,  3] = data_tempdata_valid_season[*, *,  3] = (data_temp gt 0.0)endif;年均值year_i = fix(strmid(file_basename(filelist[layer_i]),19,4)) - year_start; 从文件名当中读取年份，减去2017，就只有0或者1（仅两年：2016、2017），且不能是string格式，需要强制类型转换data_total_year[*, *,  year_i] = data_tempdata_valid_year[*, *,  year_i] = (data_temp gt 0.0)endfordata_valid_month = (data_valid_month gt 0.0)*data_valid_month + (data_valid_month eq 0.0)*(1.0)data_valid_season = (data_valid_season gt 0.0)*data_valid_season + (data_valid_season eq 0.0)*(1.0)   data_valid_year = (data_valid_year gt 0.0)*data_valid_year + (data_valid_year eq 0.0)*(1.0)  data_avr_month = data_total_month / data_valid_monthdata_avr_season = data_total_season / data_valid_season   data_avr_year = data_total_year / data_valid_year   geo_info = {$MODELPIXELSCALETAG:[0.25,0.25,0.0],$;x、y、z方向的像元分辨率MODELTIEPOINTTAG:[0.0,0.0,0.0,-180.0,90.0,0.0],$;坐标转换信息，前三个0.0代表栅格图像上的第0,0,0个像元位置（z方向一般不存在），后面-180.0代表x方向第0个位置对应的经度是-180.0度，90.0代表y方向第0个位置对应的纬度是90.0度GTMODELTYPEGEOKEY:2,$GTRASTERTYPEGEOKEY:1,$GEOGRAPHICTYPEGEOKEY:4326,$GEOGCITATIONGEOKEY:'GCS_WGS_1984',$GEOGANGULARUNITSGEOKEY:9102,$GEOGSEMIMAJORAXISGEOKEY:6378137.0,$GEOGINVFLATTENINGGEOKEY:298.25722}for mon_i = 0, 11 do beginout_name = out_path + 'month_avr_' + strcompress(string(mon_i),/remove_all) + '.tiff'write_tiff,out_name, data_avr_month[*, *,mon_i],/float, geotiff = geo_info endforfor season_i = 0, 3 do beginout_name = out_path + 'season_avr_' + strcompress(string(season_i),/remove_all) + '.tiff'write_tiff,out_name, data_avr_month[*, *,season_i],/float, geotiff = geo_infoendforfor year_i = 0, 1 do beginout_name = out_path + 'year_avr_' + strcompress(string(year_i),/remove_all) + '.tiff'write_tiff,out_name, data_avr_month[*, *,year_i],/float, geotiff = geo_infoendforend_time = systime(1)print,'Processing is end, the total time consuming is : ' + strcompress(string(end_time - start_time)) + 's.'; /removeallend

1月份数据在第0层，同理，第N月数据，在第N-1层。

1.strcomptress  函数可以 仅保留一个空格 把多余空格剔除掉。
2.加关键字，/remove_all ，可以把前面所有空格剔除，一个也不保留
3.输出路径不要把最后的/漏掉

程序完成后的提示一定要注意，若出现：

% Program caused arithmetic error: Floating illegal operand

可能是出现了0/0这样的非法操作，这样的原因可能是没有有效数据，导致存出去的结果有NAN的数值存在。所以当分母为0时候，可以进行一些操作，把分母换成其他的数字就可以了

systime这个函数，是系统内置的时间函数，默认输出的是当前时间，systime（1）代表相对于过去了多少秒，可以计算一个程序运行了多少秒。前后均添加该函数，就可以计算时间了。