ICESat-2 ATL03光子数据读取
ICESat-2数据处理的方式一般为将光子数据投影到沿轨隔断和高程的二维空间。如下图：https://img2024.cnblogs.com/blog/2207877/202404/2207877-20240427135657124-96037553.png
ATL03数据读取
H5是一种数据存储布局，读取原理就是按照该布局获取数据，这里给出两种读取方式。
ATL03的数据字典：ATL03 Product Data Dictionary (nsidc.org)
使用pandas
import warnings
import pandas as pd
def read_hdf5_atl03_beam_pandas(filename, beam, verbose=False):
# 打开HDF5文件进行读取
h5_store = pd.HDFStore(filename, mode='r')
root = h5_store.root
# 为ICESat-2 ATL03变量和属性分配python字典
atl03_mds = {}
# 读取文件中每个输入光束
# beams = ', k))]
beams = ['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']
if beam not in beams:
print('请填入正确的光束代码')
return
atl03_mds['heights'] = {}
atl03_mds['geolocation'] = {}
# -- 获取每个HDF5变量
with warnings.catch_warnings():
warnings.simplefilter("ignore")
# -- ICESat-2 Heights Group
heights_keys = ['dist_ph_across', 'dist_ph_along', 'h_ph', 'lat_ph', 'lon_ph', 'signal_conf_ph']
for key in heights_keys:
atl03_mds['heights'] = root['heights'][:]
geolocation_keys = ['ref_elev', 'ph_index_beg', 'segment_id', 'segment_ph_cnt', 'segment_dist_x', 'segment_length']
# -- ICESat-2 Geolocation Group
for key in geolocation_keys:
atl03_mds['geolocation'] = root['geolocation'][:]
h5_store.close()
return atl03_mds使用h5py
import os
import h5py
import re
def read_hdf5_atl03_beam_h5py(filename, beam, verbose=False):
"""
ATL03 原始数据读取
Args:
filename (str): h5文件路径
beam (str): 光束
verbose (bool): 输出HDF5信息
Returns:
返回ATL03光子数据的heights和geolocation信息
"""
# 打开HDF5文件进行读取
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 ATL03变量和属性分配python字典
atl03_mds = {}
# 读取文件中每个输入光束
beams = ', k))]
if beam not in beams:
print('请填入正确的光束代码')
return
atl03_mds['heights'] = {}
atl03_mds['geolocation'] = {}
atl03_mds['bckgrd_atlas'] = {}
# -- 获取每个HDF5变量
# -- ICESat-2 Measurement Group
for key, val in file_id['heights'].items():
atl03_mds['heights'] = val[:]
# -- ICESat-2 Geolocation Group
for key, val in file_id['geolocation'].items():
atl03_mds['geolocation'] = val[:]
for key, val in file_id['bckgrd_atlas'].items():
atl03_mds['bckgrd_atlas'] = val[:]
return atl03_mds重建沿轨道隔断
在读取ICESat-2 ATL03数据后，我们必要根据分段信息重建每个光子的沿轨道隔断，已知数据如下：
[*]ATL03每个分段内的光子都有一个沿轨道隔断（dist_ph_along），该隔断始于当前分段。
[*]ATL03每个分段有一个沿轨道隔断（segment_dist_x），该隔断始于赤道交叉口，即真正的沿轨道隔断。
[*]当前分段内每个光子的相对沿轨道隔断 + 当前分段的沿轨道隔断 = 每个光子的真正沿轨道隔断。
代码如下：
import numpy as np
def get_atl03_x_atc(atl03_mds):
val = atl03_mds
# 初始化
val['heights']['x_atc'] = np.zeros_like(val['heights']['h_ph']) + np.NaN
val['heights']['y_atc'] = np.zeros_like(val['heights']['h_ph']) + np.NaN
val['geolocation']['ref_elev_all'] = np.zeros_like(val['heights']['h_ph'])
# -- ATL03 Segment ID
segment_id = val['geolocation']['segment_id']
# -- 分段中的第一个光子（转换为基于0的索引）
segment_index_begin = val['geolocation']['ph_index_beg'] - 1
# -- 分段中的光子事件数
segment_pe_count = val['geolocation']['segment_ph_cnt']
# -- 每个ATL03段的沿轨道距离
segment_distance = val['geolocation']['segment_dist_x']
# -- 每个ATL03段的轨道长度
segment_length = val['geolocation']['segment_length']
# -- 对ATL03段进行迭代，以计算40m的平均值
# -- 在ATL03中基于1的索引：无效==0
# -- 此处为基于0的索引：无效==-1
segment_indices, = np.nonzero((segment_index_begin[:-1] >= 0) &
(segment_index_begin >= 0))
for j in segment_indices:
# -- j 段索引
idx = segment_index_begin
# -- 分段中的光子数（使用2个ATL03分段）
c1 = np.copy(segment_pe_count)
c2 = np.copy(segment_pe_count)
cnt = c1 + c2
# -- 沿轨道和跨轨道距离
# -- 获取当前段光子列表，idx当前段(j)第一个光子数量，c1当前段光子数量，idx+c1当前段长度
distance_along_x = np.copy(val['heights']['dist_ph_along'])
ref_elev = np.copy(val['geolocation']['ref_elev'])
# -- 给当前段的光子加上当前段沿轨道距离
distance_along_x[:c1] += segment_distance
distance_along_x += segment_distance
distance_along_y = np.copy(val['heights']['dist_ph_across'])
val['heights']['x_atc'] = distance_along_x
val['heights']['y_atc'] = distance_along_y
val['geolocation']['ref_elev_all'] += ref_elevATL03数据截取
在处理ATL03时，我们一般都会获取经过研究区域内的光子数据，因此必要对数据举行截取操纵，代码如下：
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, 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) elif mask_lon is None and mask_lon is not None: subset1 = subset1 & (d3['heights']['x_atc'] = min(mask_lon)) & (d3['heights']['x_atc'] = min(mask_lat)) & (d3['lat'] = min(mask_lon)) & (d3['lon']
页:
[1]