Functions for finding points and distances along multi-segment lines.

profile_tools.py

'''Functions for finding points and distances along multi-segment lines.

All angles are mathematical convention.

'''
import numpy as np


def distance(point1, point2):
    '''Return distance between *point1* and *point2*.'''
    return np.sqrt((point2[0] - point1[0])**2 + (point2[1] - point1[1])**2)

    
def bearing_func(point1, point2):
    return np.arctan2(point2[1] - point1[1], point2[0] - point1[0])
    
    
normal_unit_vector = lambda point, angle: [np.sin(angle), -1 * np.cos(angle)]


def projected_point(point, angle, distance):
    '''Return the point lying at some *distance* at some *angle* from a *point*.'''
    return [point[0] + distance * np.cos(angle),
            point[1] + distance * np.sin(angle)]

            
def angle_to_azimuth(angle):
    '''Convert math *angle* (in degrees) to an azimuth (in degrees).'''
    return 90 - np.rad2deg(angle)

    
def azimuth_to_angle(angle):
    '''Convert azimuth (in degrees) to a math *angle* (in degrees).'''
    return 90 - np.rad2deg(angle)


def get_distance_locations(line_vertices, spacing=100, debug=0):
    '''Return a list of (*distances, locations*) for points spaced *spacing*
    units apart along the line described by *line_vertices*.'''
    distances = np.arange(0, total_distance(line_vertices) + spacing, spacing)
    return get_distance_location_for_distances(distances, line_vertices, debug=debug)
    
    
def get_distance_location_for_distances(distances, line_vertices, debug=0):
    '''Return a list of (*distances, locations*) for points at *distances*
    along the line described by *line_vertices*.'''
    verts = line_vertices
    locations = [verts[0]]
    distances = [0]
    current_bearing = bearing_func(verts[0], verts[1])
    last_passed_vertex_index = 0
    while last_passed_vertex_index + 1 < len(verts):
        last_passed_vertex = verts[last_passed_vertex_index]
        next_unpassed_vertex = verts[last_passed_vertex_index + 1]
        candidate_location = projected_point(locations[-1], current_bearing, spacing)
        potential_segment_distance = distance(last_passed_vertex, next_unpassed_vertex)
        candidate_location_distance = distance(last_passed_vertex, candidate_location)
        if candidate_location_distance > potential_segment_distance:
            temp_moved = distance(locations[-1], next_unpassed_vertex)
            temp_to_move = spacing - temp_moved
            if len(verts) > last_passed_vertex_index + 2:
                new_bearing = bearing_func(next_unpassed_vertex, verts[last_passed_vertex_index + 2])
            else:
                new_bearing = current_bearing
            new_bearing_candidate_location = projected_point(next_unpassed_vertex, new_bearing, temp_to_move)
            distances.append(distances[-1] + temp_moved + temp_to_move)
            locations.append(tuple(new_bearing_candidate_location))
            current_bearing = new_bearing
            last_passed_vertex_index += 1
        else:
            locations.append(tuple(candidate_location))
            distances.append(distances[-1] + spacing)
    return distances, locations


def distance_to_line_segment(point, line_segment_vertices, debug=0):
    '''Return the distance from *point* to the line segment described by
    *line_segment_vertices*.'''
    vertex1 = line_segment_vertices[0]
    vertex2 = line_segment_vertices[1]
    segment_length = distance(vertex2, vertex1)
    if debug: print 'distance_to_line_segment point', point, 'segment', line_segment_vertices, 'segment_length', segment_length
    if segment_length == 0:
        d = distance(vertex1, point)
        if debug: print '\tsegment length zero,', d
        return d
    t = ((point[0] - vertex1[0]) * (vertex2[0] - vertex1[0]) + (point[1] - vertex1[1]) * (vertex2[1] - vertex1[1])) / (segment_length ** 2)
    if debug: print '\tt', t
    if t < 0:
        d = distance(vertex1, point)
        if debug: print '\tpoint past ', vertex1, 'vertex1, returning', d
        return d
    elif t > 1:
        d = distance(vertex2, point)
        if debug: print '\tpoint past ', vertex2, 'vertex2, returning', d
        return d
    else:
        d = distance((vertex1[0] + t * (vertex2[0] - vertex1[0]),
                      vertex1[1] + t * (vertex2[1] - vertex1[1])), point)
        if debug: print '\tpoint is perpendicular to line segment, returning', d
        return d
    
            
def projected_distance_on_line_segment(point, line_segment_vertices):
    '''Return the distance along the line segment described by 
    *line_segment_vertices* to the point which is closest to *point*.'''
    vertex1 = line_segment_vertices[0]
    vertex2 = line_segment_vertices[1]
    d = distance_to_line_segment(point, line_segment_vertices)
    vx, vy = normal_unit_vector(point, bearing_func(vertex1, vertex2))
    dx = (point[0] + vx * d) - vertex1[0]
    dy = (point[1] + vy * d) - vertex1[1]
    return np.sqrt(dx ** 2 + dy ** 2)
    
    
def line_segments(line_vertices):
    '''Return a list of two vertices for each segment in the line described by
    *line_vertices*.
    
    e.g. 
    
        >>> line_segments([[0, 0], [0, 1], [1, 5], [6, 6]])
        [[[0, 0], [0, 1]],
         [[0, 1], [1, 5]],
         [[1, 5], [6, 6]]]
         
    '''
    segment_vertices = []
    for i in range(1, len(line_vertices)):
        segment_vertices.append(line_vertices[i-1: i+1])
    return segment_vertices

    
def projected_distance_on_line(point, line_vertices, debug=False):
    '''Return the distance along the line to the point on the profile described
    by *line_vertices* that is closest to *point*.'''
    if debug: print 'point', point
    candidate_distances = []
    segment_vertices_all = line_segments(line_vertices)
    if debug: print 'segment_vertices_all', segment_vertices_all
    for segment_vertices in segment_vertices_all:
        candidate_distances.append(distance_to_line_segment(point, segment_vertices, debug=debug))
    if debug: print 'candidate_distances', candidate_distances
    j = np.argsort(candidate_distances)[0]
    cumulative_distance = 0
    for i in range(j + 1):
        segment_vertices = segment_vertices_all[i]
        if debug: print 'loop 2, i', i, 'segment_vertices', segment_vertices
        if i < j:
            if debug: print 'loop 2, i', i, 'adding full segment length', distance(*segment_vertices)
            cumulative_distance += distance(*segment_vertices)
        elif i == j:
            p = projected_distance_on_line_segment(point, segment_vertices)
            if debug: print 'loop 2, i', i, 'j', j, 'adding projected length', p
            cumulative_distance += p
        if debug: print 'loop 2, i', i, 'cumulative_distance', cumulative_distance
    return cumulative_distance
    
    
def total_distance(line_vertices):
    '''Return total distance along the line described by *line_vertices*.'''
    return sum((distance(*segment_vertices)
                for segment_vertices in line_segments(line_vertices)))
    
    
def point_on_line(dist, line_vertices):
    '''Return *x, y* of a point at distance *dist* along the profile described
    by *line_vertices*.'''
    cumulative_distance = 0
    for segment_vertices in line_segments(line_vertices):
        segment_distance = distance(*segment_vertices)
        if cumulative_distance + segment_distance > dist:
            return projected_point(segment_vertices[0], bearing_func(*segment_vertices), dist - cumulative_distance)
        elif cumulative_distance + segment_distance == dist:
            return segment_vertices[1]
        elif cumulative_distance + segment_distance < dist:
            cumulative_distance += segment_distance
            continue
    return projected_point(segment_vertices[1], bearing_func(*segment_vertices), dist - cumulative_distance)
    
    
def segment_distances_and_azimuths(line_vertices):
    '''Return *(distances, bearings)* for all the line segments in the profile
    described by *line_vertices*.'''
    distances = []
    bearings = []
    for segment_vertices in line_segments(line_vertices):
        distances.append(distance(*segment_vertices))
        bearings.append(angle_to_azimuth(bearing_func(*segment_vertices)))
    return distances, bearings


def pformat(line_vertices):
    '''Return summary multi-line string describing the profile.'''
    segdandaz = segment_distances_and_azimuths(line_vertices)
    return('Total length:\t\t{total:.0f} m\n'
           'Overall azimuth:\t{azimuth:.0f} deg\n'
           'Segments:\n{segments}'
           .format(total=total_distance(line_vertices),
                   azimuth=angle_to_azimuth(bearing_func(
                                line_vertices[0], line_vertices[-1])),
                   segments='\n'.join(['\t\t%.0f m at %.0f deg' % (d, az)
                                       for d, az in zip(*segdandaz)])
                  ))

def pprint(line_vertices):
    '''Print summary multi-line string describing the profile.'''
    print(pformat(line_vertices))