process_lidar_data.py

import json
import math
import os

import laspy
import numpy as np

from server.data_processing.area_processing import calculate_corners, define_gridareas, inArea, find_height, \
    closest_points, position_relative_to_pointcloud, define_grid_lidardata


# Info about data
def about_laz_file(path):
    """
    write info about the lidar's raw data

    :param path: path to the lidar data
    :return: list with info about version, point_count, scale in meters and offset in meters
    """
    with laspy.open(path) 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)

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

    return [las.header.version, las.header.point_count, las.header.scale, las.header.offset]

def find_folder_files(direcory):
    """
    find all files in a directory

    :param direcory: path to directory
    :return: list with all file paths within directory
    """
    files = []
    for root, _, fileNames in os.walk(direcory):
        for fileName in fileNames:
            files.append(os.path.join(root, fileName))
    return files

# find the height of an area based on the coordinates of it's center
# and it's affiliations (subId and groupId) (soon to be implemented
def calculate_area_data(center, body_of_water, path):
    """
    calculate area

    :param center: center point of area to be calculated
    :param body_of_water: body of water area belongs to
    :param path: path to lidar data
    :return: the calculated ice thickness in designated area
    """
    # container for all the heights in area
    area_heights = []

    # zone coords taken
    taken_coords = []

    # read json data with path to data for specific water body
    file_name = "server/map_handler/lake_relations/" + body_of_water + "_div.json"
    with open(file_name) as data:
        map_data = json.load(data)

    # grid cell offset
    # NB: make these a global var
    grid_size = 4
    #cell_x, cell_y = map_data['tile_width'], map_data['tile_height']
    #cell_x, cell_y = (400*grid_size, 400*grid_size)
    cell_x, cell_y = (2000, 1000)

    # convert the offset to meter
    #cell_x = 111.320 * cell_x * 1000
    #cell_y = (40075 * math.cos(cell_y) * 1000) / 360
    cell_x = 111.320 * cell_x
    cell_y = (111.320 * math.cos(60)*cell_y)

    # set the limit of the area compared to local coordinates
    area_limit = calculate_corners(center[0], center[1], (cell_x, cell_y))
    #area_limit = calculate_corners(sub_center[0], sub_center[1], (cell_x/4, cell_y/4))

    # grid data
    map_data = map_data['features']
    map_zones = [area_limit[2], area_limit[0]]
    map_data = list(filter(lambda point_position: inArea((point_position['properties']['sub_div_center'][0], point_position['properties']['sub_div_center'][1], 0.0), map_zones), map_data))

    # Refactor lidar data to a readable format
    iceOver = laspy.read(path)

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

    max_point = max(ice_points)
    min_point = min(ice_points)

    # define all the sub-areas within the area, local coordinates
    grid_sub_area = define_gridareas(center[0], center[1], (cell_x, cell_y),grid_size)
    # define all the sub-areas within the area, lidar coordinates
    grid_area_lidar_heights = define_grid_lidardata((min_point, max_point), grid_size, ice_points)

    sub_area_heights = list(zip(grid_sub_area, grid_area_lidar_heights))

    # find the heights of each sub-area => area-heights
    if len(map_data) > 0:
        for sub_area in sub_area_heights:
            start = (sub_area[0][2])
            end = (sub_area[0][0])

            # calculate map zones height
            ys, xs = start
            ye, xe = end
            sub_center = ((xs + xe)/2, (ys + ye)/2)
            # check if area is part of water body
            part_of_subarea_of_waterbody = list(filter(lambda pos: inArea((pos['properties']['sub_div_center'][0], pos['properties']['sub_div_center'][1], 0.0), [start,end]), map_data))

            if(len(part_of_subarea_of_waterbody) > 0):
                current_map_zone = closest_points(sub_center, part_of_subarea_of_waterbody, taken_coords)
                sub_center = current_map_zone['properties']['sub_div_center']
                taken_coords.append(sub_center)
            else:
                print("sub area not found on map")
                continue

            current_zone_id = current_map_zone['properties']['sub_div_id']

            if current_map_zone is not None:
                #                     sub_id        center         heights
                area_heights.append((current_zone_id, sub_center, sub_area[1]))
        return area_heights
    else:
        return [] # return [0] if no data collected from lidar


# print(calculate_area_data((61, 11), 'mj\u00f8sa', "server\\lidar_data\\mj\u00f8sa\\measurement_id_2.laz"))