#!/usr/bin/perl

=head1 NAME

osmtrace - add a gps trace to an Open Street Map map

=head1 SYNOPSIS

osmtrace [ --lines ] [ --speed-max kmh ] [--mark_time seconds] [--mark_segments] [--tiles n] tracefile

=head1 DESCRIPTION

osmtrace reads a gps trace file in tangoGPS, GPSJinny or GPX format, downloads map tiles
from the OSM server and produces a PNG image with the trace drawn on the
map.

Options

=over

=item --lines

Connect the points with lines. 

=item --speed-max kmh

The position is color-coded with the speed, ranging from black (0) through
green (0.2 * max), yellow (0.4 * max), red (0.6 * max), purple (0.8 *
max) to blue (max).  This option sets the maximum speed in km/h. Default
is 200.

=item --mark_time seconds

Put a time mark every I<seconds> seconds.

=item --mark_segments

Put a time mark at the start end end of each segment.

=item --tiles n

Scale to use at least n tiles (1 tile = 256 pixels) in the longer
direction.

=cut

use warnings;
use strict;

use Math::Trig;
use POSIX qw(floor ceil);
use LWP::UserAgent;
use Imager;
use Geo::Distance;
use HTTP::Date;
use List::Util qw(min max);
use Pod::Usage;
use Getopt::Long;
use POSIX qw(strftime);
use XML::Simple;

my $tiles = 3;

my ($from, $to);

sub readtrace {
    my ($filename) = @_;
    open (my $trace_fh, '<', $filename) or die "cannot open $filename: $!";
    if ($filename =~ m/\.gpx$/i) {
	return read_gpx_trace($trace_fh);
    } else {
	my $magic;
	if (read($trace_fh, $magic, 4) == 4 && $magic eq "GJT\0") {
	    seek($trace_fh, 0, 0);
	    return read_gpsjinny_trace($trace_fh);
	} else {
	    seek($trace_fh, 0, 0);
	    return read_tangogps_trace($trace_fh);
	}
    }
    return;
}

sub read_tangogps_trace {
    my ($trace_fh) = @_;
    my ($min_lat, $max_lat, $min_long, $max_long);
    my @points;
    while (<$trace_fh>) {
	chomp;
	my ($lat, $long, $elev, undef, $head, undef, $time) = split(/,/);
	push @points, {
			latitude  => $lat,
			longitude => $long,
			elevation => $elev,
			timestr   => $time,
		      };
	$min_lat = $lat if (!defined $min_lat || $lat < $min_lat);
	$max_lat = $lat if (!defined $max_lat || $lat > $max_lat);
	$min_long = $long if (!defined $min_long || $long < $min_long);
	$max_long = $long if (!defined $max_long || $long > $max_long);
    }
    return {
		min_latitude  => $min_lat,
		max_latitude  => $max_lat,
		min_longitude => $min_long,
		max_longitude => $max_long,
		points        => \@points,
	   };
}

sub read_gpsjinny_trace {
    my ($fh) = @_;
    my $buf;

    read($fh, $buf, 4);
    die unless $buf eq "GJT\0";

    read($fh, $buf, 4);
    my $version = unpack("N", $buf);
    die unless $version == 1;

    my ($min_lat, $max_lat, $min_long, $max_long);
    my @points;
    while (read($fh, $buf, 16) == 16) {
	my ($ts, $x, $y, $z) = unpack("NNNN", $buf);
	my $long = $x / 1E7 - 90;
	my $lat  = $y / 1E7 - 180;
	my $elev = $z / 1E2;
	#print STDERR "$ts, $long $lat $elev\n";
	push @points, {
			latitude  => $lat,
			longitude => $long,
			elevation => $elev,
			timenum   => $ts,
		      };
	$min_lat = $lat if (!defined $min_lat || $lat < $min_lat);
	$max_lat = $lat if (!defined $max_lat || $lat > $max_lat);
	$min_long = $long if (!defined $min_long || $long < $min_long);
	$max_long = $long if (!defined $max_long || $long > $max_long);
    }
    return {
		min_latitude  => $min_lat,
		max_latitude  => $max_lat,
		min_longitude => $min_long,
		max_longitude => $max_long,
		points        => \@points,
	   };
}

my $geo;
sub read_gpx_trace {
    my ($fh) = @_;
    my $xml = XMLin($fh, ForceArray => [ 'trkseg', 'trkpt' ]);
    my @points;
    my ($min_lat, $max_lat, $min_long, $max_long);
    my $last_point;
    for my $trkseg (@{ $xml->{trk}{trkseg} }) {
        my $first_in_segment = 1;
	for (@{ $trkseg->{trkpt} }) {
            my $point = {
			    latitude  => $_->{lat},
			    longitude => $_->{lon},
			    elevation => $_->{ele},
			    timestr   => $_->{time},
			};
            if ($from || $to) {
                $point->{timenum} = timenum($point);
                next if defined $from && $point->{timenum} < $from;
                next if defined $to   && $point->{timenum} > $to;
            }
            if ($first_in_segment) {
                $point->{segment_start} = 1;
                $first_in_segment = 0;
            }
	    push @points, $point;
	    $min_lat = $_->{lat} if (!defined $min_lat || $_->{lat} < $min_lat);
	    $max_lat = $_->{lat} if (!defined $max_lat || $_->{lat} > $max_lat);
	    $min_long = $_->{lon} if (!defined $min_long || $_->{lon} < $min_long);
	    $max_long = $_->{lon} if (!defined $max_long || $_->{lon} > $max_long);

            if (defined $last_point) {
		$geo = Geo::Distance->new() unless $geo;
                my $dist = $geo->distance('meter',
                                          $last_point->{longitude}, $last_point->{latitude},
                                          $point->{longitude}, $point->{latitude},
                                         );
                $point->{cum_dist} = $last_point->{cum_dist} + $dist;

            } else {
                $point->{cum_dist} = 0;
            }
            $last_point = $point;
	}
        $last_point->{segment_end} = 1;
    }
    return {
		min_latitude  => $min_lat,
		max_latitude  => $max_lat,
		min_longitude => $min_long,
		max_longitude => $max_long,
		points        => \@points,
	   };
}

sub x_from_long_deg {
    my ($long_deg, $zoom) = @_;
    return $long_deg / 360 * (1 << $zoom) + (1 << $zoom) / 2;
}

sub y_from_lat_deg {
    my ($lat_deg, $zoom) = @_;
    my $lat_rad = $lat_deg * pi / 180;
    my $lat_xf = log(abs(tan(($lat_rad+pi/2)/2)));
    return (1 - $lat_xf / pi) / 2 * (1 << $zoom);
}


sub get_map {
    my ($tr) =  @_;
    my $zoom = 0;
    my $min_x;
    my $max_x;
    my $max_y;
    my $min_y;

    my $tileserver = "http://b.tile.openstreetmap.org/";
    for (;;) {
	$min_x = floor(x_from_long_deg($tr->{min_longitude}, $zoom));
	$max_x = floor(x_from_long_deg($tr->{max_longitude}, $zoom));
	$max_y = floor(y_from_lat_deg($tr->{min_latitude}, $zoom));	# sic!
	$min_y = floor(y_from_lat_deg($tr->{max_latitude}, $zoom));	# sic!
	last if ($max_y - $min_y >= $tiles ||
		 $max_x - $min_x >= $tiles ||
		 $zoom >= 17);
	   
	$zoom++;
    }

    my $ua = LWP::UserAgent->new;
    $ua->env_proxy;

    my $img = Imager->new(xsize    => ($max_x - $min_x + 1) * 256,
			  ysize    => ($max_y - $min_y + 1) * 256,
			  channels => 3);

    for my $y (floor($min_y) .. floor($max_y)) {
	for my $x (floor($min_x) .. floor($max_x)) {
	    my $url =
	    sprintf("%s/%d/%d/%d.png",
		    $tileserver, $zoom, $x, $y);
	    print "url      : $url\n";
	    my $response = $ua->get($url);
	    my $tile = Imager->new();
	    $tile->read(data => $response->content, type => 'png');
	    $img->paste(left => ($x - $min_x) * 256,
			top  => ($y - $min_y) * 256,
			img  => $tile);
	}
    }
    my @map = map { int($_/2) } 0..255;
    $img->map(all => \@map);
    return ($img, $zoom, $min_x, $min_y);
}

sub speed {
    my ($tr, $i) = @_;
    my $i0 = max(0, $i - 5);
    my $i1 = min($#{ $tr->{points} }, $i + 5);

    my $p0 = $tr->{points}[$i0];
    my $p1 = $tr->{points}[$i1];
    $geo = Geo::Distance->new() unless $geo;
    my $dist = $geo->distance('meter',
			      $p0->{longitude}, $p0->{latitude},
			      $p1->{longitude}, $p1->{latitude},
			     );
    my $t0 = $p0->{timenum} // str2time($p0->{timestr});
    my $t1 = $p1->{timenum} // str2time($p1->{timestr});
    if ($t1 > $t0) {
	my $speed = $dist / ($t1 - $t0) * 3.6;
	#print STDERR "$speed\n";
	return $speed;
    } else {
	print STDERR "non-positive time difference: [$i0 $t0: $p0->{longitude}, $p0->{latitude}] -> [$i1 $t1: $p1->{longitude}, $p1->{latitude}]\n";
	return 0;
    }
}


my $speed_max = 200;
sub color_from_speed {
    my ($speed) = @_;
    if ($speed < 1/5 * $speed_max) {
	my $ls = 0/5 * $speed_max;
	my $hs = 1/5 * $speed_max;
	my $rs = ($speed - $ls) / ($hs - $ls);
	return Imager::Color->new(0, $rs * 255, 0);
    } elsif ($speed < 2/5 * $speed_max) {
	my $ls = 1/5 * $speed_max;
	my $hs = 2/5 * $speed_max;
	my $rs = ($speed - $ls) / ($hs - $ls);
	return Imager::Color->new($rs * 255, 255, 0);
    } elsif ($speed < 3/5 * $speed_max) {
	my $ls = 2/5 * $speed_max;
	my $hs = 3/5 * $speed_max;
	my $rs = ($speed - $ls) / ($hs - $ls);
	return Imager::Color->new(255, 255 * (1-$rs), 0);
    } elsif ($speed < 4/5 * $speed_max) {
	my $ls = 3/5 * $speed_max;
	my $hs = 4/5 * $speed_max;
	my $rs = ($speed - $ls) / ($hs - $ls);
	return Imager::Color->new(255, 0, 255 * $rs);
    } elsif ($speed < 5/5 * $speed_max) {
	my $ls = 3/5 * $speed_max;
	my $hs = 4/5 * $speed_max;
	my $rs = ($speed - $ls) / ($hs - $ls);
	return Imager::Color->new(255 * (1-$rs), 0, 255);
    } else {
	return Imager::Color->new(0, 0, 255);
    }
}

sub timestr {
    my ($p) = @_;
    return $p->{timestr} // strftime("%Y-%m-%dT%H:%M:%S", localtime($p->{timenum}));
}

sub timenum {
    my ($p) = @_;
    return $p->{timenum} // str2time($p->{timestr});
}


sub add_trace {
    my ($img, $zoom, $min_x, $min_y, $tr, $lines,
        $mark_time, $mark_distance, $mark_segments,
        $circle_radius) = @_;
    my $white = Imager::Color->new(255, 255, 255);
    my $black = Imager::Color->new(0, 0, 0);
    my ($last_x, $last_y);
    my $nextmark = timenum($tr->{points}[0]);
    #my $font = Imager::Font->new(file => '/usr/share/fonts/truetype/freefont/FreeMono.ttf');
    my $font = Imager::Font->new(file => '/usr/share/texmf/fonts/opentype/public/lm/lmsansdemicond10-regular.otf');
    for my $i (0 .. $#{$tr->{points}}) {
	my $p = $tr->{points}[$i];

	my $x = x_from_long_deg($p->{longitude}, $zoom);
	my $y = y_from_lat_deg($p->{latitude}, $zoom);

	my $color = color_from_speed(speed($tr, $i));
	if ($lines) {
	    if (defined $last_x) {
		$img->line(x2    => ($x - $min_x) * 256, 
			   y2    => ($y - $min_y) * 256,
			   x1    => ($last_x - $min_x) * 256,
			   y1    => ($last_y - $min_y) * 256,
			   color => $color);
	    }
	    $last_x = $x;
	    $last_y = $y;
        } elsif ($circle_radius) {
            $img->circle(x     => ($x - $min_x) * 256, 
			 y     => ($y - $min_y) * 256,
                         r     => $circle_radius,
			 color => $color);

	} else {
	    $img->setpixel(x     => ($x - $min_x) * 256, 
			   y     => ($y - $min_y) * 256,
			   color => $color);
	}
        my $txt;
	if ($mark_time && timenum($p) >= $nextmark) {
	    $txt = timestr($p);
	    if ($mark_distance) {
		$txt .= sprintf(" %.0f", $p->{cum_dist});
	    }
	    $nextmark += $mark_time;
        }
        if ($mark_segments && ($p->{segment_start} || $p->{segment_end})) {
            $txt = timestr($p) . sprintf(" %.0f", $p->{cum_dist});
        }
        if ($txt) {
	    my $x0 = ($x - $min_x) * 256;
	    my $y0 = ($y - $min_y) * 256;
	    for my $dy (-1 .. 1) {
		for my $dx (-1 .. 1) {

		    $img->string(x      => $x0 + $dx,
				 y      => $y0 + $dy,
				 color  => $white,
				 font   => $font,
				 string => $txt,
				 size   => 15,
				 aa     => 1,
				);
		}
	    }
	    $img->string(x      => $x0,
			 y      => $y0,
			 color  => $black,
			 font   => $font,
			 string => $txt,
			 size   => 15,
			 aa     => 1,
			);
	}
    }
}

my $lines = 0;
my $mark_time;
my $mark_distance;
my $mark_segments;
my $circle_radius;
GetOptions(
    'lines'           => \$lines,
    'speed-max=f'     => \$speed_max,
    'mark_time=f'     => \$mark_time,
    'mark_distance'   => \$mark_distance,
    'mark_segments'   => \$mark_segments,
    'tiles=i'         => \$tiles,
    'circle_radius=f' => \$circle_radius,
    'from=s'          => \$from,
    'to=s'            => \$to,
) or pod2usage(verbose => 1);

pod2usage(verbose => 1) unless (@ARGV == 1);

for ($from, $to) {
    if ($_) {
        $_ = str2time($_);
    }
}

my $tr = readtrace($ARGV[0]);
print "latitude : $tr->{min_latitude} .. $tr->{max_latitude}\n";
print "longitude: $tr->{min_longitude} .. $tr->{max_longitude}\n";
print "time: ",
      timestr($tr->{points}[0]),
      " .. ",
      timestr($tr->{points}[-1]),
      "\n";
print "distance: ", $tr->{points}[-1]{cum_dist}, " meters\n";

my ($img, $zoom, $min_x, $min_y) = get_map($tr);

add_trace($img, $zoom, $min_x, $min_y, $tr, $lines,
    $mark_time, $mark_distance, $mark_segments,
    $circle_radius);

my $imgfilename = $ARGV[0];
$imgfilename =~ s/(\.[^.]+)?$/.png/;
$img->write(file => $imgfilename, type => 'png');

exit(0);
# sw=4 ts=8 expandtab