#!/usr/bin/env python
#
# Copyright (C) 2011  W. Trevor King
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

"""Print route distances and plot routes on a 2D map.

>>> from StringIO import StringIO
>>> stream = StringIO('\\n'.join([
...     '#lat  lng',
...     '41.293595 -73.552322',
...     '44.998213 -73.258667',
...     '45.021515 -74.752808',
...     '44.168809 -76.290894',
...     '43.463286 -76.433716',
...     '43.231597 -79.004517',
...     '42.846168 -78.894654',
...     '42.255357 -79.762573',
...     '41.994610 -79.773560',
...     '42.019101 -75.422974',
...     '41.478130 -74.994507',
...     '40.996795 -73.901138',
...     '']))
>>> latlngs = read_latlngs(stream)
>>> print('%.2f m' % route_distance(latlngs))
1679848.96 m
>>> plot_route(latlngs)
>>> _pylab.show()
"""

import numpy as _numpy
import pylab as _pylab
import pyproj as _pyproj


__version__ = '0.1'


def read_latlngs(stream):
    "Read in data from space-separated lat/lng lines with '#' comments."
    return _numpy.loadtxt(stream)

def route_distance(latlngs, geod=_pyproj.Geod(ellps='WGS84')):
    """Calculate the geodesic distance along an array of lat/lng pairs.

    If `geod` is not given, it defaults to the WGS 84 reference
    ellipsoide.
    """
    az12s,az21s,dists = geod.inv(latlngs[:-1,1], latlngs[:-1,0],
                                 latlngs[1:,1], latlngs[1:,0])
    return dists.sum()

def plot_route(latlngs, proj=None):
    """Plot an array of lat/lng pairs using the projection `proj`.

    If `proj` is not given, it defaults to a projection from the WGS
    84 reference ellipsoide (for lat/lng) to the closest UTM zone (for
    x/y).  The closes UTM zone is determined by averaging the input
    longitudes.

    Plots the line on the current pylab axis.
    """
    if proj is None:
        mlng = latlngs[:,1].mean()
        zone = int(round((mlng + 180) / 6.))
        proj = _pyproj.Proj(proj='utm', zone=zone, ellps='WGS84')
    xs,ys = proj(latlngs[:,1], latlngs[:,0])
    lines = _pylab.plot(xs, ys, '.-')

if __name__ == '__main__':
    from argparse import ArgumentParser

    p = ArgumentParser(
        description=__doc__.splitlines()[0],)
    p.add_argument(
        '--version', action='version', version='%(prog)s '+__version__)
    p.add_argument('coordfiles', metavar='COORDFILE', type=str, nargs='+',
                   help='Space-separated lat/lng line file')

    args = p.parse_args()

    _pylab.hold(True)

    for filename in args.coordfiles:
        latlngs = read_latlngs(filename)
        dist = route_distance(latlngs)
        print('%s distance: %g m (%g mi)' % (filename, dist, dist/1.6e3))
        plot_route(latlngs)

    _pylab.axis('equal')
    _pylab.show()