ICESat-2 ATL03數(shù)據(jù)和ATL08數(shù)據(jù)的分段距離不一致，ATL08在ATL03的基礎(chǔ)上重新分段，并對分段內(nèi)的數(shù)據(jù)做處理得到一系列的結(jié)果，詳情見數(shù)據(jù)字典：

ATL08 Product Data Dictionary (nsidc.org)

ATL08使用DRAGANN算法對ATL03數(shù)據(jù)做了去噪處理，并使用分類算法對每個光子進行分類

ATL08使用ph_segment_id和classed_pc_indx可以和ATL03對應(yīng)起來。基于此，可從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('請?zhí)钊胝_的光束代碼')
        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分類數(shù)據(jù)到ATL03中
    Args:
        filepath_08: ATL08數(shù)據(jù)文件位置
        beam: 波束，與ATL03保持一致
        atl03_mod: ATL03數(shù)據(jù)

    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']

    # -- 分段中的第一個光子（轉(zhuǎn)換為基于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內(nèi)的特定光子。
    classed_pc_index = val_08['signal_photons']['classed_pc_indx']
    # 每個光子的陸地植被ATBD分類標(biāo)志為噪聲、地面、樹冠和樹冠頂部。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數(shù)據(jù)代碼見：http://www.rzrgm.cn/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):
    """
    選擇數(shù)據(jù)范圍
    Args:
        atl03_data: 所有數(shù)據(jù)
        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數(shù)據(jù)，根據(jù)維度截取數(shù)據(jù)
    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)
    # 添加沿軌道距離到數(shù)據(jù)中
    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: 分類數(shù)據(jù)
    :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)