1. 实验内容
这里我就不写第三周的课后作业，直接连着这个实验一起写了。
整体思路：
1. 获取Lon、Lat、NO2数据集的数据并对NO2数据做一些处理(填充值处理、单位换算、南北极调换)
2. 求和并计数
3. 求均值
2. 知识储备
由于OMI产品是HDF5文件，需要使用HDF5的相关函数，另外还会使用一些文件处理的函数。
获取HDF5数据的一般步骤：
1. 打开HDF5文件，获取文件id(使用h5f_open()函数)
2. 获取数据集名称(这里你可以通过编程，但是实在没有必要，使用HDF Explorer软件可以查看数据集、属性、还可以查看各个内容数据等，甚至可视化数据)
3. 通过数据集名称以及前面获取的文件id获取数据集的id(通过h5d_open()函数)
4. 获取数据集的数据(通过h5d_read()函数)
5. 关闭数据集和文件(通过h5f_close()、h5d_close()函数)
1、打开HDF5文件，获取文件id
2、获取数据集名
3、将数据集名转换为数据集id
4、从数据集id处读取数据
5、关闭数据集和文件
3. 编程
function get_hdf5_ds, file_path, ds_name ; 构建函数: 获取数据集的数据 ; 获取数据集所在文件的id file_id = h5f_open(file_path) ; 传入数据集所在文件的路径 ; 获取数据集的id ds_id = h5d_open(file_id, ds_name) ; 传入数据集所在文件的id， 传入数据集的名称 ; 获取数据集的数据 ds_data = h5d_read(ds_id) ; 传入数据集的id ; 关闭数据集以及文件 h5d_close, ds_id h5f_close, file_id ; 返回获取数据集的数据 return, ds_data end pro week_three_test start = systime(1) ; 计算NO2的月均值、季均值、年均值,并以Geotiff格式输出、输出单位是mol/km2. ; 要求:程序需要一次完成所有均值的计算 ; 所有文件的路径 in_path = 'D:/IDL_program/experiment_data/chapter_2/NO2' ; 输出文件的路径 out_path = 'D:/IDL_program/experiment_data/chapter_2/NO2/output/' if file_test(out_path, /directory) eq 0 then begin ; 检查该目录是否存在, 不存在则创建 file_mkdir, out_path endif ; 数据集所在group在文件中的路径 group_path = '/HDFEOS/GRIDS/ColumnAmountNO2/Data Fields/' ; 这里使用\后面会报错，应该使用/ ; 数据集的名称 ds_name = 'ColumnAmountNO2TropCloudScreened' ; 数据集在文件中的路径 ds_path = group_path + ds_name ; 获取in_path内所有文件的路径 file_path_array = file_search(in_path, '*.he5', count=file_count) ; count=文件数目 ; 从每一个文件路径获取文件名称 file_name_array = file_basename(file_path_array, '.he5') ; 从文件名称中获取该文件的年份 file_year_array = fix(strmid(file_name_array, 19, 4)) ; 获取起始年份以及年份数 start_year = min(file_year_array) end_year = max(file_year_array) year_count = end_year - start_year + 1 ; 各种数据初始化 ; 存储月份数据集初始化 data_total_month = fltarr(1440, 720, 12) ; 每一个数据集都是(1440, 720), 一年有十二个月份 data_valid_month = fltarr(1440, 720, 12) ; 有的数据集上的NO2数据是无效的，那么我们求平均就不是/ 12,而是看有几次有效次数了 data_aver_month = fltarr(1440, 720, 12) ; 上面的data_total_month / data_valid_month即可得到 ; 存储季节数据集的初始化 data_total_season = fltarr(1440, 720, 4) ; 一年4季 data_valid_season = fltarr(1440, 720, 4) ; 有效次数 data_aver_season = fltarr(1440, 720, 4) ; 平均值 ; 存储年数据集的初始化 data_total_year = fltarr(1440, 720, year_count) ; 一共有year_count年 data_valid_year = fltarr(1440, 720, year_count) data_aver_year = fltarr(1440, 720, year_count) ; 进入循环 for file_i = 0, file_count - 1 do begin ; 检查该文件的NO2数据集是否存在 flag = 0 file_path = file_path_array[file_i] ; 获取该循环下的文件路径 file_id = h5f_open(file_path) ; 获取文件id ds_count = h5g_get_nmembers(file_id, group_path) ; 传入文件id以及group路径返回该group下的数据集个数 for ds_i = 0, ds_count - 1 do begin ds_i_name = h5g_get_member_name(file_id, group_path, ds_i) if ds_i_name eq ds_name then begin flag = 1 continue endif endfor if flag ne 0 then begin ; 提示该文件存在数据集 print, 'Exit NO2 data set>>> ' + file_basename(file_path) ; 根据文件的名称得到该文件的年月日 file_year = fix(strmid(file_basename(file_path), 19, 4)) file_month = fix(strmid(file_basename(file_path), 24, 2)) ; 获取数据集的数据 ds_data = get_hdf5_ds(file_path, ds_path) ; 填充值处理 ds_data = (ds_data gt 0) * ds_data ; 负数(填充值)变成0 ; 单位换算(分子数/cm2 ————》mol/km2) ds_data = (ds_data * 10.0 ^ 10.0) / !const.NA ; 将南北极位置调换 ds_data = rotate(ds_data, 7) ; 存储数据集————》月 data_total_month[*, *, file_month - 1] = data_total_month[*, *, file_month - 1] + ds_data data_valid_month[*, *, file_month - 1] = data_valid_month[*, *, file_month - 1] + (ds_data gt 0) ; 存储数据集————》季 season_index = [3, 3, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3] data_total_season[*, *, season_index[file_month - 1]] = data_total_season[*, *, season_index[file_month - 1]] + ds_data data_valid_season[*, *, season_index[file_month - 1]] = data_valid_season[*, *, season_index[file_month - 1]] + (ds_data gt 0) ; 存储数据集————》年 data_total_year[*, *, file_year - start_year] = data_total_year[*, *, file_year - start_year] + ds_data data_valid_year[*, *, file_year - start_year] = data_valid_year[*, *, file_year - start_year] + (ds_data gt 0) endif else begin print, 'Not exit NO2 data set>> ' + file_basename(file_path) endelse endfor ; 年月日的数据集的均值处理 ; 月处理 data_valid_month = (data_valid_month gt 0) * data_valid_month + (data_valid_month eq 0) * 1.0 data_aver_month = data_total_month / data_valid_month ; 季处理 data_valid_season = (data_valid_season gt 0) * data_valid_season + (data_valid_month eq 0) * 1.0 data_aver_season = data_total_season / data_valid_season ; 年处理 data_valid_year = (data_valid_year gt 0) * data_valid_year + (data_valid_year eq 0) * 1.0 data_aver_year = data_total_year / data_valid_year ; geotiff设置 geoinfo = {$ MODELPIXELSCALETAG:[0.25,0.25,0.0],$ MODELTIEPOINTTAG:[0.0,0.0,0.0,-180.0,90.0,0.0],$ GTMODELTYPEGEOKEY:2,$ GTRASTERTYPEGEOKEY:1,$ GEOGRAPHICTYPEGEOKEY:4326,$ GEOGCITATIONGEOKEY:'GCS_WGS_1984',$ GEOGANGULARUNITSGEOKEY:9102,$ GEOGSEMIMAJORAXISGEOKEY:6378137.0,$ GEOGINVFLATTENINGGEOKEY:298.25722} ; 输出设置 month_out = ['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'] season_out = ['spring', 'summer', 'autumn', 'winter'] ; 依次输出月、季、年均值 for month_i = 0, 11 do begin out_month_path = out_path + 'month_aver_' + month_out[month_i] + '.tiff' write_tiff, out_month_path, data_aver_month[*, *, month_i], /float, geotiff=geoinfo print, file_basename(out_month_path) + ' have completed' endfor for season_i = 0, 3 do begin out_season_path = out_path + 'season_aver_' + month_out[season_i] + '.tiff' write_tiff, out_season_path, data_aver_season[*, *, season_i], /float, geotiff=geoinfo print, file_basename(out_season_path) + ' have completed' endfor for year_i = 0, year_count - 1 do begin out_year_path = out_path + 'year_aver_' + strcompress(string(year_i + start_year), /REMOVE_ALL) + '.tiff' write_tiff, out_year_path, data_aver_year[*, *, year_i], /float, geotiff=geoinfo print, file_basename(out_year_path) + ' have completed' endfor ; 这里再将年结果输出为txt文件 out_txt_path = out_path for year_i = 0, year_count - 1 do begin out_year_path = out_txt_path + 'year_aver_' + strcompress(string(year_i + start_year), /REMOVE_ALL) + '.txt' openw, 2, out_year_path printf, 2, 'lon lat NO2' ; 每一个像元的最中心，一个像元是0.25*0.25 for row_i = 0, 719 do begin lat = 89.875 - 0.25 * row_i for column_i = 0, 1439 do begin lon = -179.875 + column_i * 0.25 printf, 2, lon, lat, data_aver_year[column_i, row_i, year_i], format='(3(f0.3, :, ","))' endfor endfor print, file_basename(out_year_path) + ' have completed' free_lun, 2 endfor stop = systime(1) print, 'Spend time : ' + strcompress(string(stop - start)) end
编译运行结果展示：
生成的tiff文件展示：
使用ENVI打开确实没有问题：
另外需要说一点，这里报错，浮点数非法！但是貌似不影响结果的输出。
虽然但是，我还是心情不舒服。
还有一点就是前面的tiff文件输出是比较快的，但是最后面的txt文件输出很慢，我差点都以为没有运行成功，报时91秒左右。
