• ENVI_IDL: 如何读取HDF5文件再优美地做个均值处理?(以OMI数据为例)

    目录

    01 阅读提示

    02 HDF5函数说明

    2.1 常用的一些函数

    2.2 使用频率不高的函数 

    03 读取HDF5文件数据集输出为TIFF格式

    04 读取HDF5文件数据集并进行均值处理而后输出为GEOTIFF格式 

    01 阅读提示

    我将会以两个案例作为分析，第一个案例即是单纯地读取HDF5文件然后输出为Tiff格式，第二个案例即是读取HDF5文件并进行月、季、年地均值处理并输出为GeoTiff格式。（使用的数据是OMI Level 2G（Grid）数据， 是全球范围的关于NO2的数据）

    当然，如果你觉得阅读案例比较繁琐，你可以选择直接阅读下方关于HDF5文件的相关操作函数的说明。

    02 HDF5函数说明

    2.1 常用的一些函数

    前提说明：

    1. IDL对大小写是不敏感的，所以下面的函数既可以大写也可以小写，我个人比较偏好小写。

    2. 以下说明的函数中h5f(file)表示对hdf5文件进行的相关操作的函数；h5d(dataset)表示对hdf5文件中的数据集进行的相关操作的函数；h5a(attribute)表示对hdf5文件中的数据集的属性或者对hdf5文件中的全局属性进行的相关操作的函数。h5g(Group)表示对hdf5文件中的组进行的相关操作的函数。

    file_id = h5f_open（file_path，/write）

    功能：获取已经存在的hdf5文件，返回该hdf5文件的内存地址/id。

    解释：

    传入hdf5文件在硬盘上的路径；目前还传入了/write参数表示打开hdf5文件的方式是写入，如果不传入/write默认是以只读方式打开hdf5文件。

     dataset_id = h5d_open（file_id，dataset_path）

    功能：获取hdf5文件中的某一个数据集的id

    解释：

    如果你已经学习了如何打开hdf4文件，那么你可能存有疑虑：

    第一个参数传入数据集所在hdf5文件的id这个可以理解，但是第二个参数为什么是传入数据集的路径而不是传入数据集的名称呢？

    这个怎么说呢，我只能说这是人家的规定，如果你觉得不好，你可以自己尝试写一个idl，没有办法。

    这里解释一下第二个参数传入的数据集路径是什么？

    在hdf5文件中也有类似像硬盘上的一个个文件夹的东西，不过IDL将其称之为Geoup（组），不过大同小异，我们的数据集路径就是传入形如Group/Group/Group.../数据集名称 的格式。

    这里有两点需要注意：

    第一个Group你不能传入数据集所在hdf5文件的名称，因为它不是Group。

    第二个在IDL中，硬盘上的路径既可以以反斜杠<\>作为分隔符，也可以以斜杆作为分隔符，但是像Group这种路径则必须使用斜杆作为分隔符。

    dataset_data = h5d_read（dataset_id） 

    功能：获取数据集内的主体数据。

    解释：传入数据集的id号即可。这里的主体数据即是字面意思，它不包含属性内容。

     attr_id = h5a_open_name（dataset_id，attr_name）

    功能：获取属性的id/内存地址。

    解释：第一个参数传入属性所在数据集的id；如果属性是全局属性，那么传入文件id。第二个参数传入属性的名称。

    attr_data = h5a_read（attr_id）

    功能：获取属性内容。

    解释：传入属性的id号即可。这里需要说明一下，一般返回的属性内容比较简短，但是它大部分仍是以数组(array)的形式返回，所以利用属性进行一些数组之间的运算时需要特别注意，否则你会得到意想不到的结果。譬如全部都是数字2的3乘3矩A，我想让矩阵中的每一个数字都乘以属性b值，但是我们直接用获得到的数组形式的属性[a]，那么是[a] * A，最终得到的值是只含有一个值的数组[2a]；如果需要实现上述要求，那么应该是将属性值从数组中取出：[a][0] * A。这里不再赘述，可以自行敲代码验证。

     h5f_close，file_id

    功能：关闭已经打开的hdf5文件。

    解释：传入hdf5文件的id。但是这里需要说明该文件关闭之后，你将无法对该hdf5文件中任何资源进行访问，包括全局属性，数据集属性，数据集内容，组（Group）等等，即使你已经获取到了它们的内存地址（id）。

    其实对于数据集也有对应的关闭函数，但是这里我们一般都是把文件中所有需要获取的内容全部获取到后在直接关闭hdf5文件，而不会操作完某一个数据集关闭该数据集，操作完一个属性关闭该属性。总之一句话，一般我们只需要执行文件关闭的命令而累赘地执行数据集和属性地关闭命令。

    该函数会将hdf5文件在内存中占用的资源全部回收，连同内部的数据集等等，一旦资源被回收，那么你将无法通过id等去操作该文件内部的任何数据。除非你将该文件重新打开，将其从硬盘中重新加载到内存中。

    另外，其它的诸如h5d_close，h5a_close函数也是类似的使用方法，这里不再赘述。

     

    2.2 使用频率不高的函数 

    接下来的仍然是hdf5文件的相关操作函数，但是使用的频率低，只是略提一下。

    h5f_list，file_path，OUTPUT=变量

    功能：获取hdf5文件的相关基本信息。（看下图）

    解释：第一个参数传入hdf5文件的路径；后面的OUTPUT可传可不传，区别在于不传那么结果输出在屏幕上，传入一个变量，那么输出结果赋值给变量而不会输出到屏幕上。 

    attr_count = h5a_get_num_attrs（dataset_id） 

    功能：获取数据集下的属性个数，或者获取hdf5文件中的全局属性个数。

    解释：目前传入的是数据集的id，那么返回数据集下的属性个数，如果传入hdf5文件的id，那么返回hdf5文件下的全局属性的个数。

     attr_id = h5a_open_idx（dataset_id，index）

    功能：获取数据集下索引（index）为index的属性的id，或者获取hdf5文件中索引为index的全局属性的id。

    解释：目前第一个参数传入属性所在数据集的id，第二个参数传入属性的索引 ==> 获取的数据集的索引为index的属性的id；如果需要获取全局属性的id，那么需要第一个参数传入hdf5文件的id，第二个参数仍是传入索引。

    这个需要说明一下，索引和id是不一样的，索引是从0开始。如果有三个属性a，b，c，那么一次索引是0,1,2。但是这种函数一般使用的比较少，因为获取属性的名字我们可以通过软件HDF Explorer查询。

    attr_name = h5a_get_name（attr_id） 

    功能：获取属性的名称

    解释：传入属性的id，函数返回属性的名称。

     dataset_count = h5g_get_members（file_id，group_path）

    功能：获取某一个组下面所有数据集的个数。

    解释：传入hdf5文件的id，那么group_path就需要从最初的group开始写。如果传入的是Group的id，那么group_path可以从该Group之后往下写。

    dataset_name = h5g_get_obj_name_by_idx（group_id，dataset_index）

    功能：获取数据集的名称

    解释：传入数据集所在组的id，传入数据集的index。

    我所说明的函数只是常见或者较为常见的，还有其它函数大家可以自行查看书籍或者官方文档并敲代码实践 ==> 实践出真知。

    03 读取HDF5文件数据集输出为TIFF格式

    1. ; 本函数用于打开HDF5文件并获取数据集
    2. function get_hdf5_ds, file_path, ds_path  ; 输入数据集所在文件的路径以及 --数据集的路径(关于这个看主程序的注释)
    3.   ; 获取文件的id————》通过h5f_open()函数获取
    4.   file_id = h5f_open(file_path)
    5.   ; 传入文件的路径
    6.   
    7.   ; 获取数据集的id————》通过h5d_open()函数获取(另外需要说一下， h5d表示hdf(简称h)5(简称5)文件中的data set(简称d)的相关操作)
    8.   ds_id = h5d_open(file_id, ds_path)
    9.   ; 传入数据集所在文件的id，传入数据集的名称
    10.   
    11.   ; 获取数据集的数据————》通过h5d_read()函数获取
    12.   data = h5d_read(ds_id)
    13.   ; 传入数据集的id
    14.   
    15.   ; 打开的文件以及数据集需要全部关闭当不用的时候，这是习惯和态度问题(其实是为了后续你再使用这些文件避免一些错误和问题)
    16.   h5d_close, ds_id  ; 这里是h5d, 下面是h5f，两个函数是不一样的，需要注意一下
    17.   h5f_close, file_id
    18.   
    19.   ; 返回数据集的数据
    20.   return, data
    21.   
    22.   ; 销毁数据变量(避免空占内存——但是我是存疑的，return之后函数内部的所有变量不都会被销毁吗？还需要自己操作吗)
    23.   ; 现在不需要存疑了，确实不会销毁，需要自己销毁，当然就算你不销毁，只要该程序整个运行完毕后，这些资源会被回收，除非你将其保存在硬盘中。
    24.   ; 但是IDL中似乎对于程序运行占用的内存好像是由限制的，但是具体多少我忘了，只要你过多的几个变量存储同样大量的数据那么会有内存溢出的错误，导致程序运行失败，不过现在你还没有遇到。
    25.   ; 但是养成这样子的习惯不好吗。未雨绸缪嘛
    26.   data = !null  ; 其实这个我感觉就是类似于重新赋值，那么原来的变量对应的数据就会被销毁。
    27. end
    28.  
    29.  
    30. pro week_three_study1
    31.   ; 本程序解决如何读取OMI数据以及它的输出(OMI数据这里是HDF5格式,实际上是学习如何读取和输出HDF5文件)
    32.   
    33.   ; 获取文件路径以及数据集的路径(通过hdf explorer软件轻松获取)
    34.   file_path = 'D:\IDL_program\experiment_data\chapter_2\NO2\2017\OMI-Aura_L3-OMNO2d_2017m0101_v003-2018m0627t042221.he5'
    35.   ds_path = '/HDFEOS/GRIDS/ColumnAmountNO2/Data Fields/ColumnAmountNO2TropCloudScreened'
    36.   ; 这里说明一下，hdf4不需要数据集在文件里面的目录，但是hdf5需要输入数据集在文件里面的目录
    37.   
    38.   ; 获取输出的文件夹路径
    39.   out_dir_path = 'D:\IDL_program\experiment_data\chapter_2\output'
    40.   ; 这里假定你不知道路径中output文件夹是否存在，且你并不愿意亲自查看，可以进行如下操作
    41.   if file_test(out_dir_path, /directory) eq 0 then begin
    42.     ; 这里有一个函数需要认识一下————》file_test()——》检查文件是否存在(默认检查文件，需要检查文件夹需要加上参数/directory)
    43.     ; 其实上面的/directory说的不够准确，应该说是file_test()用来检查文件或者文件夹是否是否存在，加上/directory表示检查文件是否存在并且该文件是否是文件夹 ==> 也就是说即使一个文件存在但是他不是文件夹那么也会返回0
    44.     ; 该函数会返回0或者1，0表示不存在，1表示存在
    45.     
    46.     ; 创建文件夹
    47.     file_mkdir, out_dir_path
    48.   endif
    49.   
    50.   ; 获取输出的文件的名称
    51.   out_name = file_basename(file_path, '.he5') + '.tiff'  ; 输出为tiff格式可以进行查看
    52.   ; 这里file_basename()传入file_path已经可以将文件名称输出，但是我们不需要后缀名.he5
    53.   ; 所以传入字符串'.he5'函数会将已经获取的文件名称中这个字符串删除并返回
    54.   
    55.   ; 输出的路径 = 输出文件夹的路径 + 输出文件的名称
    56.   out_path = out_dir_path + '\' + out_name
    57.   
    58.   ; 获取数据集的数据
    59.   data = get_hdf5_ds(file_path, ds_path)
    60.   
    61.   ; 输出
    62.   write_tiff, out_path, data, /float
    63.   ;传入输出的路径、输出的数据、另外这里默认是输出int型，但是我们通过hdf explorer发现该数据集的type是float，所需要指定一下(这很重要)
    64. end

     

    04 读取HDF5文件数据集并进行均值处理而后输出为GEOTIFF格式 

    这个案例举个可能不够好，因为它并不适合刚学完上述函数的学习者。但是相信你仍可以从中获取一些有用的信息。

    1. ; 获取数据集内容
    2. function get_he5_dataset, file_path, dataset_path
    3.   ; 获取he5文件的id
    4.   file_id = H5F_OPEN(file_path)  ; 传入需要打开的hdf5文件的绝对路径, 不传入打开方式默认以只读方式打开, 传入/write表示可读可写
    5.   
    6.   ; 获取数据集的id
    7.   dataset_id = H5D_OPEN(file_id, dataset_path)  ; 传入参数: 数据集所在文件的id 数据集的名称
    8.   
    9.   ; 获取数据集的内容
    10.   dataset_data = H5D_READ(dataset_id)
    11.   
    12.   ; 关闭打开的hdf5文件
    13.   H5F_CLOSE, file_id
    14.   
    15.   ; 返回数据集内容
    16.   return, dataset_data
    17. end
    18.  
    19. ; 获取数据集属性
    20. function get_he5_attr, file_path, dataset_path, attr_name
    21.   ; 获取he5文件的id
    22.   file_id = H5F_OPEN(file_path)  ; 传入需要打开的hdf5文件的绝对路径, 不传入打开方式默认以只读方式打开, 传入/write表示可读可写
    23.   
    24.   ; 获取数据集的id
    25.   dataset_id = H5D_OPEN(file_id, dataset_path)  ; 传入参数: 数据集所在文件的id 数据集的名称
    26.  
    27.   ; 获取数据集的属性的id
    28.   attr_id = H5A_OPEN_NAME(dataset_id, attr_name)
    29.  
    30.   ; 获取数据集属性的内容
    31.   attr_data = H5A_READ(attr_id)
    32.  
    33.   h5a_  
    34.   ; 关闭打开的hdf5文件
    35.   H5F_CLOSE, file_id
    36.   
    37.   ; 返回数据集属性的内容
    38.   return, attr_data
    39. end
    40.  
    41. ; 判断目标数据集是否存在
    42. function exist_dataset, file_path, group_path, dataset_name
    43.   ; 标记是否存在, 1为存在, 0为不存在
    44.   exit_falg = 0
    45.   
    46.   ; 获取hdf5文件的id
    47.   file_id = h5f_open(file_path)  ; 不传入读写方式默认只读
    48.   
    49.   ; 获取目标数据集所在组的数据集个数
    50.   dataset_num = H5G_GET_NMEMBERS(file_id, group_path)  ; 传入文件的id, 传入组的路径
    51.   
    52.   ; 获取组id
    53.   group_id = H5G_OPEN(file_id, group_path)  ; 传入hdf5文件的id, 传入组的路径
    54.   
    55.   ; 循环搜查, 查看是否存在目标数据集
    56.   for dataset_i = 0, dataset_num - 1 do begin
    57.     ; 获取当前循环下的数据集名称
    58.     dataset_i_name = H5G_GET_OBJ_NAME_BY_IDX(group_id, dataset_i)
    59.     
    60.     ; 查看是否为传入的dataset_name一致
    61.     if dataset_name eq dataset_i_name then begin
    62.       ; 标记赋值
    63.       exit_falg = 1
    64.       ; 结束循环
    65.       break
    66.     endif
    67.   endfor
    68.   
    69.   
    70.   ; 关闭刚刚打开的hdf5文件 ==> 注意, 这里你关闭了整个hdf5文件，就不需要再去关闭什么数据集，组什么的了，因为它们都在hdf5文件中，文件都关闭了，里面的数据集，组都无法访问，即便你已经得到了它们的id
    71.   ; 所以，只有再获取到了文件里面你所需要的所有数据，你才可以去关闭hdf5文件，关闭之后内存中该hdf5文件的资源被回收，如果还需要访问该文件，你需要重新打开该文件将其从硬盘中重新加载到内存中才可以继续访问
    72.   H5F_CLOSE, file_id
    73.   
    74.   ; 返回值
    75.   return, exit_falg
    76. end
    77.  
    78. pro omi_no2_average
    79.   ; 记录一下开始时间
    80.   start_time = SYSTIME(1)
    81.   
    82.   ; 所有he5文件所在的总路径
    83.   files_path = 'D:\IDL_program\experiment_data\chapter_2\NO2'
    84.   
    85.   ; 获取所有he5文件的路径列表
    86.   file_path_list = FILE_SEARCH(files_path, '*NO2*.he5')  ; 限制一下, 需要取文件名称里面含有NO2以及后缀名是.he5的文件, 否则该文件不是我们寻找的对象
    87.   ; 总共的he5文件的数量
    88.   file_path_num = N_ELEMENTS(file_path_list)
    89.   
    90.   ; 需要获取的数据集的路径 --> 这个路径指的是在hdf5里面group/group/数据集名称这样子写, 而不能像hdf4那样子只写数据集名称
    91.   ; 数据集所在的组的路径 ==> 这里可能需要解释一下, 就是反斜杠\和斜杠/, 在上面这种正常的文件路径里面你可以使用\或者/, 但是在hdf5文件里面的路径是只能\的, 否则报错
    92.   group_path = 'HDFEOS/GRIDS/ColumnAmountNO2/Data Fields'
    93.   ; 数据集名称
    94.   no2_name = 'ColumnAmountNO2TropCloudScreened'
    95.   ; 组合得到数据集在hdf5文件内的路径
    96.   no2_path = group_path + '/' + no2_name
    97.   ; 获取ColumnAmountNO2TropCloudScreened数据集的_FillValue属性内容
    98.   no2_fv_name = '_FillValue'
    99.   
    100.   ; 获取年的最大值和最小值以及年数
    101.   ; 获取所有文件的名称列表
    102.   file_name_list = file_basename(file_path_list)
    103.   ; 获取所有年的列表
    104.   years_list = strmid(file_name_list, 19, 4)
    105.   ; 获取最大年(str)
    106.   max_year = max(years_list)
    107.   ; 获取最小年(str)
    108.   min_year = min(years_list)
    109.   ; 获取最大年(int)
    110.   max_year = fix(max_year)
    111.   ; 获取最小年(int)
    112.   min_year = fix(min_year)
    113.   ; 获取年数
    114.   year_num = max_year - min_year + 1  ; 这里注意要加一, 自己也要算进去
    115.   
    116.   ; 所有月数据的事先存储箱子
    117.   months_total_box = fltarr(1440, 720, 12)  ; 1440列, 720行, 12个月
    118.   ; 所有月数据的有效次数的事先存储箱子
    119.   months_valid_box = fltarr(1440, 720, 12)
    120.   ; 月均值的事先存储箱子
    121.   months_aver_box = fltarr(1440, 720, 12)
    122.   ; 所有季数据的事先存储箱子
    123.   seasons_total_box = fltarr(1440, 720, 4)  ; 1440列 720行 4季
    124.   ; 所有季数据的有效次数的事先存储箱子
    125.   seasons_valid_box = fltarr(1440, 720, 4)
    126.   ; 季均值的事先存储箱子
    127.   seasons_aver_box = fltarr(1440, 720, 4)
    128.   ; 所有年数据的事先存储箱子
    129.   years_total_box = fltarr(1440, 720, year_num)
    130.   ; 所有年数据的有效次数的事先存储箱子
    131.   years_valid_box = fltarr(1440, 720, year_num)
    132.   ; 年均值的事先存储箱子
    133.   years_aver_box = fltarr(1440, 720, year_num)
    134.   
    135.   ; 由于获取的数据集的全球数据集, 但是hdf5文件中没有经纬度数据集, 所以需要我们自己构造
    136.   ; 创建经纬度数据的空箱子
    137.   latitude_box = fltarr(1440, 720)
    138.   longitude_box = fltarr(1440, 720)
    139.   ; 往箱子里面存放数据
    140.   for column_i = 0, 1439 do begin
    141.     for row_i = 0, 719 do begin
    142.       latitude_box[column_i, row_i] = 89.875 - row_i * 0.25  ; 为什么是89.875而不是90呢？因为我们想要取每一个像素的最中间的经纬度, 下面的-179.875也是这个原因
    143.       longitude_box[column_i, row_i] = -179.875 + column_i * 0.25
    144.     endfor
    145.   endfor
    146.  
    147.   ; 循环每一个路径，对每一个he5文件进行读取获取其中NO2云柱量数据集的数据并进行均值处理
    148.   for file_path_i = 0, file_path_num - 1 do begin
    149.     ; 获取当前循环下的路径
    150.     file_path = file_path_list[file_path_i]
    151.     
    152.     ; 这里需要说明一下, 就是可能有的hdf5文件中目标数据集可能缺失没有, 所以需要判断一下, 如果没有那么结束本次循环跳转到下一次循环中
    153.     if ~exist_dataset(file_path, group_path, no2_name) then begin
    154.       ; 注意这里我加了一个~表示否定, 和python的not, C语言的!是一个意思. 这里的用法和C语言类似, 在这里 <和> 可以使用and, 也可以使用&&; <或> 可以使用or, 也可以使用||
    155.       ; 结束当前循环
    156.       print, '>>> 文件 ' + file_path + '不存在目标数据集'
    157.       continue
    158.     endif
    159.     
    160.     ; 获取数据集ColumnAmountNO2TropCloudScreened数据集的数据(数组)
    161.     no2_data = get_he5_dataset(file_path, no2_path)
    162.     ; 将数据进行颠倒  ==> 因为我们的OMI数据是全球数据, 而且上面是南极, 下面是北极, 所以需要颠倒
    163.     no2_data = rotate(no2_data, 7)  ; 7表示现在的x等于原来的x, 现在的y等于原来的-y, 也就是上下颠倒, 左右不颠倒
    164.     ; 获取数据集no2_data中的属性_FillValue(array)
    165.     no2_fv = get_he5_attr(file_path, no2_path, no2_fv_name)
    166.     ; 获取数据集no2_data中的属性_FillValue(float)
    167.     no2_fv = no2_fv[0]
    168.     
    169.     ; 通过he5文件的路径获取该文件的日期
    170.     ; 获取he5文件的名称
    171.     file_name = file_basename(file_path, '.he5')  ; '.he5'表示获取最后一个斜杆之后的字符串之后,再将获取到的字符串中的字符串.he5删除(实际上是为了去除后缀名)
    172.     ; 获取日期(str)
    173.     file_year = strmid(file_name, 19, 4)
    174.     file_month = strmid(file_name, 24, 2)
    175.     file_day = strmid(file_name, 26, 2)
    176.     ; 获取日期(int)
    177.     file_year = fix(file_year)  ; 没有int()这个函数, 就是fix()表示转化为int型
    178.     file_month = fix(file_month)
    179.     file_day = fix(file_day)
    180.     
    181.     ; 单位换算: 由于数据的单位是molc(分子数)/cm平方, 需要转化为mol/km平方, 1mol=6.02*10^23=!const.NA
    182.     no2_data = ((no2_data ne no2_fv) * no2_data / !CONST.NA) * (10.0 ^ 10.0)
    183.     ; 严重注意上面的10.0 ^ 10.0, 这里如果是10 ^ 10, 那么你的均值运算一定会出现问题, 就是10 ^ 10 是整数相乘, 但是需要注意, int的范围是比较小的, 所以计算会失败, 得出结果为-7168, 所以换成小数进行^是正确的
    184.     ; 上面的式子让数值为no2_fv变为0, 数值不为该值的值进行单位换算.
    185.     
    186.     ; 存储月数据
    187.     months_total_box[*, *, file_month - 1] += no2_data
    188.     ; 存储月数据的有效次数
    189.     months_valid_box[*, *, file_month - 1] += (no2_data ne 0.0)
    190.     
    191.     ; 存储季数据和季数据的有效次数
    192.     season_list = [3, 3, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3]
    193.     season_i = season_list[file_month - 1]
    194.     ; 开始存储季数据和有效次数
    195.     seasons_total_box[*, *, season_i] += no2_data
    196.     seasons_valid_box[*, *, season_i] += no2_data ne 0.0
    197.     
    198.     ; 存储年数据
    199.     years_total_box[*, *, file_year - min_year] += no2_data
    200.     ; 存储年数据的有效次数
    201.     years_valid_box[*, *, file_year - min_year] += no2_data ne 0.0
    202.   endfor
    203.   
    204.   ; 这里假使有的数据的有效次数为0, 那么待会儿total / valid(求均值)就会出现除数为0的情况, 但是这种计算也是有结果的, 结果为NA(表示无效值)
    205.   months_valid_box = (months_valid_box eq 0.0) + (months_valid_box ne 0.0) * months_valid_box
    206.   seasons_valid_box = (seasons_valid_box eq 0.0) + (seasons_valid_box ne 0.0) * seasons_valid_box
    207.   years_valid_box = (years_valid_box eq 0.0) + (years_valid_box ne 0.0) * years_valid_box
    208.   ; 计算月均值
    209.   for month_i = 0, 11 do begin
    210.     months_aver_box[*, *, month_i] = months_total_box[*, *, month_i] / months_valid_box[*, *, month_i]
    211.   endfor
    212.   ; 计算季均值
    213.   for season_i = 0, 3 do begin
    214.     seasons_aver_box[*, *, season_i] = seasons_total_box[*, *, season_i] / seasons_valid_box[*, *, season_i]
    215.   endfor
    216.   ; 计算年均值
    217.   for year_i = 0, year_num - 1 do begin
    218.     years_aver_box[*, *, year_i] = years_total_box[*, *, year_i] / years_valid_box[*, *, year_i]
    219.   endfor
    220.   
    221.   ; Geotiff的设置
    222.   geo_info={$
    223.     MODELPIXELSCALETAG:[0.25,0.25,0.0],$
    224.     MODELTIEPOINTTAG:[0.0,0.0,0.0,-180.0,90.0,0.0],$
    225.     GTMODELTYPEGEOKEY:2,$
    226.     GTRASTERTYPEGEOKEY:1,$
    227.     GEOGRAPHICTYPEGEOKEY:4326,$
    228.     GEOGCITATIONGEOKEY:'GCS_WGS_1984',$
    229.     GEOGANGULARUNITSGEOKEY:9102,$
    230.     GEOGSEMIMAJORAXISGEOKEY:6378137.0,$
    231.     GEOGINVFLATTENINGGEOKEY:298.25722}
    232.   
    233.   ; 输出路径
    234.   out_path = files_path + '\average'
    235.   ; 判断该路径是否存在, 如果不存在就创建
    236.   if FILE_TEST(out_path, /DIRECTORY) eq 0 then begin
    237.     ; 创建该目录
    238.     file_mkdir, out_path
    239.   endif
    240.   ;  输出月均值
    241.   for month_i = 0, 11 do begin
    242.     out_month_path = out_path + '\month_aver_' + STRCOMPRESS(string(month_i + 1), /REMOVE_ALL) + '.tiff'
    243.     WRITE_TIFF, out_month_path, months_aver_box[*, *, month_i], GEOTIFF=geo_info, /float
    244.     ; 注意/float不要忘记了, 因为默认按照浮点数写入tiff文件, 但是我们的数据是浮点的, 所以需要指定/float告诉write_tiff函数要以float形式写入数据， 下面/float同理
    245.   endfor
    246.   print, '>>> 输出月均值完毕'
    247.   ; 输出季均值
    248.   season_name = ['spring', 'summer', 'autumn', 'winter']
    249.   for season_i = 0, 3 do begin
    250.     out_season_path = out_path + '\season_aver_' + STRCOMPRESS(season_name[season_i], /REMOVE_ALL) + '.tiff'
    251.     WRITE_TIFF, out_season_path, seasons_aver_box[*, *, season_i], GEOTIFF=geo_info, /float
    252.   endfor
    253.   print, '>>> 输出季均值完毕'
    254.   ; 输出年均值
    255.   for year_i = 0, year_num - 1 do begin
    256.     out_year_path = out_path + '\year_aver_' + STRCOMPRESS(string(year_i + min_year), /REMOVE_ALL) + '.tiff'
    257.     WRITE_TIFF, out_year_path, years_aver_box[*, *, year_i], GEOTIFF=geo_info, /float
    258.   endfor
    259.   print, '>>> 输出年均值完毕'
    260.   
    261.   ; 记录一下结束时间
    262.   end_time = SYSTIME(1)
    263.   ; 程序总计用时
    264.   spend_time = strcompress(string(end_time - start_time), /REMOVE_ALL)  ; 里面是为了str化, 外边这个函数是为了将多余的空格删除清空
    265.   ; 提示
    266.   print, '>>> 用时 ' + spend_time + ' s'
    267. end

    我是炒茄子

  • 相关阅读:
    python——GIL锁详解
    防御第六次作业-笔记整理
    ssm基于Springboot的大学生竞赛辅导网站系统java-maven项目
    主机加固如何应对数据世界的绑匪
    OpenGL入门（一）之认识OpenGL和创建Window
    征途服务端安装
    ubuntu16.04安装驱动记录
    SpringCloud | 单体商城项目拆分（微服务）
    软考知识点之需求管理
    【MAPBOX基础功能】21、mapbox绘制自动旋转的图标
  • 原文地址：https://blog.csdn.net/m0_63001937/article/details/126382993