ICESat-2 从ATL08中获取ATL03分类效果

ICESat-2 ATL03数据和ATL08数据的分段距离不一致，ATL08在ATL03的基础上重新分段，并对分段内的数据做处置惩罚得到一系列的效果，详情见数据字典：
ATL08 Product Data Dictionary (nsidc.org)
ATL08使用DRAGANN算法对ATL03数据做了去噪处置惩罚，并使用分类算法对每个光子进行分类
标志值标志含义-1未分类0噪声1地面2冠层3冠顶ATL08使用ph_segment_id和classed_pc_indx可以和ATL03对应起来。基于此，可从ATL08中获取ATL03每个光子的分类信息。

读取ATL08

1. import os
2. import h5py
3. import re
6. def read_hdf5_atl08(filename, beam, verbose=False):
7.     file_id = h5py.File(os.path.expanduser(filename), 'r')
9.     # 输出HDF5文件信息
10.     if verbose:
11.         print(file_id.filename)
12.         print(list(file_id.keys()))
13.         print(list(file_id['METADATA'].keys()))
14.     # 为ICESat-2 ATL08变量和属性分配python字典
15.     atl08_mds = {}
17.     # 读取文件中每个输入光束
18.     beams = [k for k in file_id.keys() if bool(re.match('gt\\d[lr]', k))]
19.     if beam not in beams:
20.         print('请填入正确的光束代码')
21.         return
23.     atl08_mds['signal_photons'] = {}
24.     # -- ICESat-2 Geolocation Group
25.     for key, val in file_id[beam]['signal_photons'].items():
26.         atl08_mds['signal_photons'][key] = val[:]
28.     return atl08_mds

复制代码

映射ATL08

将 ATL08 映射到 ATL03

1. def get_atl08_mapping(atl03_ph_index_beg, atl03_segment_id, atl08_classed_pc_indx,
2.                       atl08_classed_pc_flag, atl08_segment_id):
3.     """
4.     Function to map ATL08 to ATL03 class photons
5.     Args:
6.         atl03_ph_index_beg:
7.         atl03_segment_id:
8.         atl08_classed_pc_indx:
9.         atl08_classed_pc_flag:
10.         atl08_segment_id:
12.     Returns:
14.     """
15.     # Get ATL03 data
16.     indsNotZero = atl03_ph_index_beg != 0
17.     atl03_ph_index_beg = atl03_ph_index_beg[indsNotZero]
18.     atl03_segment_id = atl03_segment_id[indsNotZero]
20.     # Find ATL08 segments that have ATL03 segments
21.     atl03SegsIn08TF, atl03SegsIn08Inds = ismember(atl08_segment_id, atl03_segment_id)
23.     # Get ATL08 classed indices and values
24.     atl08classed_inds = atl08_classed_pc_indx[atl03SegsIn08TF]
25.     atl08classed_vals = atl08_classed_pc_flag[atl03SegsIn08TF]
27.     # Determine new mapping into ATL03 data
28.     atl03_ph_beg_inds = atl03SegsIn08Inds
29.     atl03_ph_beg_val = atl03_ph_index_beg[atl03_ph_beg_inds]
30.     newMapping = atl08classed_inds + atl03_ph_beg_val - 2
32.     # Get max size of output array
33.     sizeOutput = newMapping[-1]
35.     # Pre-populate all photon classed array with zeroes
36.     allph_classed = (np.zeros(sizeOutput + 1)) - 1
38.     # Populate all photon classed array from ATL08 classifications
39.     allph_classed[newMapping] = atl08classed_vals
41.     # Return all photon classed array
42.     return allph_classed

复制代码

添加分类信息

1. def add_atl08_classed_flag(filepath_08, beam, atl03_mod):
2.     """
3.     添加ATL08分类数据到ATL03中
4.     Args:
5.         filepath_08: ATL08数据文件位置
6.         beam: 波束，与ATL03保持一致
7.         atl03_mod: ATL03数据
9.     Returns:
10.     携带ATL08分类信息
11.     """
12.     val_03 = atl03_mod
13.     val_08 = read_hdf5_atl08(filepath_08, beam)
15.     # val_03['classed_pc_flag'] = np.zeros_like(val_03['heights']['h_ph']) + np.NaN
16.     atl03_heights = val_03['heights']['h_ph']
18.     # -- 分段中的第一个光子（转换为基于0的索引）
19.     segment_index_begin = val_03['geolocation']['ph_index_beg']
20.     segment_id = val_03['geolocation']['segment_id']
22.     # 追踪到ATL03上特定20m Segment_ID的光子的段ID
23.     ph_segment_id = val_08['signal_photons']['ph_segment_id']
25.     # 该索引追溯到ATL03上20m segment_id内的特定光子。
26.     classed_pc_index = val_08['signal_photons']['classed_pc_indx']
27.     # 每个光子的陆地植被ATBD分类标志为噪声、地面、树冠和树冠顶部。0=噪音，1=地面，2=冠层，或3=冠层顶部
28.     classed_pc_flag = val_08['signal_photons']['classed_pc_flag']
30.     # Map ATL08 classifications to ATL03 Photons
31.     all_ph_classed = get_atl08_mapping(segment_index_begin, segment_id,
32.                                        classed_pc_index, classed_pc_flag, ph_segment_id)
34.     if len(all_ph_classed) < len(atl03_heights):
35.         n_zeros = len(atl03_heights) - len(all_ph_classed)
36.         zeros = np.zeros(n_zeros)
37.         all_ph_classed = np.append(all_ph_classed, zeros)
39.     val_03['classed_pc_flag'] = all_ph_classed

复制代码

使用姿势

读取ATL03数据代码见：https://www.cnblogs.com/sw-code/p/18161987

1. from glob import glob
3. import numpy as np
4. from matplotlib import pyplot as plt
5. from matplotlib.ticker import MultipleLocator
7. from readers.add_atl08_info import add_atl08_classed_flag
8. from readers.get_ATL03_x_atc import get_atl03_x_atc
9. from readers.read_HDF5_ATL03 import read_hdf5_atl03_beam_h5py
12. def select_atl03_data(atl03_data, mask):
13.     """
14.     选择数据范围
15.     Args:
16.         atl03_data: 所有数据
17.         mask (list): 维度范围
18.     Returns:
19.     """
20.     # 选择范围
21.     d3 = atl03_data
22.     subset1 = (d3['heights']['lat_ph'] > min(mask)) & (d3['heights']['lat_ph'] < max(mask))
24.     x_act = d3['heights']['x_atc'][subset1]
25.     h = d3['heights']['h_ph'][subset1]
26.     signal_conf_ph = d3['heights']['signal_conf_ph'][subset1]
27.     lat = d3['heights']['lat_ph'][subset1]
28.     lon = d3['heights']['lon_ph'][subset1]
29.     classed_pc_flag = d3['classed_pc_flag'][subset1]
31.     return x_act, h, signal_conf_ph, lat, lon, classed_pc_flag
34. def get_atl03_data(filepath, beam):
35.     """
36.     读取ATL03数据，根据维度截取数据
37.     Args:
38.         filepath (str): h5文件路径
39.         beam (str): 光束
40.     Returns:
41.         返回沿轨道距离，高程距离，光子置信度
42.     """
43.     atl03_file = glob(filepath)
44.     is2_atl03_mds = read_hdf5_atl03_beam_h5py(atl03_file[0], beam=beam, verbose=False)
45.     # 添加沿轨道距离到数据中
46.     get_atl03_x_atc(is2_atl03_mds)
47.     return is2_atl03_mds
50. def show_classification(x_origin, y_origin, classification, clz):
51.     """
52.     :param clz: -1:未分类, 0:噪声, 1:地形, 2:冠层, 3:冠顶, 4:海洋
53.     :param classification: 分类数据
54.     :param y_origin:
55.     :param x_origin:
56.     """
57.     plt.subplots(num=1, figsize=(24, 6))
58.     ax = plt.gca()
59.     ax.get_xaxis().get_major_formatter().set_useOffset(False)
60.     plt.xticks(rotation=270)
61.     ax.set_xlabel('x_atc, km')
62.     ax.set_ylabel('h, m')
63.     ax.xaxis.set_major_locator(MultipleLocator(100))
64.     colors = ['red', 'black', 'green', 'violet', 'blue', 'grey']
65.     for flag in clz:
66.         idx = np.where(classification == flag)
67.         plt.scatter(x_origin[idx], y_origin[idx], s=5, c=colors[flag])
69.     plt.show()
72. if __name__ == '__main__':
73.     data = {
74.         'filepath': 'D:\\Users\\SongW\\Documents\\ICESat-2 Data\\ATL03\\ATL03_20200620024106_13070701_005_01.h5',
75.         'filepath_08': 'D:\\Users\\SongW\\Documents\\ICESat-2 Data\\ATL08\\ATL08_20200620024106_13070701_005_01.h5',
76.         'beam': 'gt2l',
77.         'mask': [19.6468, 19.6521]
78.     }
79.     atl03_data = atl03_data = get_atl03_data(data['filepath'], data['beam'])
80.     add_atl08_classed_flag(data['filepath_08'], data['beam'], atl03_data)
82.     x_origin, y_origin, conf, lat, lon, classed_pc_flag = select_atl03_data(atl03_data, data['mask'])
84.     show_classification(x_origin, y_origin, classed_pc_flag, [-1, 0, 1, 2, 3])

复制代码

项目源码

sx-code - icesat-2-atl03 (github.com)

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。