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): 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): files = [] for root, _, fileNames in os.walk(direcory): for fileName in fileNames: files.append(os.path.join(root, fileName)) return files print(find_folder_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): # test dataset laz_root_path = "server\\lidar_data\\" + body_of_water laz_data_paths = find_folder_files(laz_root_path) # 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(laz_data_paths[0]) iceUnder = laspy.read(laz_data_paths[1]) # add two files together temporary test data(soon to be removed) ice_points = list(zip(iceOver.X,iceOver.Y,iceOver.Z)) + list(zip(iceUnder.X,iceUnder.Y,iceUnder.Z)) print(len(ice_points)) ice_points = list(filter(lambda point_position: inArea((point_position[0],point_position[1],0.0),[(-3300000, 4500000), (-3200000, 5000000)]), ice_points)) print(len(ice_points)) # only for visualizating the testing data #print("max",max(ice_points)) #print("min",min(ice_points)) #x = [points[0] for points in ice_points] #y = [points[1] for points in ice_points] # Plot the points #plt.plot(x, y, 'bo') # 'bo' specifies blue color and circle markers #plt.xlabel('X axis') #plt.ylabel('Y axis') #plt.title('Plotting Points') #plt.grid(True) # Add grid #plt.show() 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) print('grid size: ', grid_sub_area) print('grid size: ', grid_area_lidar_heights) sub_area_heights = list(zip(grid_sub_area, grid_area_lidar_heights)) print('sub_area_heights: ', sub_area_heights[0]) # 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]) #test data # zone coordinates sett to be relative to the lidar data's point cloud #print("l1 = ", min_point, " l2 = ", max_point) #print("p1 = ", start, " p2 = ", end) #print("center_lim = ", center) #print("max_lim = ", area_limit) #areazone = position_relative_to_pointcloud(min_point, max_point, start, end, center, area_limit) #print("area",areazone) # 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)) print('sub_area: ', sub_area) print('part_of_subarea_of_waterbody: ', part_of_subarea_of_waterbody) 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) print("item added", sub_center, " len ", len(taken_coords)) else: print("sub area not found on map") continue current_zone_id = current_map_zone['properties']['sub_div_id'] # filter data within sub-area zones #print(areazone[0][0]," < ",ice_points[0][0]," < ",areazone[1][0]) #print(areazone[0][1]," > ",ice_points[0][1]," > ",areazone[1][1]) #points_in_area = list(filter(lambda point_position: inArea(point_position, areazone), ice_points)) 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'))