【ICESat-2处理】ATL03的去噪（python）

我最近整理了一点ICESat-2的处理代码，我看现在网上的分享内容还挺少的，所以就想把自己整理的发出来吧，方便大家用。我是一个粗人，我不怎么会写代码，不喜欢只分享处理的关键函数，所以我的代码都是从数据读取到最后的可视化都有，直接拿去运行的。请读者们根据自己的需求挑选吧。

前言

我不懂原理，我也不想讲原理，总之这一节的处理就是把ICESat-2 ATL03的光子数据的噪点去除，只留下其中密度最大的一群光子点，我试了很多种方法，我最推荐密度域值法，好用！

ps：白天的数据噪点真的特别多，数据量还特别大，能用晚上的还是用晚上的数据吧。

ATL03的去噪代码

运行前，设置代码中的参数

包括输入和输出文件名路径：

【一定注意】输入是之前讲过的，ATL03读取之后的CSV文件啊！如下所示：

我这里有五个字段，分别是dist_ph_along（沿轨距离），dist_ph_across（跨轨距离），h_ph（高度），lat_ph，lon_ph（纬度和经度）

还包括两个处理函数（基于统计的去噪和基于半径的去噪）的域值

输出结果：去噪后的CSV和一副去噪后的示意图，CSV会保留输入数据的字段不变，就是把噪点行剔除了

完整代码

import pandas as pd
import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt

# 【提示】噪声和数据的时间也有关系，最好下载晚上的数据，白天的数据噪声多且数据量大（因为都是噪点）
# 该代码在处理后的新csv会仅存去噪后的点云，噪点被删去

inputCSVpath = "E:/ICEsat2/ATL03_L08-6/outputtest/L03_55.csv"
outputCSVpath = "E:/ICEsat2/ATL03_L08-6/outputtest/L03_55_inliers_radius.csv"

def statistical_outlier_removal(pcd, neighbors=20, ratio=2.0):
    cl, ind = pcd.remove_statistical_outlier(nb_neighbors=neighbors, std_ratio=ratio)
    inlier_cloud = pcd.select_by_index(ind)
    outlier_cloud = pcd.select_by_index(ind, invert=True)
    return ind, outlier_cloud, inlier_cloud

def radius_outlier_removal(pcd, nb_points=16, radius=10):
    cl, ind = pcd.remove_radius_outlier(nb_points=nb_points, radius=radius)
    inlier_cloud = pcd.select_by_index(ind)
    outlier_cloud = pcd.select_by_index(ind, invert=True)
    return ind, outlier_cloud, inlier_cloud

# 读取ATL03 CSV数据
atlo3_data = pd.read_csv(inputCSVpath)

# 将数据转换为点云
point_cloud = o3d.geometry.PointCloud()
#这里横纵坐标也可以变
point_cloud.points = o3d.utility.Vector3dVector(atlo3_data[['dist_ph_along','dist_ph_across', 'h_ph']].values)


# 【基于半径阈值的方法，后面还有基于统计的方法】
# 对点云进行半径异常点移除
inlier_indices_radius, outlier_cloud_radius, inlier_cloud_radius = radius_outlier_removal(point_cloud)

# 提取半径异常点移除后的内点对应的行
inlier_data_radius = atlo3_data.iloc[inlier_indices_radius]

# 将内点数据保存为新的CSV文件
inlier_data_radius.to_csv(outputCSVpath, index=False)

# 可视化原始点云和去噪后的点云
# 转换为numpy数组以便可视化
outlier_points_radius = np.asarray(outlier_cloud_radius.points)
inlier_points_radius = np.asarray(inlier_cloud_radius.points)

plt.figure(figsize=(10, 8))

# 绘制半径异常点移除结果
plt.scatter(outlier_points_radius[:, 0], outlier_points_radius[:, 2], c='gray', label='Radius Outliers', s=0.5)
plt.scatter(inlier_points_radius[:, 0], inlier_points_radius[:, 2], c='red', label='Radius Inliers', s=1)

# 设置图例
plt.legend()

# 设置坐标轴标签
plt.xlabel('Distance along the track(m)')
plt.ylabel('height(m)')

# 显示图形
plt.show()


# 【基于统计的方法，我感觉在噪声多的时候没有基于半径的方法好用】
# 对点云进行统计异常点移除
# inlier_indices_statistical, outlier_cloud_statistical, inlier_cloud_statistical = statistical_outlier_removal(point_cloud)
# 提取统计异常点移除后的内点对应的行
# inlier_data_statistical = atlo3_data.iloc[inlier_indices_statistical]
# 转换为numpy数组以便可视化
# outlier_points_statistical = np.asarray(outlier_cloud_statistical.points)
# inlier_points_statistical = np.asarray(inlier_cloud_statistical.points)
# 将内点数据保存为新的CSV文件
# inlier_points_statistical.to_csv("D:/ICEsat2/ATL03_L08-6/L03_55_inliers_statistical.csv", index=False)
# plt.figure(figsize=(10, 8))
# 绘制统计异常点移除结果
# plt.scatter(outlier_points_statistical[:, 0], outlier_points_statistical[:, 2], c='red', label='Statistical Outliers', s=1)
# plt.scatter(inlier_points_statistical[:, 0], inlier_points_statistical[:, 2], c='gray', label='Statistical Inliers', s=1)
# 设置图例
# plt.legend()
#
# # 设置坐标轴标签
# plt.xlabel('Distance along the track(m)')
# plt.ylabel('height(m)')
#
# # 显示图形
# plt.show()