area_processing.py

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


# check if lidar points is within range of the area selected
def inArea(position, areaRange):
    """
    finds out if position is in area's range

    :param position: current checked position in
    :param areaRange: range of checked position, containing a touple with max and min limit
    :return: True if position is in range otherwise False
    """
    x, y, _ = position # position to be checked
    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):
    """
    calculates the distance between two points

    :param point1: first point
    :param point2: second points
    :return: the distance between two points
    """
    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):
    """
    finds the closest points within a list of points and returns the closest point that has not been taken

    :param point: position to check for closest points
    :param list: data about each section of a lake map, content used = center of cells
    :param coords_taken: list of coordinates that has been taken
    :return: the closest to point that has not been 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['properties']['sub_div_center'] 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):
    """
    finds the heights of a list of points that has grouped together points with identical XY coordinates

    :param points: a list of points in an area
    :return: the heights within a list of 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 z-position
        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

    :param lat -- center latitude
    :param lng -- center longitude
    :return: a list of the corners of the polygon based on a center coordinate
    """
    # From https://stackoverflow.com/questions/7477003/calculating-new-longitude-latitude-from-old-n-meters
    # Formulas:
    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 - (offset_lat * METER)
    lng_neg = lng - (offset_lng * METER) / cos(lat * (pi / 180))

    return [
        (lat_neg, lng_pos), # top left
        (lat_pos, lng_pos), # top right
        (lat_pos, lng_neg), # bottom right
        (lat_neg, lng_neg)  # bottom left
    ]

# separate the zones into smaller area, into a grid
def define_gridareas(lat, lng, area_offset, grid_size):
    """
    Define the grid areas of a zones center and offset

    :param lat: center latitude of the grid
    :param lng: center longitude of the grid
    :param area_offset: offset between center to corner
    :param grid_size: number of zone divisions into grid
    :return: list of corners of each grid, based on the number of grid divisions and offset of the zone
    """
    # container for an area turned into a grid of areas
    grided_area = []

    # find the main area's corner positions
    main_area = calculate_corners(lat, lng, area_offset)

    # 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_lng = offset_lng/grid_size
    subarea_offset_lat = offset_lat/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 # relative y position on grid
        for x in (range(grid_size)):
            relative_size_lng = x / grid_size # relative x position on grid
            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]

            # 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(corners)

    return grided_area

# separate the zones into smaller area, into a grid
def define_grid_lidardata(max_area, grid_size, points):
    """
    Define the grid of the lidar data and divide each point to its designated area,
    based on maximum limited area, grid size and points position on the grid

    :param max_area: corner coordinates of its maximum area
    :param grid_size: number of area divisions into grids
    :param points: list of points
    :return: list of grid with all the points inside its designated grid
    """
    # container for an area turned into a grid of areas
    grided_area = []

    # smallest and biggest corner in 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])
    subarea_offset_lat = (top_right[1]-bottom_left[1])

    # find subareas size relative to the full size of data
    subarea_hypotenuse_lng = subarea_offset_lng / grid_size
    subarea_hypotenuse_lat = subarea_offset_lat / grid_size

    # organize all coordinates into their corresponding area
    # their own 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_sub_corner = (subarea_offset_lng * relative_size_lng + bottom_left[0],
                                 subarea_offset_lat * relative_size_lat + bottom_left[1])


            top_right_sub_corner = (subarea_offset_lng * relative_size_lng + bottom_left[0] + subarea_hypotenuse_lng,
                               subarea_offset_lat * relative_size_lat + bottom_left[1] + subarea_hypotenuse_lat)

            # smallest and highest points in area
            area_range = [bottom_left_sub_corner, top_right_sub_corner]

            # filter out points that is outside of area range
            area = list(filter(lambda point_position: inArea(point_position, area_range), points))

            # find the ice thickness from all the points in area
            thickness_in_area = find_height(area)

            # store all ice thickness in area
            grided_area.append(thickness_in_area)

    return grided_area