ICESat-2 ATL03光子数据读取

ICESat-2数据处理的方式一般为将光子数据投影到沿轨隔断和高程的二维空间。如下图：

ATL03数据读取

H5是一种数据存储布局，读取原理就是按照该布局获取数据，这里给出两种读取方式。
ATL03的数据字典：ATL03 Product Data Dictionary (nsidc.org)
使用pandas

1. import warnings
2. import pandas as pd
5. def read_hdf5_atl03_beam_pandas(filename, beam, verbose=False):
6.     # 打开HDF5文件进行读取
7.     h5_store = pd.HDFStore(filename, mode='r')
8.     root = h5_store.root
9.     # 为ICESat-2 ATL03变量和属性分配python字典
10.     atl03_mds = {}
12.     # 读取文件中每个输入光束
13.     # beams = [k for k in file_id.keys() if bool(re.match('gt\\d[lr]', k))]
14.     beams = ['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']
15.     if beam not in beams:
16.         print('请填入正确的光束代码')
17.         return
19.     atl03_mds['heights'] = {}
20.     atl03_mds['geolocation'] = {}
22.     # -- 获取每个HDF5变量
23.     with warnings.catch_warnings():
24.         warnings.simplefilter("ignore")
25.         # -- ICESat-2 Heights Group
26.         heights_keys = ['dist_ph_across', 'dist_ph_along', 'h_ph', 'lat_ph', 'lon_ph', 'signal_conf_ph']
27.         for key in heights_keys:
28.             atl03_mds['heights'][key] = root[beam]['heights'][key][:]
30.         geolocation_keys = ['ref_elev', 'ph_index_beg', 'segment_id', 'segment_ph_cnt', 'segment_dist_x', 'segment_length']
31.         # -- ICESat-2 Geolocation Group
32.         for key in geolocation_keys:
33.             atl03_mds['geolocation'][key] = root[beam]['geolocation'][key][:]
35.     h5_store.close()
36.     return atl03_mds

复制代码

使用h5py

1. import os
2. import h5py
3. import re
5. def read_hdf5_atl03_beam_h5py(filename, beam, verbose=False):
6.     """
7.     ATL03 原始数据读取
8.     Args:
9.         filename (str): h5文件路径
10.         beam (str): 光束
11.         verbose (bool): 输出HDF5信息
13.     Returns:
14.         返回ATL03光子数据的heights和geolocation信息
15.     """
17.     # 打开HDF5文件进行读取
18.     file_id = h5py.File(os.path.expanduser(filename), 'r')
20.     # 输出HDF5文件信息
21.     if verbose:
22.         print(file_id.filename)
23.         print(list(file_id.keys()))
24.         print(list(file_id['METADATA'].keys()))
26.     # 为ICESat-2 ATL03变量和属性分配python字典
27.     atl03_mds = {}
29.     # 读取文件中每个输入光束
30.     beams = [k for k in file_id.keys() if bool(re.match('gt\\d[lr]', k))]
31.     if beam not in beams:
32.         print('请填入正确的光束代码')
33.         return
35.     atl03_mds['heights'] = {}
36.     atl03_mds['geolocation'] = {}
37.     atl03_mds['bckgrd_atlas'] = {}
39.     # -- 获取每个HDF5变量
40.     # -- ICESat-2 Measurement Group
41.     for key, val in file_id[beam]['heights'].items():
42.         atl03_mds['heights'][key] = val[:]
44.     # -- ICESat-2 Geolocation Group
45.     for key, val in file_id[beam]['geolocation'].items():
46.         atl03_mds['geolocation'][key] = val[:]
48.     for key, val in file_id[beam]['bckgrd_atlas'].items():
49.         atl03_mds['bckgrd_atlas'][key] = val[:]
51.     return atl03_mds

复制代码

重建沿轨道隔断

在读取ICESat-2 ATL03数据后，我们必要根据分段信息重建每个光子的沿轨道隔断，已知数据如下：

• ATL03每个分段内的光子都有一个沿轨道隔断（dist_ph_along），该隔断始于当前分段。
  
• ATL03每个分段有一个沿轨道隔断（segment_dist_x），该隔断始于赤道交叉口，即真正的沿轨道隔断。
  
• 当前分段内每个光子的相对沿轨道隔断 + 当前分段的沿轨道隔断 = 每个光子的真正沿轨道隔断。
  

代码如下：

1. import numpy as np
4. def get_atl03_x_atc(atl03_mds):
5.     val = atl03_mds
7.     # 初始化
8.     val['heights']['x_atc'] = np.zeros_like(val['heights']['h_ph']) + np.NaN
9.     val['heights']['y_atc'] = np.zeros_like(val['heights']['h_ph']) + np.NaN
10.     val['geolocation']['ref_elev_all'] = np.zeros_like(val['heights']['h_ph'])
12.     # -- ATL03 Segment ID
13.     segment_id = val['geolocation']['segment_id']
14.     # -- 分段中的第一个光子（转换为基于0的索引）
15.     segment_index_begin = val['geolocation']['ph_index_beg'] - 1
16.     # -- 分段中的光子事件数
17.     segment_pe_count = val['geolocation']['segment_ph_cnt']
18.     # -- 每个ATL03段的沿轨道距离
19.     segment_distance = val['geolocation']['segment_dist_x']
20.     # -- 每个ATL03段的轨道长度
21.     segment_length = val['geolocation']['segment_length']
23.     # -- 对ATL03段进行迭代，以计算40m的平均值
24.     # -- 在ATL03中基于1的索引：无效==0
25.     # -- 此处为基于0的索引：无效==-1
26.     segment_indices, = np.nonzero((segment_index_begin[:-1] >= 0) &
27.                                   (segment_index_begin[1:] >= 0))
28.     for j in segment_indices:
29.         # -- j 段索引
30.         idx = segment_index_begin[j]
31.         # -- 分段中的光子数（使用2个ATL03分段）
32.         c1 = np.copy(segment_pe_count[j])
33.         c2 = np.copy(segment_pe_count[j + 1])
34.         cnt = c1 + c2
36.         # -- 沿轨道和跨轨道距离
37.         # -- 获取当前段光子列表，idx当前段(j)第一个光子数量，c1当前段光子数量，idx+c1当前段长度
38.         distance_along_x = np.copy(val['heights']['dist_ph_along'][idx: idx + cnt])
39.         ref_elev = np.copy(val['geolocation']['ref_elev'][j])
40.         # -- 给当前段的光子加上当前段沿轨道距离
41.         distance_along_x[:c1] += segment_distance[j]
42.         distance_along_x[c1:] += segment_distance[j + 1]
43.         distance_along_y = np.copy(val['heights']['dist_ph_across'][idx: idx + cnt])
45.         val['heights']['x_atc'][idx: idx + cnt] = distance_along_x
46.         val['heights']['y_atc'][idx: idx + cnt] = distance_along_y
47.         val['geolocation']['ref_elev_all'][idx: idx + c1] += ref_elev

复制代码

ATL03数据截取

在处理ATL03时，我们一般都会获取经过研究区域内的光子数据，因此必要对数据举行截取操纵，代码如下：
[code]from glob import globfrom readers.get_ATL03_x_atc import get_atl03_x_atcfrom readers.read_HDF5_ATL03 import read_hdf5_atl03_beam, read_hdf5_atl03_coordinatedef read_data(filepath, beam, mask_lat, mask_lon):    """    读取数据，返回沿轨道隔断和高程隔断    :param filepath: h5文件路径    :param beam: 轨道光束    :param mask_lat: 维度范围    :param mask_lon: 经度范围    :return:    """    atl03_file = glob(filepath)    is2_atl03_mds = read_hdf5_atl03_beam(atl03_file[0], beam=beam, verbose=False)    # 添加沿轨道隔断到数据中    get_atl03_x_atc(is2_atl03_mds)    # 选择范围    d3 = is2_atl03_mds    subset1 = (d3['heights']['lat_ph'] >= min(mask_lat)) & (d3['heights']['lat_ph'] = mask_lon[0])        elif mask_lon[0] is None and mask_lon[1] is not None:            subset1 = subset1 & (d3['heights']['x_atc'] = min(mask_lon)) & (d3['heights']['x_atc'] = min(mask_lat)) & (d3[beam]['lat'] = min(mask_lon)) & (d3[beam]['lon']