#!/usr/bin/env perl
use strict;
use warnings;
use FindBin;
use lib ("$FindBin::RealBin/..", "$FindBin::RealBin/../lib");

use Encode qw(decode);
use I18N::Langinfo qw(langinfo CODESET);
use LWP::UserAgent;
use XML::LibXML;
use Getopt::Long;
use List::Util qw(min);
use POSIX qw(strftime);

my $codeset = langinfo(CODESET());
$codeset = lc $codeset; # 'UTF-8' is not recognized by emacs, but 'utf-8' is
binmode STDERR, ":encoding($codeset)";

# Command line args
my ($relation_id, $start_coord, $end_coord, $output, $backward, $cache, $in_xml, $title);
GetOptions(
    'id=i'     => \$relation_id,
    'start=s'  => \$start_coord,
    'end=s'    => \$end_coord,
    'out=s'    => \$output,
    'backward!'=> \$backward,
    'cache!'   => \$cache,
    'in-xml=s' => \$in_xml,
    'simplify=i' => \my $simplify,
    'title=s'  => \$title,
) or die "Usage: $0 --id RELATION_ID [--start lon,lat] [--end lon,lat] [--out file.gpx] [--backward] [--cache] [--in-xml /path/to/overpassresult.xml] [--simplify number]\n";
die "OSM relation ID required" unless $relation_id;

for ($title) {
    if (defined $_) {
	$_ = decode($codeset, $_);
    }
}

# API query (default: Overpass)
my $overpass_url = "https://overpass-api.de/api/interpreter";
my $query = <<EOF;
[out:xml][timeout:60];
relation($relation_id);
(._;>;);
out meta;
EOF

my $xml;
if (defined $in_xml) {
    $xml = do { open my $fh, $in_xml or die $!; local $/; <$fh> };
} else {
    my $ua = LWP::UserAgent->new(keep_alive => 1);

    my $res = $ua->post($overpass_url, Content => $query);
    die "Error fetching data: " . $res->status_line unless $res->is_success;
    $xml = $res->decoded_content(charset => 'none');
    if ($cache) {
	open my $ofh, ">", "/tmp/osmrel2gpx_response_${relation_id}.xml" or die $!;
	print $ofh $xml;
	close $ofh or die $!;
    }
}

# Parse XML
my $parser = XML::LibXML->new();
my $doc    = $parser->parse_string($xml);

# Extract nodes
my %nodes;
for my $node ($doc->findnodes('//node')) {
    my $id = $node->getAttribute('id');
    my $lat = $node->getAttribute('lat');
    my $lon = $node->getAttribute('lon');
    $nodes{$id} = [$lat, $lon];
}

# Extract ways and build ordered node list
my %ways;
for my $way ($doc->findnodes('//way')) {
    my $id = $way->getAttribute('id');
    my @nds = map { $_->getAttribute('ref') } $way->findnodes('./nd');
    $ways{$id} = \@nds;
}

# Build way order based on relation's member sequence (forward, backward)
my @members = $doc->findnodes('//relation/member[@type="way"]');

if ($backward) {
    @members = reverse @members;
}

my @ordered_nodes;
my $last_node;
my $last_way;

my $get_w_nodes = sub {
    my($mem_i, $last_node) = @_;

    my $mem = $members[$mem_i];
    my $role = $mem->getAttribute('role') || '';
    my $rid  = $mem->getAttribute('ref');
    my @w_nodes = @{ $ways{$rid} };

    # note: handling of role=backward used to be here, but often the
    # role is missing so in the end check manully if the way needs to
    # be reversed
    if (defined $last_node && $last_node == $w_nodes[-1]) {
	@w_nodes = reverse @w_nodes;
    }

    return ($rid, \@w_nodes);
};

WAY_LOOP:
for(my $mem_i=0; $mem_i <= $#members; $mem_i++) {
    my($rid, $w_nodes) = $get_w_nodes->($mem_i, $last_node);
    my @w_nodes = @$w_nodes;

    # First way: initialize list, but we need to decide the right direction first
    if (!defined $last_node) {
	my($next_rid, $next_w_nodes) = $get_w_nodes->($mem_i+1);
	if ($w_nodes[0] eq $next_w_nodes->[0] || $w_nodes->[0] eq $next_w_nodes->[-1]) {
	    # obviously need to reverse
	    @ordered_nodes = reverse @w_nodes;
	} else {
	    # probably starts in right direction (but may be wrong)
	    @ordered_nodes = @w_nodes;
	}
        $last_node = $ordered_nodes[-1];
        $last_way = $rid;
        next;
    }

    if (defined $last_node && $w_nodes[0] eq $w_nodes[-1]) {
    SEARCH_START_NODE: {
	    for(my $i=0; $i<=$#w_nodes; $i++) {
		if ($last_node eq $w_nodes[$i]) {
		    @w_nodes = (@w_nodes[$i..$#w_nodes], @w_nodes[1..$i-1]);
		    warn "Roundabout hack (1): reorder way $rid nodes after with new start index $i (nodes now: @w_nodes)\n";
		    if ($mem_i < $#members) {
			my($next_rid, $next_w_nodes) = $get_w_nodes->($mem_i+1);
		    SEARCH_END_NODE: {
			    for(my $j=1; $j<=$#w_nodes; $j++) {
				if ($next_w_nodes->[0] eq $w_nodes[$j] || $next_w_nodes->[-1] eq $w_nodes[$j]) {
				    @w_nodes = @w_nodes[0..$j];
				    warn "Roundabout hack (2): quit roundabout at (reordered) index $i (nodes now: @w_nodes)\n";
				    last SEARCH_END_NODE;
				}
			    }
			    warn "Roundabout hack (2): did not find roundabout quit node, expect problems (expected $next_w_nodes->[0] or $next_w_nodes-[-1], not found in @$w_nodes)\n";
			}
		    }
		    last SEARCH_START_NODE;
		}
	    }
	    warn "Roundabout hack (1): cannot find node in roundabout, expect problems\n";
	}
    }

    # Decide if we use this or the next member
    # (splitted route merged again)

    # If next way's start node == last_node, append (avoid duplicate node)
    if ($w_nodes[0] eq $last_node) {
	if ($mem_i < $#members) {
	    my($next_rid, $next_w_nodes) = $get_w_nodes->($mem_i+1, $last_node);
	    if ($next_w_nodes->[0] eq $last_node) {
		warn "Do not choose way $rid/node $w_nodes[0] (probably wrong direction), choose $next_rid/$next_w_nodes->[0]\n";
		@w_nodes = @$next_w_nodes;
		$rid = $next_rid;
		$mem_i++;
	    }
	}
        shift @w_nodes; # remove duplicate node
        push @ordered_nodes, @w_nodes;
        $last_node = $ordered_nodes[-1];
        $last_way = $rid;
        next;
    }

    # Otherwise skip this way (no connection)
    warn "Skipping way $rid (after last way $last_way); does not connect to last node $last_node\n";
    # Do not update last_node or last_way, continue with next member
}

# Extract relation name if available
my $relation_name = "OSM Relation $relation_id";
my ($rel_node) = $doc->findnodes('//relation[@id="'.$relation_id.'"]');
if ($rel_node) {
    my ($name_node) = $rel_node->findnodes('./tag[@k="name"]');
    $relation_name = $name_node->getAttribute('v') if $name_node;
}

# Optionally crop to start/end coordinates
sub haversine {
    my ($lat1,$lon1,$lat2,$lon2) = map { $_ * 3.14159265358979 / 180 } @_;
    my $R = 6371000; # Earth radius (m)
    my $dlat = $lat2-$lat1;
    my $dlon = $lon2-$lon1;
    my $a = sin($dlat/2)**2 + cos($lat1) * cos($lat2) * sin($dlon/2)**2;
    my $c = 2 * atan2(sqrt($a), sqrt(1-$a));
    return $R * $c;
}
sub find_nearest_node {
    my ($coord, $nodeids) = @_;
    my ($lon0, $lat0) = split /,/, $coord;
    my $min_dist = 1e38; my $best_id;
    for my $nid (@$nodeids) {
        my ($lat, $lon) = @{ $nodes{$nid} || [] };
        next unless defined $lat && defined $lon;
        my $dist = haversine($lat0, $lon0, $lat, $lon);
        if ($dist < $min_dist) {
            $min_dist = $dist; $best_id = $nid;
        }
    }
    return $best_id;
}

if ($start_coord) {
    my $start_id = find_nearest_node($start_coord, \@ordered_nodes);
    # Crop before start node
    my ($start_idx) = grep {$ordered_nodes[$_] eq $start_id} 0..$#ordered_nodes;
    @ordered_nodes = @ordered_nodes[$start_idx .. $#ordered_nodes] if defined $start_idx;
}
if ($end_coord) {
    my $end_id = find_nearest_node($end_coord, \@ordered_nodes);
    # Crop after end node
    my ($end_idx) = grep {$ordered_nodes[$_] eq $end_id} 0..$#ordered_nodes;
    @ordered_nodes = @ordered_nodes[0 .. $end_idx] if defined $end_idx;
}

my @lonlats = map { +{x => $nodes{$_}->[1], y => $nodes{$_}->[0]} } grep { $nodes{$_} } @ordered_nodes;
if (defined $simplify && @lonlats > $simplify) {
    require Route::LineSimplification;
    warn "Points before line simplification: " . scalar(@lonlats) . "\n";
    @lonlats = @{ Route::LineSimplification::visvalingam_whyatt(\@lonlats, ['elements', 1000]) };
    warn "Points after line simplification: " . scalar(@lonlats) . "\n";
}
@lonlats = map { [$_->{x}, $_->{y}] } @lonlats;

$title = $relation_name if !defined $title;

# Generate GPX
my $gpx = <<EOH;
<?xml version="1.0" encoding="UTF-8"?>
<gpx version="1.1" creator="rel2gpx.pl" xmlns="http://www.topografix.com/GPX/1/1">
<trk>
  <name>@{[ XML::LibXML::Text->new($title)->toString ]}</name>
  <trkseg>
EOH

for my $lonlat (@lonlats) {
    my ($lon, $lat) = @$lonlat;
    $gpx .= qq{    <trkpt lat="$lat" lon="$lon" />\n};
}

$gpx .= <<EOF;
  </trkseg>
</trk>
</gpx>
EOF

if ($output) {
    open my $fh, ">", $output or die "Cannot write $output: $!\n";
    binmode $fh, ':utf8';
    print $fh $gpx;
    close $fh;
} else {
    binmode STDOUT, ':utf8';
    print $gpx;
}

=head1 NAME

osmrel2gpx - create a GPX track file from a OSM relation

=head1 SYNOPSIS

    osmrel2gpx --id RELATION_ID [--start lon,lat] [--end lon,lat] [--out file.gpx] [--backward] [--cache] [--in-xml /path/to/overpassresult.xml] [--simplify number]

=head1 DESCRIPTION

Take a OpenStreetMap relation id, download the necessary data and
create a GPX track file out of it.

=head1 EXAMPLES

Generate a GPX file for Berlin-Usedom Radweg, starting in Berlin:

    osmrel2gpx --id 27727 --out berlin_usedom.gpx

Generate a GPX file for this relation, starting on Usedom:

    osmrel2gpx --id 27727 --backward --out usedom_berlin.gpx

Generate a shortened track using the given coordinates (here: start is
Bernau, goal is Prenzlau):

    osmrel2gpx --id 27727 --start 13.592188,52.676209 --end 13.861777,53.313756 --out bernau_prenzlau.gpx

Simplify the generated track to the given number of track points:

    osmrel2gpx --id 27727 --simplify 1024 --out berlin_usedom_simplified.gpx

=head1 BUGS AND LIMITATIONS

Theoretically it would be possible to determine the route direction
automatically if C<--start> and C<--end> are both specified. This is
not implemented, so it's necessary to manually specify C<--backward>
if needed.

Success is highly dependent on the quality of the relation. It is
likely that a generated track will be incomplete if
L<https://ra.osmsurround.org/analyzeRelation> detects problems for a
given relation.

=cut
