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
import os
import h5py
import re
def read_hdf5_atl08(filename, beam, verbose=False):
file_id = h5py.File(os.path.expanduser(filename), 'r')
# 输出HDF5文件信息
if verbose:
print(file_id.filename)
print(list(file_id.keys()))
print(list(file_id['METADATA'].keys()))
# 为ICESat-2 ATL08变量和属性分配python字典
atl08_mds = {}
# 读取文件中每个输入光束
beams = [k for k in file_id.keys() if bool(re.match('gt\\d[lr]', k))]
if beam not in beams:
print('请填入正确的光束代码')
return
atl08_mds['signal_photons'] = {}
# -- ICESat-2 Geolocation Group
for key, val in file_id[beam]['signal_photons'].items():
atl08_mds['signal_photons'][key] = val[:]
return atl08_mds
映射ATL08
将 ATL08 映射到 ATL03
def get_atl08_mapping(atl03_ph_index_beg, atl03_segment_id, atl08_classed_pc_indx,
atl08_classed_pc_flag, atl08_segment_id):
"""
Function to map ATL08 to ATL03 class photons
Args:
atl03_ph_index_beg:
atl03_segment_id:
atl08_classed_pc_indx:
atl08_classed_pc_flag:
atl08_segment_id:
Returns:
"""
# Get ATL03 data
indsNotZero = atl03_ph_index_beg != 0
atl03_ph_index_beg = atl03_ph_index_beg[indsNotZero]
atl03_segment_id = atl03_segment_id[indsNotZero]
# Find ATL08 segments that have ATL03 segments
atl03SegsIn08TF, atl03SegsIn08Inds = ismember(atl08_segment_id, atl03_segment_id)
# Get ATL08 classed indices and values
atl08classed_inds = atl08_classed_pc_indx[atl03SegsIn08TF]
atl08classed_vals = atl08_classed_pc_flag[atl03SegsIn08TF]
# Determine new mapping into ATL03 data
atl03_ph_beg_inds = atl03SegsIn08Inds
atl03_ph_beg_val = atl03_ph_index_beg[atl03_ph_beg_inds]
newMapping = atl08classed_inds + atl03_ph_beg_val - 2
# Get max size of output array
sizeOutput = newMapping[-1]
# Pre-populate all photon classed array with zeroes
allph_classed = (np.zeros(sizeOutput + 1)) - 1
# Populate all photon classed array from ATL08 classifications
allph_classed[newMapping] = atl08classed_vals
# Return all photon classed array
return allph_classed
添加分类信息
def add_atl08_classed_flag(filepath_08, beam, atl03_mod):
"""
添加ATL08分类数据到ATL03中
Args:
filepath_08: ATL08数据文件位置
beam: 波束，与ATL03保持一致
atl03_mod: ATL03数据
Returns:
携带ATL08分类信息
"""
val_03 = atl03_mod
val_08 = read_hdf5_atl08(filepath_08, beam)
# val_03['classed_pc_flag'] = np.zeros_like(val_03['heights']['h_ph']) + np.NaN
atl03_heights = val_03['heights']['h_ph']
# -- 分段中的第一个光子（转换为基于0的索引）
segment_index_begin = val_03['geolocation']['ph_index_beg']
segment_id = val_03['geolocation']['segment_id']
# 追踪到ATL03上特定20m Segment_ID的光子的段ID
ph_segment_id = val_08['signal_photons']['ph_segment_id']
# 该索引追溯到ATL03上20m segment_id内的特定光子。
classed_pc_index = val_08['signal_photons']['classed_pc_indx']
# 每个光子的陆地植被ATBD分类标志为噪声、地面、树冠和树冠顶部。0=噪音，1=地面，2=冠层，或3=冠层顶部
classed_pc_flag = val_08['signal_photons']['classed_pc_flag']
# Map ATL08 classifications to ATL03 Photons
all_ph_classed = get_atl08_mapping(segment_index_begin, segment_id,
classed_pc_index, classed_pc_flag, ph_segment_id)
if len(all_ph_classed) < len(atl03_heights):
n_zeros = len(atl03_heights) - len(all_ph_classed)
zeros = np.zeros(n_zeros)
all_ph_classed = np.append(all_ph_classed, zeros)
val_03['classed_pc_flag'] = all_ph_classed
使用姿势
读取ATL03数据代码见：https://www.cnblogs.com/sw-code/p/18161987
from glob import glob
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.ticker import MultipleLocator
from readers.add_atl08_info import add_atl08_classed_flag
from readers.get_ATL03_x_atc import get_atl03_x_atc
from readers.read_HDF5_ATL03 import read_hdf5_atl03_beam_h5py
def select_atl03_data(atl03_data, mask):
"""
选择数据范围
Args:
atl03_data: 所有数据
mask (list): 维度范围
Returns:
"""
# 选择范围
d3 = atl03_data
subset1 = (d3['heights']['lat_ph'] > min(mask)) & (d3['heights']['lat_ph'] < max(mask))
x_act = d3['heights']['x_atc'][subset1]
h = d3['heights']['h_ph'][subset1]
signal_conf_ph = d3['heights']['signal_conf_ph'][subset1]
lat = d3['heights']['lat_ph'][subset1]
lon = d3['heights']['lon_ph'][subset1]
classed_pc_flag = d3['classed_pc_flag'][subset1]
return x_act, h, signal_conf_ph, lat, lon, classed_pc_flag
def get_atl03_data(filepath, beam):
"""
读取ATL03数据，根据维度截取数据
Args:
filepath (str): h5文件路径
beam (str): 光束
Returns:
返回沿轨道距离，高程距离，光子置信度
"""
atl03_file = glob(filepath)
is2_atl03_mds = read_hdf5_atl03_beam_h5py(atl03_file[0], beam=beam, verbose=False)
# 添加沿轨道距离到数据中
get_atl03_x_atc(is2_atl03_mds)
return is2_atl03_mds
def show_classification(x_origin, y_origin, classification, clz):
"""
:param clz: -1:未分类, 0:噪声, 1:地形, 2:冠层, 3:冠顶, 4:海洋
:param classification: 分类数据
:param y_origin:
:param x_origin:
"""
plt.subplots(num=1, figsize=(24, 6))
ax = plt.gca()
ax.get_xaxis().get_major_formatter().set_useOffset(False)
plt.xticks(rotation=270)
ax.set_xlabel('x_atc, km')
ax.set_ylabel('h, m')
ax.xaxis.set_major_locator(MultipleLocator(100))
colors = ['red', 'black', 'green', 'violet', 'blue', 'grey']
for flag in clz:
idx = np.where(classification == flag)
plt.scatter(x_origin[idx], y_origin[idx], s=5, c=colors[flag])
plt.show()
if __name__ == '__main__':
data = {
'filepath': 'D:\\Users\\SongW\\Documents\\ICESat-2 Data\\ATL03\\ATL03_20200620024106_13070701_005_01.h5',
'filepath_08': 'D:\\Users\\SongW\\Documents\\ICESat-2 Data\\ATL08\\ATL08_20200620024106_13070701_005_01.h5',
'beam': 'gt2l',
'mask': [19.6468, 19.6521]
}
atl03_data = atl03_data = get_atl03_data(data['filepath'], data['beam'])
add_atl08_classed_flag(data['filepath_08'], data['beam'], atl03_data)
x_origin, y_origin, conf, lat, lon, classed_pc_flag = select_atl03_data(atl03_data, data['mask'])
show_classification(x_origin, y_origin, classed_pc_flag, [-1, 0, 1, 2, 3])
项目源码
sx-code - icesat-2-atl03 (github.com)
本文来自博客园，作者：sw-code，转载请注明原文链接：https://www.cnblogs.com/sw-code/p/18161991