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