AgRS-FE/c_Extract_PC_Feature.py at main · George-Horus/AgRS-FE · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
import laspy
import numpy as np
import pandas as pd
from scipy.stats import skew, kurtosis

def extract_pointcloud_stats_by_coords(
    las_path,
    coord_excel_path,
    output_excel_path="pointcloud_stats_by_location.xlsx",
    radius=1.0,
    x_col="x",
    y_col="y"
):
    """
    根据坐标表提取.las点云文件指定邻域的统计特征，并输出为Excel。
    :param las_path: 点云.las文件路径
    :param coord_excel_path: Excel文件路径，需有x/y或lon/lat列
    :param output_excel_path: 输出统计结果的Excel文件名
    :param radius: 搜索半径（单位与点云坐标一致，通常为米）
    :param x_col: Excel中x坐标列名
    :param y_col: Excel中y坐标列名
    """
    print("读取点云文件...")
    las = laspy.read(las_path)
    x = las.x
    y = las.y
    z = las.z

    print("读取Excel坐标表...")
    excel = pd.read_excel(coord_excel_path)
    if {x_col, y_col}.issubset(excel.columns):
        coords = excel[[x_col, y_col]].to_numpy()
    else:
        raise ValueError(f'Excel必须包含 "{x_col}, {y_col}" 列')

    results = []
    print(f"共{len(coords)}个坐标点，开始逐点统计...")
    for idx, (x0, y0) in enumerate(coords):
        # 计算每个点到当前坐标点的距离，筛选在半径范围内的点
        dist = np.sqrt((x - x0) ** 2 + (y - y0) ** 2)
        mask = dist <= radius
        z_sub = z[mask]
        if len(z_sub) < 1:
            # 区域内无点，全部置为NaN或0
            stats = dict(
                max_height=np.nan, mean_height=np.nan, std_height=np.nan,
                iqr_height=np.nan, p25=np.nan, p75=np.nan,
                skewness=np.nan, kurtosis=np.nan, point_density=0
            )
        else:
            # 区域面积 = π*r²
            area = np.pi * radius ** 2
            point_density = len(z_sub) / area
            stats = dict(
                max_height=np.max(z_sub),
                mean_height=np.mean(z_sub),
                std_height=np.std(z_sub),
                iqr_height=np.percentile(z_sub, 75) - np.percentile(z_sub, 25),
                p25=np.percentile(z_sub, 25),
                p75=np.percentile(z_sub, 75),
                skewness=skew(z_sub),
                kurtosis=kurtosis(z_sub),
                point_density=point_density
            )
        # 添加原始坐标信息
        stats.update({x_col: x0, y_col: y0})
        results.append(stats)
        if (idx + 1) % 10 == 0 or (idx + 1) == len(coords):
            print(f"已完成 {idx + 1}/{len(coords)}")

    # 保存统计结果为Excel
    df_out = pd.DataFrame(results)
    df_out.to_excel(output_excel_path, index=False)
    print(f"统计结果已保存到 {output_excel_path}")

if __name__ == "__main__":

    extract_pointcloud_stats_by_coords(
        las_path="exp_data/0317_PC.las",
        coord_excel_path="exp_data/0317_coordinate.xlsx",
        output_excel_path="c_point_cloud.xlsx",
        radius=2.0,
        x_col="lon",
        y_col="lat"
    )