line.geojson

line_dist.py

# -*- coding: utf-8 -*-
import json
from itertools import izip

from vincenty import vincenty
from shapely.geometry import shape


def line_length(line):
    # our line is in EPSG:4326, shapely uses euacledian distance for .lengtg
    # so we have to use thg vincenty formula
    coords = line.coords
    return sum([vincenty((a[1], a[0]), (b[1], b[0]))
                for a, b in izip(coords, coords[1:])]) * 1000


def open_file(filename):
    with open(filename, 'r') as f:
        return f.read()


def read_geoms(filename):
    data = json.loads(open_file(filename))
    geoms = []
    for feature in data['features']:
        geoms.append(shape(feature['geometry']))
    return geoms


if __name__ == '__main__':

    line = read_geoms('line.geojson')[0]
    length = line_length(line)
    print 'line length: %s m' % length
    points = read_geoms('points.geojson')
    for point in points:
        print 'computing point ', point
        snapped = line.interpolate(line.project(point))
        fraction = line.project(point, normalized=True)
        print 'snapped to ', snapped
        from_start = length * fraction
        from_end = length - from_start
        print 'distance from start: %s m' % from_start
        print 'distance from end  : %s m' % from_end

points.geojson