ICESat-2数据处理的方式一般为将光子数据投影到沿轨距离和高程的二维空间。如下图:
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 for k in file_id.keys() if bool(re.match('gt\\d[lr]', 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'][key] = root[beam]['heights'][key][:]
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'][key] = root[beam]['geolocation'][key][:]
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 for k in file_id.keys() if bool(re.match('gt\\d[lr]', 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[beam]['heights'].items():
atl03_mds['heights'][key] = val[:]
# -- ICESat-2 Geolocation Group
for key, val in file_id[beam]['geolocation'].items():
atl03_mds['geolocation'][key] = val[:]
for key, val in file_id[beam]['bckgrd_atlas'].items():
atl03_mds['bckgrd_atlas'][key] = 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[1:] >= 0))
for j in segment_indices:
# -- j 段索引
idx = segment_index_begin[j]
# -- 分段中的光子数(使用2个ATL03分段)
c1 = np.copy(segment_pe_count[j])
c2 = np.copy(segment_pe_count[j + 1])
cnt = c1 + c2
# -- 沿轨道和跨轨道距离
# -- 获取当前段光子列表,idx当前段(j)第一个光子数量,c1当前段光子数量,idx+c1当前段长度
distance_along_x = np.copy(val['heights']['dist_ph_along'][idx: idx + cnt])
ref_elev = np.copy(val['geolocation']['ref_elev'][j])
# -- 给当前段的光子加上当前段沿轨道距离
distance_along_x[:c1] += segment_distance[j]
distance_along_x[c1:] += segment_distance[j + 1]
distance_along_y = np.copy(val['heights']['dist_ph_across'][idx: idx + cnt])
val['heights']['x_atc'][idx: idx + cnt] = distance_along_x
val['heights']['y_atc'][idx: idx + cnt] = distance_along_y
val['geolocation']['ref_elev_all'][idx: idx + c1] += ref_elev
ATL03数据截取
在处理ATL03时,我们一般都会获取经过研究区域内的光子数据,因此需要对数据进行截取操作,代码如下:
from glob import glob
from readers.get_ATL03_x_atc import get_atl03_x_atc
from readers.read_HDF5_ATL03 import read_hdf5_atl03_beam, read_hdf5_atl03_coordinate
def 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'] <= max(mask_lat))
if mask_lon is not None:
if mask_lon[0] is not None and mask_lon[1] is None:
subset1 = subset1 & (d3['heights']['x_atc'] >= mask_lon[0])
elif mask_lon[0] is None and mask_lon[1] is not None:
subset1 = subset1 & (d3['heights']['x_atc'] <= mask_lon[1])
else:
subset1 = subset1 & (d3['heights']['x_atc'] >= min(mask_lon)) & (d3['heights']['x_atc'] <= max(mask_lon))
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]
ref_elev = d3['geolocation']['ref_elev_all'][subset1]
del d3, subset1
return x_act, h, signal_conf_ph, lat, lon, ref_elev
def read_all_beam_coordinate(filepath, mask_lat, mask_lon):
"""
读取所有波束的数据
:param filepath:
:param mask_lat:
:param mask_lon:
:return:
"""
atl03_file = glob(filepath)
is2_atl03_mds = read_hdf5_atl03_coordinate(atl03_file[0])
# 禁止加载全部数据
# if mask_lat is None or len(mask_lat) == 0 or mask_lon is None or len(mask_lon) == 0:
# return False
d3 = is2_atl03_mds
if mask_lon is None and mask_lat is None:
# 加载全部数据
return d3
for beam in is2_atl03_mds.keys():
subset1 = (d3[beam]['lat'] >= min(mask_lat)) & (d3[beam]['lat'] <= max(mask_lat))
subset1 = subset1 & (d3[beam]['lon'] >= min(mask_lon)) & (d3[beam]['lon'] <= max(mask_lon))
d3[beam]['lat'] = d3[beam]['lat'][subset1]
d3[beam]['lon'] = d3[beam]['lon'][subset1]
return d3
数据可视化
使用沿轨道距离和高程数据绘制散点图,示例代码如下:
def save2file(act, h, conf, lat, lon):
"""
保存研究区域的一下数据
:param act: act,沿轨道距离
:param h: h,高程
:param conf: 置信度
:param lat: 维度
:param lon: 经度
"""
points = list(zip(act, h, lat, lon, conf))
data = pd.DataFrame(points, columns=['沿轨道距离', '高程', '维度', '经度', '置信度'])
data.to_csv('result/points_origin.csv', mode='w', index=False)
if __name__ == '__main__':
filepath = r'D:\Users\SongW\Downloads\ATL03_20190222135159_08570207_005_01.h5'
beam = 'gt3l'
mask_lat = [16.533, 16.550]
act, h, conf, lat, lon, ref_elev = read_data(filepath, beam, mask_lat, None)
save2file(act, h, conf, lat, lon)
plt.scatter(act, h)
plt.show()
输出图像如下:
文章来源:https://www.toymoban.com/news/detail-859924.html
项目源码
sx-code - icesat-2-atl03 (github.com)文章来源地址https://www.toymoban.com/news/detail-859924.html
到了这里,关于ICESat-2 ATL03光子数据读取的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!