process_lidar_data.py

import numpy as np
import laspy
from itertools import groupby

lazData_path = ["server/example_lidar_data/ot_N_000005_1.laz", "server/example_lidar_data/ot_N_000033_1.laz"]

# Info about data
with laspy.open(lazData_path[0]) as fh:
    # Print metadata properties
    print("File Version:", fh.header.version)
    print("Point Count:", fh.header.point_count)
    print("Scale Factors:", fh.header.scale)
    print("Offset:", fh.header.offset)

    print('Points from Header:', fh.header.point_count)
    las = fh.read()
    print(las)
    print('Points from data:', len(las.points))
    ground_pts = las.classification == 2
    bins, counts = np.unique(las.return_number[ground_pts], return_counts=True)
    print('Ground Point Return Number distribution:')
    for r, c in zip(bins, counts):
        print('    {}:{}'.format(r, c))

# Limit the data to specific areas
# should be within this range (10.616913,60.742712) - (10.825653,60.940206)
# this is only temporary data for coordinates in arctic
areas = [
    [(-3671212, 7898422), (-3200175, 4699978)],
    [(-3671212, 7898422), (-3200175, 4699978)]
]

def inArea(position, areaRange):
    x, y, _ = position
    if (areaRange[0][0] < x < areaRange[1][0]) and (areaRange[0][1] > y > areaRange[1][1]):
        return True
    else:
        return False

# Calculation of height gotten from las files
def find_height(points):
    sorted_coords = sorted(points, key=lambda coord: (coord[0],coord[1]))

    groupCoords = [list(group) for key, group in groupby(sorted_coords, key=lambda coord: (coord[0], coord[1]))]

    height_differences = []
    groupss = []

    for groups in groupCoords:
        if len(groups) < 2 or groups[0] == groups[1]:
            continue

        min_height = min(coords[2] for coords in groups)
        max_height = max(coords[2] for coords in groups)

        difference = max_height - min_height
        height_differences.append(difference)
        groupss.append(groups)

    return min(height_differences)

print(lazData_path[0])

# Refactor data format
iceOver = laspy.read(lazData_path[0])
iceUnder = laspy.read(lazData_path[1])

ice_points = list(zip(iceOver.X,iceOver.Y,iceOver.Z)) + list(zip(iceUnder.X,iceUnder.Y,iceUnder.Z))

area_heights = []

# areas
for area in areas:
    ice_points = list(filter(lambda position: inArea(position, area), ice_points))
    area_heights.append(find_height(ice_points))

print(area_heights)