import math
import fiona
from shapely.geometry import shape, LineString, Point
from collections import OrderedDict

def dist(P3D1, P3D2):
    sqr = (P3D1[0] - P3D2[0])*(P3D1[0] - P3D2[0]) 
    sqr += (P3D1[1] - P3D2[1])*(P3D1[1] - P3D2[1])
    sqr += (P3D1[2] - P3D2[2])*(P3D1[2] - P3D2[2])
    return math.sqrt(sqr)

def pairs(seq):
    i = iter(seq)
    prev = item = i.next()
    for item in i:
        yield prev, item
        prev = item

def getnearest(point, line, distancefunc=None):
    """
    Get nearest index and point on a line for a given point external to the line
    
    point: some object representing a point in space 
    line: iterable composed of objects of the same type as point
    distancefunc: a function taking two points and returning a number
    
    returns: index, point on line that is nearest to the input point
    
    Not efficient for very long lines, because it searches along the line.
    """
    if not distancefunc:
        try:
            distancefunc = lambda x, y: (x - y) * (x - y)
        except:
            # if no distance function provided and objects can't be subtracted and multiplied, use constant
            distancefunc = lambda x, y: 1
    # This would be inefficient for long lines, but ok for simple applications
    # making sure we have an iterator; if line is an iterator, this doesn't make a difference
    xline = iter(line)
    # set nearest to first line element
    nearest = next(xline)
    nidx = 0
    # iterate through the rest of the line
    for idx, pt in enumerate(xline):
        if distancefunc(pt, point) < distancefunc(nearest, point):
            nearest = pt
            nidx = idx
    return nidx, nearest    

data_track = fiona.open('Iditarod_routes_AKAlbers.shp', 'r')
data_checkpoints = fiona.open('Iditarod_chkp_AKAlbers.shp', 'r')

print type(data_track)
print len(data_track)

for record in data_track:
    print record['geometry']['type']
    print record['geometry']['coordinates'][0], record['geometry']['coordinates'][-1]

fullinecoords = []
for line in data_track:
    fullinecoords.extend(line['geometry']['coordinates'])

fulline = {
    'geometry': {
        'coordinates': fullinecoords,
        'type': 'LineString'    
    },
    'id': 0,
    'properties': OrderedDict([(u'Name', u'Iditarod_northern'), (u'descriptio', None), (u'timestamp', None), (u'begin', None), (u'end', None), (u'altitudeMo', None), (u'tessellate', 1), (u'extrude', -1), (u'visibility', -1)])    
}
fulltrack = shape(fulline['geometry'])
fulltrack.geom_type

chkp_ontrack = {}
for checkpoint in data_checkpoints:
    name = checkpoint['properties']['Name']
    coords = checkpoint['geometry']['coordinates']
    nearestidx, nearestpoint = getnearest(coords, fulltrack.coords, dist)
    chkp_ontrack[name] = {}
    chkp_ontrack[name]['idx'] = nearestidx
    chkp_ontrack[name]['coords'] = nearestpoint
chkp_ontrack = OrderedDict(sorted(chkp_ontrack.items(), key=lambda item: item[1]['idx']))
print chkp_ontrack

totaldist = 0
for chkpt, nextchkpt in pairs(chkp_ontrack):
    idx, nextidx = chkp_ontrack[chkpt]['idx'], chkp_ontrack[nextchkpt]['idx']
    partialtrack = LineString(fulltrack.coords[idx:nextidx+1])
    totaldist = totaldist + partialtrack.length
    print '%s --> %s' % (chkpt, nextchkpt)
    print '  Distance in km: %s' % format(partialtrack.length/1000, '.1f')
    print '  Distance in miles: %s' % format(partialtrack.length/1600, '.1f')
print
print 'Total distance summed up: %s km -- %s miles' % (format(totaldist/1000, '.1f'), format(totaldist/1600, '.1f'))
print 'Total distance from full track: %s km -- %s miles' % (format(fulltrack.length/1000, '.1f'), format(fulltrack.length/1600, '.1f'))

len(fulltrack.coords)

waypoints = {}
import itertools
for checkpoint in data_checkpoints:
    name = checkpoint['properties']['Name']
    chkptshape = Point(checkpoint['geometry']['coordinates'])
    a = fulltrack.project(chkptshape)
    waypoints[name] = a
waypoints = OrderedDict(sorted(waypoints.items(), key=lambda x: x[1]))
print waypoints

prev = 0
for pt1, pt2 in pairs(waypoints):
    print '%s --> %s' % (pt1, pt2)
    print '  Distance in km: %s' % format((waypoints[pt2] - waypoints[pt1])/1000, '.1f')
    print '  Distance in miles: %s' % format((waypoints[pt2] - waypoints[pt1])/1600, '.1f')
print
print 'Total distance: %s km -- %s miles' % (format(waypoints[pt2]/1000, '.1f'), format(waypoints[pt2]/1600, '.1f'))