Skip to content
Snippets Groups Projects
Commit 74888a04 authored by Hoa Ben The Nguyen's avatar Hoa Ben The Nguyen
Browse files

Merge branch 'connect_data_to_map' of gitlab.stud.idi.ntnu.no:sarasdj/prog2900

parents 6579433c 4fb33e9a
No related branches found
No related tags found
Loading
Showing
with 310 additions and 108 deletions
import math
from itertools import groupby
from math import pi, cos
EARTH = 6378.137 # Radius of the earth in kilometer
METER = (1 / ((2 * pi / 360) * EARTH)) / 1000 # 1 meter in degree
# find real world coordinates to lidar
#
# p - area using real world coordinate
# l - area using lidar coordinate(pointcloud)
# max_limit - maximum area limit of a lidar scann in real world coordinates
def position_relative_to_pointcloud(l1, l2, p1, p2, center_lim, max_lim):
center_l = tuple((a+b)/2 for a, b in zip(l1,l2))
return [
((p1[0] - center_lim[0]) * ((l1[0] - l2[0]) / (max_lim[2][0] - max_lim[0][0])) + center_l[0],
(p1[1] - center_lim[1]) * ((l1[1] - l2[1]) / (max_lim[2][1] - max_lim[0][1])) + center_l[1]),
((p2[0] - center_lim[0]) * ((l1[0] - l2[0]) / (max_lim[2][0] - max_lim[0][0])) + center_l[0],
(p2[1] - center_lim[1]) * ((l1[1] - l2[1]) / (max_lim[2][1] - max_lim[0][1])) + center_l[1]),
]
#print(position_relative_to_pointcloud((-20,20), (-10,30), (1,3), (2,2), (2.5,2.5), [(5,5),(0,5),(0,0),(5,0)]))
#print(position_relative_to_pointcloud((-3299999, 4608018), (-3200001, 4687153), (61.47620866851029, 8.961138281887507), (61.95241733702057, 6.8508373935926645), (61, 11), [(61.95241733702057, 15.125349549679886), (60.04758266297943, 15.125349549679886), (60.04758266297943, 6.8746504503201145), (61.95241733702057, 6.8746504503201145)]))
# check if lidar points is within range of the area selected
def inArea(position, areaRange):
x, y, _ = position # position to be checked
#print((areaRange[0][0])," < ",x," < ",(areaRange[1][0])," and ",(areaRange[0][1])," < ",(y)," < ",(areaRange[1][1])," ",((areaRange[0][0]) < x < (areaRange[1][0])) and ((areaRange[0][1]) < (y) < (areaRange[1][1])))
if ((areaRange[0][0]) < x < (areaRange[1][0])) and ((areaRange[0][1]) < y < (areaRange[1][1])):
return True
else:
return False
# find distance between two points
def distance(point1, point2):
y1, x1 = point1
y2, x2 = point2
return math.sqrt(abs(y2 - y1)**2 + abs(x2 - x1)**2)
# find the closest point in json list
def closest_points(point, list, coords_taken):
closest_point_found = None
closest_dist = float('inf')
for current_point in list:
dist = distance(point, current_point['properties']['sub_div_center'])
if dist < closest_dist and current_point not in coords_taken:
closest_dist = dist
closest_point_found = current_point
return closest_point_found
# Calculation of height in an area from las files
def find_height(points):
if len(points) == 0:
return [0]
height_differences = [] # final container for height data
# sort the points
sorted_coords = sorted(points, key=lambda coord: (coord[0], coord[1]))
# group the sorted points that has the same xy- coordinates together
groupCoords = [list(group) for key, group in groupby(sorted_coords, key=lambda coord: (coord[0], coord[1]))]
# loop through the groups to find the difference in height
for group in groupCoords:
if len(group) < 2 or group[0] == group[1]:
continue # jump over iteration if there is only one coordinate or lidar registered the same point twice
# find max and min height
min_height = min(coords[2] for coords in group)
max_height = max(coords[2] for coords in group)
# difference between them
difference = max_height - min_height
height_differences.append(difference)
# list of thickness in an area
return height_differences
def calculate_corners(lat, lng, area_offset):
"""Calculate corners of polygon based on a center coordinate
......@@ -12,11 +84,12 @@ def calculate_corners(lat, lng, area_offset):
"""
# From https://stackoverflow.com/questions/7477003/calculating-new-longitude-latitude-from-old-n-meters
# Formulas:
lat_pos = lat + (area_offset * METER)
lng_pos = lng + (area_offset * METER) / cos(lat * (pi / 180))
offset_lng, offset_lat = area_offset
lat_pos = lat + (offset_lat * METER)
lng_pos = lng + (offset_lng * METER) / cos(lat * (pi / 180))
lat_neg = lat - (area_offset * METER)
lng_neg = lng - (area_offset * METER) / cos(lat * (pi / 180))
lat_neg = lat - (offset_lat * METER)
lng_neg = lng - (offset_lng * METER) / cos(lat * (pi / 180))
return [
(lat_neg, lng_pos), # top left
......@@ -32,37 +105,64 @@ def define_gridareas(lat, lng, area_offset, grid_size):
# find the main area's corner positions
main_area = calculate_corners(lat, lng, area_offset)
offset_lng, offset_lat = area_offset
# find the main area's range size
area_size = (main_area[2][0] - main_area[0][0], main_area[0][1] - main_area[2][1])
# find each subareas vector to it's center
dist_to_subcenter = (area_size[0]/(grid_size*2), area_size[1]/(grid_size*2))
# find subareas size relative to main area
subarea_offset = area_offset/grid_size
subarea_offset_lng = offset_lng/grid_size
subarea_offset_lat = offset_lat/grid_size
group_id = 0
i=0
# find the center coordinates of each area in grid to find the corner areas
for y in (range(grid_size)):
relative_size_lng = y / grid_size
relative_size_lat = y / grid_size
for x in (range(grid_size)):
relative_size_lat = x / grid_size
relative_size_lng = x / grid_size
lat_pos = main_area[3][0] + relative_size_lat * area_size[0] + dist_to_subcenter[0]
lng_pos = main_area[3][1] + relative_size_lng * area_size[1] + dist_to_subcenter[1]
# id of each subarea
sub_id = x % grid_size + grid_size * y
if (x+(y*grid_size)) % ((grid_size*2) + 1) == grid_size * 2:
group_id += 2
if (x+(y*grid_size)) % 2 == 0 and (x+(y*grid_size)) != 0:
i += 1
# use the center of sub areas to find the corner of each subarea
subarea_offset = (subarea_offset_lng, subarea_offset_lat)
corners = calculate_corners(lat_pos, lng_pos, subarea_offset)
grided_area.append((sub_id, group_id + i%2, corners))
grided_area.append(corners)
return grided_area
#print(define_gridareas(-60,10,10000,4))
# separate the zones into smaller area, into a grid
def define_grid_lidardata(max_area, grid_size, points):
# container for an area turned into a grid of areas
grided_area = []
bottom_left, top_right = max_area
# find subareas size relative to the full size of data
subarea_offset_lng = (top_right[0]-bottom_left[0])/grid_size
subarea_offset_lat = (top_right[1]-bottom_left[1])/grid_size
# find the center coordinates of each area in grid to find the corner areas
for y in (range(grid_size)):
relative_size_lat = y / grid_size
for x in (range(grid_size)):
relative_size_lng = x / grid_size
bottom_left_range = (subarea_offset_lng * (relative_size_lng - 1 / 2) + bottom_left[0],
subarea_offset_lat * (relative_size_lat - 1 / 2) + bottom_left[1])
top_right_range = (subarea_offset_lng * relative_size_lng + bottom_left[0],
subarea_offset_lat * relative_size_lat + bottom_left[1])
area_range = [bottom_left_range, top_right_range]
area = list(filter(lambda point_position: inArea(point_position, area_range), points))
height_in_area = find_height(area)
grided_area.append(height_in_area)
return grided_area
#[1,2,2,3,4,5,6,3,4,6,8,9,5,3,5.7,8,5,3]
#print(define_gridareas(60,10,1500,4))
#print(define_gridareas(3435693.5,6299200.0, 400000,4))
\ No newline at end of file
import numpy as np
import json
import math
import laspy
import utm # src: https://github.com/Turbo87/utm
from itertools import groupby
from server.data_processing.area_processing import calculate_corners, define_gridareas
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
# hard coded files for test data
# hard coded files for test datas
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:
las = fh.read()
ground_pts = las.classification == 2
bins, counts = np.unique(las.return_number[ground_pts], return_counts=True)
# check if lidar points is within range of the area selected
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 in an area from las files
def find_height(points):
height_differences = [] # final container for height data
# sort the points
sorted_coords = sorted(points, key=lambda coord: (coord[0],coord[1]))
# group the sorted points that has the same xy- coordinates together
groupCoords = [list(group) for key, group in groupby(sorted_coords, key=lambda coord: (coord[0], coord[1]))]
# loop through the groups to find the difference in height
for group in groupCoords:
if len(group) < 2 or group[0] == group[1]:
continue # jump over iteration if there is only one coordinate or lidar registered the same point twice
# find max and min height
min_height = min(coords[2] for coords in group)
max_height = max(coords[2] for coords in group)
# difference between them
difference = max_height - min_height
height_differences.append(difference)
# list of thickness in an area
return height_differences
def about_laz_file():
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)
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]
# 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):
def calculate_area_data(center, body_of_water):
# 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], 1500)
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(lazData_path[0])
......@@ -63,29 +72,82 @@ def calculate_area_data(center):
# 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_area_heights = define_gridareas(60,10, 1500,4)
# find the heights of each sub-area => area-heights
for sub_area in grid_area_heights:
start = min(sub_area[2])
end = max(sub_area[2])
#test data
# zone coordinates sett to be relative to the lidar data
areazone = [(((min_point[0] - max_point[0]) * ((start[0]-center[0])/(area_limit[1][0]-area_limit[3][0]))) + (min_point[0] - max_point[0])/2 + max_point[0],
(((min_point[1] - max_point[1]) * ((start[1]-center[1])/(area_limit[1][1]-area_limit[3][1]))) + (min_point[1] - max_point[1])/2 + min_point[1])),
((((min_point[0] - max_point[0]) * ((end[0]-center[0])/(area_limit[1][0]-area_limit[3][0]))) + (min_point[0] - max_point[0])/2 + max_point[0],
(((min_point[1] - max_point[1]) * ((end[1]-center[1])/(area_limit[1][1]-area_limit[3][1])))) + (min_point[1] - max_point[1])/2 + min_point[1]))]
# filter data within sub-area zones
points_in_area = list(filter(lambda point_position: inArea(point_position, areazone), ice_points))
area_heights.append((sub_area[0],sub_area[1],find_height(points_in_area)))
# 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)
return area_heights
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'))
......@@ -119,7 +119,7 @@ class IceHTTP(BaseHTTPRequestHandler):
def do_POST(self):
if self.path == '/new_lidar_data':
input_new_Lidar_data(self, self.cursor, 1, 'Mjosa') # hardcoded body of water must change later
input_new_Lidar_data(self, self.cursor, 1, 'mj\u00f8sa') # hardcoded body of water must change later
# Start a server on port 8443 using self defined HTTP class
......
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
import json
import os
from datetime import datetime
from server.data_processing.process_lidar_data import calculate_area_data
from server.data_processing.process_lidar_data import calculate_area_data, about_laz_file
# input_new_Lidar_data send new data gathered from the lidar and send it to the database (from the drone, most likely)
def input_new_Lidar_data(self, cursor, sensorId, bodyOfWater):
try:
print("name=",bodyOfWater)
# hard coded coordinates
latitude = 60.0
longitude = 10.0
latitude = 60.816848
longitude = 10.723823
total_measurement_average = 0 # the total average of a measurement
# data about the file read from
about_laz = about_laz_file()
scale_factor = max(about_laz[2])
# create a new measurement with the time the data is sent, sensor type, where
# and an estimate of average thickness of ice on water body
cursor.execute('''
INSERT INTO Measurement( SensorID, TimeMeasured, WaterBodyName,
WholeAverageThickness, CenterLat, CenterLon) VALUES
INSERT INTO Measurement(MeasurementID, SensorID, TimeMeasured, WaterBodyName, CenterLat, CenterLon) VALUES
(?,?,?,?,?,?);
''', (sensorId, datetime.utcnow().replace(microsecond=0), bodyOfWater, 0, latitude, longitude))
''', (1, sensorId, datetime.utcnow().replace(microsecond=0), bodyOfWater, latitude, longitude))
# auto generate new measurement id
measurement_id = cursor.lastrowid
# calculate the area of to be calculated based on the coordinates given to the calculation model
areas_data = calculate_area_data((latitude, longitude))
areas_data = calculate_area_data((latitude, longitude), bodyOfWater)
lidar_json_data = []
if (areas_data):
if(areas_data):
# store lidar data in jason formate
# calculate data for each zone within the area
for area in areas_data:
if (len(area[2]) != 0):
average = sum(area[2]) / len(area[2])
minimum_thickness = min(area[2])
# lng and lat relative to map
subId = int(area[0])
map_lat, map_lng = area[1]
heights = area[2]
if(len(heights) != 0 or sum(heights) == 0):
average = sum(heights)/len(heights)
minimum_thickness = min(heights)
else:
average = 0
minimum_thickness = 0
total_measurement_average += average
time_measured = datetime.utcnow().replace(microsecond=0)
# input the data into the database
cursor.execute('''
INSERT INTO SubDivision(MeasurementID, SubDivisionID, GroupID, MinimumThickness, AverageThickness, CenterLatitude, CenterLongitude, Accuracy) VALUES
INSERT INTO SubdivisionMeasurementData(MeasurementID, TimeMeasured, SubdivID, WaterBodyName, MinimumThickness, AverageThickness, CalculatedSafety, Accuracy) VALUES
(?,?,?,?,?,?,?,?);
''', (measurement_id, area[0], area[1], float(minimum_thickness), float(average), float(latitude),
float(longitude), float(1)))
''', (measurement_id, time_measured, subId, bodyOfWater, float(minimum_thickness), float(average), float(0.0), scale_factor))
sub_center = (map_lat, map_lng)
# set up json formate
lidar_read = {
'MeasurementId': str(measurement_id),
'SubId': str(subId),
'SubCenter': str(sub_center),
'Heights': str(heights)
}
total_measurement_average = total_measurement_average / len(areas_data)
lidar_json_data.append(lidar_read)
# input the newly generated measurement_id and whole average thickness
cursor.execute('''
UPDATE Measurement
SET measurementID = ?, WholeAverageThickness = ?
WHERE MeasurementID IS NULL AND WholeAverageThickness = 0;
''', (int(measurement_id), total_measurement_average), )
SET measurementID = ?
WHERE MeasurementID IS NULL;
''', (int(measurement_id),))
else:
print('No data found')
print('No data found, line 79')
# send the changes to the database
cursor.connection.commit()
# Send response
self.send_response(200)
self.send_header('Content-type', "application/json")
self.end_headers()
file_path = "./server/map_handler/lake_relations/"+bodyOfWater+"_lidar_data.json"
content = None
current_directory = os.getcwd()
print("Current working directory:", current_directory)
if len(lidar_json_data) > 0:
if os.path.exists(file_path):
os.remove(file_path)
# convert list of lidar data to json
content = json.dumps(lidar_json_data)
with open(file_path, "w") as file:
file.write(content)
else:
print('No data found, line 101')
content = json.dumps([])
# Write content data to response object
self.wfile.write(content.encode('utf-8'))
# error handling
except Exception as e:
print("An error occurred", e)
print("An error occurred: ", e)
# rollback in case of error
cursor.connection.rollback()
[{"MeasurementId": "7", "SubId": "36", "SubCenter": "(60.841532, 10.717878)", "Heights": "[1]"}, {"MeasurementId": "7", "SubId": "83", "SubCenter": "(60.828326, 10.982563)", "Heights": "[1, 27]"}, {"MeasurementId": "7", "SubId": "33", "SubCenter": "(60.771059, 10.698341)", "Heights": "[1]"}, {"MeasurementId": "7", "SubId": "136", "SubCenter": "(60.396856, 11.220933)", "Heights": "[1, 4, 6, 27, 7]"}]
\ No newline at end of file
[{"MeasurementId": "17", "SubId": "48", "SubCenter": "(60.816856, 10.723847)", "Heights": "[1]"}, {"MeasurementId": "17", "SubId": "100", "SubCenter": "(60.816856, 10.990419)", "Heights": "[1, 27]"}, {"MeasurementId": "17", "SubId": "36", "SubCenter": "(60.756856, 10.710419)", "Heights": "[1]"}, {"MeasurementId": "17", "SubId": "168", "SubCenter": "(60.396856, 11.230419)", "Heights": "[1, 4, 6, 27, 7]"}]
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment