【Lidar】Laspy库介绍+基础函数（读取、可视化、保存、旋转、筛选、创建点云数据）

def read_las(file=r'G:/彭俊喜/10.las'):  # 读取点云数据
    las = laspy.read(file)
    las_header = las.header
    # point = las.points  # print(point)
    # points_data = np.v stack((las.x, las.y, las.z, las.intensity, las.gps_time, las.point_source_id,
    # las.classification, las.user_data, las.red, las.green, las.blue)).transpose()
    # xyz = points_data[:, :3]  # 以数组的形式显示所有点
    # print(points)
    print('主版本号：' + str(las_header.major_version))
    print('副版本号：' + str(las_header.minor_version))
    print('最小值：%f,%f,%f' % (las_header.min[0], las_header.min[1], las_header.min[2]))
    print('最大值：%f,%f,%f' % (las_header.max[0], las_header.max[1], las_header.max[2]))
    print('比例：%f,%f,%f' % (las_header.scale[0], las_header.scale[1], las_header.scale[2]))
    print('偏移量：%f,%f,%f' % (las_header.offset[0], las_header.offset[1], las_header.offset[2]))
    print('点云数量:%d' % las_header.point_records_count)
    # 遍历点
    for i in range(las_header.point_records_count):
        # 遍历所有的点云数据
        print('x=%f, y=%f, z=%f, intensity=%d, GPStime=%f, PointSourceID=%f, '
              'Classification=%d, UserData=%d, Red=%d, Green=%d, Blue=%d'
              % (las.x[i], las.y[i], las.z[i],
                 las.intensity[i], las.gps_time[i], las.point_source_id[i],
                 las.classification[i], las.user_data[i],
                 las.red[i], las.green[i], las.blue[i]))

def filtering_las(file=r'G:/彭俊喜/10.las'):  # 通过条件筛选点云数据（颜色/回波强度等）
    las = laspy.read(file)  # 读取点云数据
    new_file_c0 = laspy.create(point_format=las.header.point_format,
                               file_version=las.header.version)

    def multiple_conditions(value1, value2, value3, value4):
        # 创建筛选条件
        condition1 = value1 < 60
        condition2 = value1 > 105
        condition3 = value2 < 125
        condition4 = value2 > 255
        condition5 = value3 < 120
        condition6 = value3 > 255
        condition7 = value4 > 100
        # 若条件都满足，返回True，否则返回False
        return ((condition1 or condition2) and (condition3 or condition4) and (condition5 or condition6)) and condition7
    result = np.vectorize(multiple_conditions)(las.red, las.green, las.blue, las.intensity)  # 获取符合条件的布尔值
    new_file_c0.points = las.points[result]  # 通过布尔值选择点云数据
    new_file_c0.write(r'G:/彭俊喜/10.las')  # 写入点云数据

def rotation_las(file=r'G:/彭俊喜/10.las'):  # 旋转点云
    las = laspy.read(file)  # 读取点云数据
    las_header = las.header
    center_x = (las_header.min[0] + las_header.max[0]) / 2
    center_y = (las_header.min[1] + las_header.max[1]) / 2  # 计算中心点（有的点云含有地理坐标）
    cos = -1  # 控制旋转角度
    sin = math.sqrt(1 - cos ** 2)
    x_all = []
    y_all = []
    xj = []  # 创建新列表存储旋转后的点云数据
    for i in range(las_header.point_records_count):
        # 遍历所有的点云数据
        x = las.x[i] - center_x
        y = las.y[i] - center_y
        x_rotation = x * cos - y * sin + center_x
        y_rotation = x * sin + y * cos + center_y
        x_all.append(x_rotation)
        y_all.append(y_rotation)
        xj.append(i)  # 将旋转后的点加入列表中
    header = laspy.LasHeader(point_format=3, version="1.2")  # 定义头文件
    header.offsets = [np.min(x_all), np.min(y_all), las.header.offsets[2]]  # 定义偏移量
    header.scales = np.array([0.0001, 0.0001, 0.0001])  # 定义缩放比例
    header.add_extra_dim(laspy.ExtraBytesParams(name="xj", type=np.int32))
    new_las = laspy.LasData(header)
    new_las.x = x_all
    new_las.y = y_all
    new_las.z = las.z
    new_las.intensity = las.intensity
    new_las.gps_time = las.gps_time
    new_las.point_source_id = las.point_source_id
    new_las.classification = las.classification
    new_las.red = las.red
    new_las.green = las.green
    new_las.blue = las.blue
    new_las.xj = xj
    new_las.write('G:/彭俊喜/10.las')

def create_las(file=r'G:/彭俊喜/10.las'):  # 创建点云
    # 0. Creating some dummy data
    my_data_xx, my_data_yy = np.meshgrid(np.linspace(-20, 20, 15), np.linspace(-20, 20, 15))
    my_data_zz = my_data_xx ** 2 + 0.25 * my_data_yy ** 2
    my_data = np.hstack((my_data_xx.reshape((-1, 1)), my_data_yy.reshape((-1, 1)), my_data_zz.reshape((-1, 1))))

    # 1. Create a new header
    header = laspy.LasHeader(point_format=3, version="1.2")
    header.add_extra_dim(laspy.ExtraBytesParams(name="random", type=np.int32))
    header.offsets = np.min(my_data, axis=0)
    header.scales = np.array([0.1, 0.1, 0.1])

    # 2. Create a Las
    las = laspy.LasData(header)

    las.x = my_data[:, 0]
    las.y = my_data[:, 1]
    las.z = my_data[:, 2]
    las.random = np.random.randint(-1503, 6546, len(las.points), np.int32)

    las.write(file)

def save_point(path, point_type):  # 保存点云
    # ---------------------------laspy库保存----------------------------
    las = laspy.read(r"Z:\Personal\彭俊喜\Lidar_try/2.las")  # read a las file
    points = las.points
    out_file = laspy.LasData(las.header)
    ground = [2, 5, 6]  # 索引
    out_file.points = points[np.array(ground)]  # extract ground points, and save it to a las file.
    out_file.write(r"Z:\Personal\彭俊喜\Lidar_try/out1.las")

def show_point(file=r"Z:\Personal\彭俊喜\Lidar_try/out.las"):  # 可视化点云
    las = laspy.read(file)  # read a las file
    points = las.points
    xyz = np.vstack((las.x, las.y, las.z)).transpose()  # extract x, y, z and put into a list  # 点云的空间位置
    x = las.x
    y = las.y
    z = las.z
    # 绘制三维散点图
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')
    ax.set_xlabel('X Axes')
    ax.set_ylabel('Y Axes')
    ax.set_zlabel('Z Axes')
    # ax.grid(None)  # 隐藏格网
    ax.scatter(x, y, z)
    plt.show()

# -*- coding: utf-8 -*-
"""
@Time ： 2023/11/30 15:30
@Auth ： RS迷途小书童
@File ：Laspy Function.py
@IDE ：PyCharm
@Purpose：laspy库的基础操作
"""
import math
import laspy
import numpy as np
import matplotlib.pyplot as plt