ICESat-2 ATL03光子数据读取

这篇具有很好参考价值的文章主要介绍了ICESat-2 ATL03光子数据读取。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

ICESat-2数据处理的方式一般为将光子数据投影到沿轨距离和高程的二维空间。如下图:

ICESat-2 ATL03光子数据读取

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()

输出图像如下:

ICESat-2 ATL03光子数据读取

项目源码

sx-code - icesat-2-atl03 (github.com)文章来源地址https://www.toymoban.com/news/detail-859924.html

到了这里,关于ICESat-2 ATL03光子数据读取的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 【数据处理】Pandas读取CSV文件示例及常用方法(入门)

    查看读取前10行数据 2067 向前填充 指定列的插值填充 使用某数据填充指定列的空值 示例: 类似切片 array([‘SE’, ‘cv’, ‘NW’, ‘NE’], dtype=object) 类似数据库查询中的groupby查询 先添加新的一列按月将数据划分 聚合,对指定的列按月划分求平均值等 min 最小值 max 最大值 sum

    2024年02月06日
    浏览(205)
  • Python点云处理(一)点云数据读取与写入

    当处理点云数据时,我们通常需要读取各种不同格式的点云文件。Python作为一种强大的编程语言,在点云处理领域提供了许多库和工具,可以帮助我们读取和处理各种格式的点云文件。本文将介绍如何使用Python读取和写入各种格式的点云文件。 LAS(Lidar Data Exchange)和LAZ(L

    2024年02月08日
    浏览(38)
  • Python处理xlsx文件(读取、转为列表、新建、写入数据、保存)

    xlsxwriter**库对于xslx表的列数不做限制, xlrd 库不能写入超过65535行,256列的数据。 由于需要处理的数据行列数较多,遇到报错才发现库的限制问题,记录一下。

    2024年02月12日
    浏览(66)
  • USSOCOM Urban3D 数据集读取与处理

    Urban3D数据集图像为正摄RGB影像,分辨率为50cm。 从SpaceNet上使用aws下载数据,文件夹结构为: 每一块.tif大小为2048*2048。 使用深度学习需要对Urban3D数据进行裁剪,这里采用 torchvision.transforms.RandomCrop 进行裁剪。 RandomCrop 可以直接作用于PIL.Image打开的文件和torch类型的数据上,但

    2024年01月23日
    浏览(31)
  • 【记录】USSOCOM Urban3D 数据集读取与处理

    Urban3D数据集图像为正摄RGB影像,分辨率为50cm。 从SpaceNet上使用aws下载数据,文件夹结构为: 每一块.tif大小为2048*2048。 使用深度学习需要对Urban3D数据进行裁剪,这里采用 torchvision.transforms.RandomCrop 进行裁剪。 RandomCrop 可以直接作用于PIL.Image打开的文件和torch类型的数据上,但

    2024年02月11日
    浏览(33)
  • Python两种读取txt与csv文件方式(利用numpy处理数据)

    一共80个数据(只截取前10个数据) 在excel中显示的内容 在pycharm中显示的内容 一共80个数据 在记事本中显示的内容 在pycharm中显示的内容 1、读取所有内容 data_pd打印结果 2、数据转为numpy data_np打印结果 1、读取所有内容 data_pd打印结果 2、数据转为numpy data_np打印结果 1、读取所

    2023年04月11日
    浏览(54)
  • NCDC气象数据的提取与处理(四):python批量读取、写入nc数据经纬度格点数值

    1.问题描述: 2.思路: 3.实现过程: 3.1格点位置匹配 3.2写入表格 4.运行效果 4.1打包站点信息 4.2读取nc文件列表 4.3提取对应格点的nc数据 4.4数据写入 NCDC的站点数据处理在之前三节里已经介绍过了,但是NCDC的就那么几种数据可能不能满足日常使用,比如说辐射数据他就没有。

    2024年02月05日
    浏览(73)
  • 【Matlab】如何读取文件夹下所有txt数据进行处理并以txt结果更名输出

    如何读取文件夹下所有txt数据进行处理并以txt结果更名输出 目录 前言 一、Matlab中fullfile函数用法 二、使用步骤 1.读取文件夹下所有txt文件并以struct存储变量 2.循环下读取每个txt文件中的数据并进行处理 总结 遇到Matlab需要大批量处理一个文件夹下所有的txt格式,经过信号处

    2024年02月07日
    浏览(68)
  • 【完整教程】在win10平台下使用d435i深度相机读取数据并保存到本地,以便进行后续3D点云处理

    进入网址:RealSense SDK 2.0 直接拉到网站最下端,在Asset下可以看到很多exe可执行软件,由于我的电脑是win10,所以选择第三个。说句题外话,鄙人曾经考英语六级时记得Asset专门指不动资产,没错,就是房子! 下载完成后文件夹内有如下图所示软件,直接安装即可。 安装完成

    2024年02月02日
    浏览(174)
  • SAP_ABAP_编程基础_文件处理(CRUD)_R3系统_打开文件 / 关闭文件 / 删除文件 / 向文件中写入数据 / 从文件中读取数据 / 使用服务器上的文件

    SAP ABAP 顾问(开发工程师)能力模型_Terry谈企业数字化的博客-CSDN博客 文章浏览阅读490次。目标:基于对SAP abap 顾问能力模型的梳理,给一年左右经验的abaper 快速成长为三年经验提供超级燃料! https://blog.csdn.net/java_zhong1990/article/details/132469977 平时在  ‘ 工地搬砖 ’,很少关

    2024年02月22日
    浏览(46)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包