#!/usr/bin/perl

## see POD doc at bottom of this file. 
## html docs, this sourcecode, and changelogs now automatically published by git-commit
## released under GPL v2, use at your own risk
## patches for corrections/improvements welcome
	
use strict;
use warnings;
use Text::ParseWords; # to replace CSV_XS
use IO::Handle;
use Time::Local;
use Getopt::Long;
use Pod::Usage;

my ($incsvfile, $outcsvfile, $outgpxfile1, $outgpxfile2, $outkmlfile1, $outkmlfile2);

my ($INFILE, $OUTFILE);

# ignore all lines with HDOP bigger than this
my $MAXHDOP = 2.5;

# approx radius of Earth in meters.  True radius varies from
# 6357km (polar) to 6378km (equatorial).
use constant EARTH_RADIUS => 6367000;
use constant PI => 4 * atan2(1, 1);
use constant EARTH_DIAMETER => EARTH_RADIUS * 2 * PI;
use constant METRE => 360 / EARTH_DIAMETER; # in degrees

my $MIN_DIST = 20; 
use constant MTIME => 9;
 
# field numbers for particular headings we need in the csv file
# turn this into a single structured variable?
my ($INDEX, $HDOP, $RCR, $DATE, $TIME, $LATITUDE, $LONGITUDE, $NS, $EW);
my (@csv_timestamp, @csv_lat, @csv_long, @csv_seconds);
# button data, same comment as above
my (@buttonpress, $buttoncount, $refpress, $reftime, $startpress, $stoppress);

my $inputlines;
my ($minlat, $minlong, $maxlat, $maxlong);

my $man = 0;
my $help = 0;
my $REVERSE;
my $DRYRUN;

GetOptions (	'help|?'	=> 	\$help,
		'man' 		=> 	\$man,
		'outgpxlong=s' 	=> 	\$outgpxfile1,
		'outgpxshort=s' => 	\$outgpxfile2,
		'outkmllong=s'	=>	\$outkmlfile1,
		'outkmlshort=s'	=>	\$outkmlfile2,
		'hdop=f'	=>	\$MAXHDOP,
		'mindist=f'	=>	\$MIN_DIST,
		'reverse!'	=>	\$REVERSE,
		'dryrun!'	=>	\$DRYRUN) or pod2usage(2);

pod2usage(-exitstatus => 0, -verbose => 2) if $man;

$incsvfile =	shift;

pod2usage(1) if $help or not defined $incsvfile;

print "loading csvfile $incsvfile\n";
open $INFILE, "<", $incsvfile or die "$incsvfile: $!";

$MIN_DIST = $MIN_DIST * METRE;

# read header
my $row = <$INFILE>;
my @fields = parse_csv($row);

# find position of the fields we are interested in
# (could we make this simpler?)
foreach my $col (0 .. --$#fields) {
  if ($fields[$col] eq "INDEX" )     {$INDEX = $col;}
  if ($fields[$col] eq "HDOP" )      {$HDOP = $col;}
  if ($fields[$col] eq "RCR" )       {$RCR = $col;}
  if ($fields[$col] eq "DATE" )      {$DATE = $col;}
  if ($fields[$col] eq "TIME" )      {$TIME = $col;}
  if ($fields[$col] eq "LATITUDE" )  {$LATITUDE = $col;}
  if ($fields[$col] eq "LONGITUDE" ) {$LONGITUDE = $col;}
  if ($fields[$col] eq "N/S" )       {$NS = $col;}
  if ($fields[$col] eq "E/W" )       {$EW = $col;}
}

# 8 million ways to die ...
die "Error - could not find INDEX at top of $incsvfile\n"     unless defined $INDEX;
die "Error - could not find HDOP at top of $incsvfile\n"      unless defined $HDOP;
die "Error - could not find RCR at top of $incsvfile\n"       unless defined $RCR;
die "Error - could not find DATE at top of $incsvfile\n"      unless defined $DATE;
die "Error - could not find TIME at top of $incsvfile\n"      unless defined $TIME;
die "Error - could not find LATITUDE at top of $incsvfile\n"  unless defined $LATITUDE;
die "Error - could not find LONGITUDE at top of $incsvfile\n" unless defined $LONGITUDE;
die "Error - could not find N/S at top of $incsvfile\n"       unless defined $NS;
die "Error - could not find E/W at top of $incsvfile\n"       unless defined $EW;

# now parse all the lines in the file, dropping those
# where HDOP is too high (too much positional error)
# the QT1000 gives *unsigned* csv_lat/csv_long out, so convert to signed.

my $count = 0;
while (my $row = <$INFILE>) {
	$inputlines++;
	my @fields = parse_csv($row);
	if ($fields[$HDOP] < $MAXHDOP) {
		$count++;			

		$fields[$LATITUDE] = -$fields[$LATITUDE] if ($fields[$NS] eq "S");
		$csv_lat[$count] = $fields[$LATITUDE];
		 
		$fields[$LONGITUDE] = -$fields[$LONGITUDE] if ($fields[$EW] eq "W");
		$csv_long[$count] = $fields[$LONGITUDE];
		
		$csv_timestamp[$count] = "$fields[$DATE]T$fields[$TIME]Z";
		$csv_timestamp[$count] =~ s/\//-/g;

		$csv_seconds[$count] = &to_seconds($csv_timestamp[$count]);
		
		if ($fields[$RCR] eq "B") {
			$buttoncount++;
			$buttonpress[$buttoncount]=$count;
		}
	}
}
close $INFILE;

if (defined $buttoncount) 	{
	$refpress = 1;
	if ($buttoncount>2) {
		$stoppress = $buttoncount;
		$startpress = 2;
	}
}


print "Ref button press was at index # $buttonpress[$refpress], timestamp was $csv_timestamp[$refpress] GMT\n" if defined $refpress;
print "Start button press was at index # $buttonpress[$startpress]\n" if defined $startpress;
print "Stop button press was at index # $buttonpress[$stoppress]\n" if defined $stoppress;
print "\n";

if (defined $outgpxfile1) {
	&out_gpx_file($outgpxfile1,1,$count,0);
}

if (defined $outgpxfile2) {
	if (defined $stoppress) {
		&out_gpx_file($outgpxfile2,$buttonpress[$startpress],$buttonpress[$stoppress],$MIN_DIST)
	}
	else {
		&out_gpx_file($outgpxfile2,1,$count,$MIN_DIST)
	}
}

if (defined $outkmlfile1) {
	&out_kml_file($outkmlfile1,1,$count,0);
}

if (defined $outkmlfile2) {
	if (defined $stoppress) {
		&out_kml_file($outkmlfile2,$buttonpress[$startpress],$buttonpress[$stoppress],$MIN_DIST)
	}
	else {
		&out_kml_file($outkmlfile2,1,$count,$MIN_DIST)
	}
}

exit;

##################### subroutines below this line #########################
sub to_seconds{
	my $date = shift;
	my $time;
	$date =~ s/ /:/g;
	$date =~ s/T/:/g;
	$date =~ s/-/:/g;
	$date =~ s/Z//g;

	my @mons = split(/\:/, $date);
	my $year = $mons[0];
	my $mon  = $mons[1];
	my $day  = $mons[2];
	my $hr   = $mons[3];
	my $min  = $mons[4];
	my $sec  = $mons[5];
	unless (($mon < 1) || ($mon > 12)) { ##Bad Dates Chk
		$time = timelocal($sec, $min, $hr, $day, $mon-1, $year-1900); #converts to Epoch
	}
	chomp($time);
	return $time;
}

sub dist_squared{
# distance in degrees^2 between two points,
# it assumes points are close together, so don't worry about being exact
# treat world as flat sheet 	
	my $n1 = shift;
	my $n2 = shift;
	
	my $y = $csv_lat[$n1] - $csv_lat[$n2];
	my $x =	$csv_long[$n1] - $csv_long[$n2];
	
	my $hypSquared = $x*$x + $y*$y;
	
	return $hypSquared;
}

sub out_gpx_file{
	my $filename = shift;
	my $start    = shift;
	my $stop     = shift;
	my $mindist  = shift;
	
	my $mindistSquared = $mindist * $mindist; ##to avoid lots of square roots later
	
	if (!defined $REVERSE) {
		print "attempting to output to $filename, from $start to $stop\n";
	}
	else {
		print "attempting to output to $filename, in reverse from $stop to $start\n";
	}

# find bounding box for 'track'.
### do we need to limit waypoints to be inside this box??	

	$minlat = 999;
	$minlong = 999;
	$maxlat =-999;
	$maxlong =-999;

	for (my $count=$start;$count<=$stop;$count++) {
		$maxlat = $csv_lat[$count] if ($csv_lat[$count] > $maxlat);
		$minlat = $csv_lat[$count] if ($csv_lat[$count] < $minlat);
		$maxlong = $csv_long[$count] if ($csv_long[$count] > $maxlong);
		$minlong = $csv_long[$count] if ($csv_long[$count] < $minlong);
	}
	
	$filename ="/dev/null" if defined $DRYRUN; 

	open $OUTFILE, ">", $filename or die "$filename: $!";
	&print_gpx_head;	
	
	if (defined $stoppress) {
		my $waypointcount;
		if (!defined $REVERSE) {
			for (my $count=$startpress;$count<=$stoppress;$count++) {
					&print_gpx_wpt($count);
					$waypointcount++;
			}
			print "$waypointcount waypoints dumped to gpx file\n";
		}
		else {
			for (my $count=$stoppress;$count>=$startpress;$count--) {
					&print_gpx_wpt($count);
					$waypointcount++;
			}
			print "$waypointcount waypoints dumped to gpx file in reverse\n";
		}
	}


	&print_gpx_trk_head;
	my $lastoutput;
	my $outputlines = 0;
	if (!defined $REVERSE) {
		for (my $count=$start;$count<=$stop;$count++) {
			if (($count == $start) or ($count == $stop) or (&dist_squared($lastoutput,$count) >= $mindistSquared)) {
				&print_gpx_trk($count);
				$lastoutput=$count;
				$outputlines++;
			}	
		}
	}
	else {
		for (my $count=$stop;$count>=$start;$count--) {
			if (($count == $start) or ($count == $stop) or (&dist_squared($lastoutput,$count) >= $mindistSquared)) {
				&print_gpx_trk($count);
				$lastoutput=$count;
				$outputlines++;
			}	
		}
	}
	&print_gpx_trk_tail;
	
	&print_gpx_tail;
	close $OUTFILE;
	my $diff = $csv_seconds[$stop]-$csv_seconds[$start];
	print "$inputlines read in, $outputlines written out\n";
	print "length of trace is $diff seconds\n\n";
}

sub print_gpx_head{
	print $OUTFILE "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n";
	print $OUTFILE "<gpx\n";
	print $OUTFILE " version=\"1.0\"\n";
	print $OUTFILE "creator=\"csv2kml - http://www.morgad.co.uk\"\n";
	print $OUTFILE "xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"\n";
	print $OUTFILE "xmlns=\"http://www.topografix.com/GPX/1/0\"\n";
	print $OUTFILE "xsi:schemaLocation=\"http://www.topografix.com/GPX/1/0 http://www.topografix.com/GPX/1/0/gpx.xsd\">\n";
	print $OUTFILE "<bounds minlat=\"$minlat\" minlon=\"$minlong\" maxlat=\"$maxlat\" maxlon=\"$maxlong\"/>\n";
}

sub print_gpx_wpt{
	my $index = shift;
	print $OUTFILE "<wpt lat=\"$csv_lat[$buttonpress[$index]]\" lon=\"$csv_long[$buttonpress[$index]]\">\n";
	print $OUTFILE "<name>waypoint $index</name>\n";
	print $OUTFILE "</wpt>\n";
}


sub print_gpx_trk_head{
	print $OUTFILE "<trk>\n";
	print $OUTFILE "<trkseg>\n";
}

sub print_gpx_trk{
	my $index = shift;
	print $OUTFILE "<trkpt lat=\"$csv_lat[$index]\" lon=\"$csv_long[$index]\">\n";
	print $OUTFILE "<time>$csv_timestamp[$index]</time>\n";
	print $OUTFILE "</trkpt>\n";
}

sub print_gpx_trk_tail{
	print $OUTFILE "</trkseg>\n";
	print $OUTFILE "</trk>\n";
}

sub print_gpx_tail{
	print $OUTFILE "</gpx>\n";
}

sub out_kml_file{
### //FIXME// - is it worth turning outgpx and outkml into one subroutine?
	my $filename = shift;
	my $start    = shift;
	my $stop     = shift;
	my $mindist  = shift;
	
	my $mindistSquared = $mindist * $mindist; ##to avoid lots of square roots later
	
	print "attempting to output to $filename, from $start to $stop\n";
	
	$minlat = 999;
	$minlong = 999;
	$maxlat =-999;
	$maxlong =-999;
	
	$filename ="/dev/null" if defined $DRYRUN; 

	open $OUTFILE, ">", $filename or die "$filename: $!";
	&print_kml_head;	
	
	if (defined $stoppress) {
		my $waypointcount;
		if (!defined $REVERSE) {
			for (my $count=$startpress;$count<=$stoppress;$count++) {
					&print_kml_wpt($count);
					$waypointcount++;
			}
			print "$waypointcount waypoints dumped to kml file\n";
		}
		else {
			for (my $count=$stoppress;$count>=$startpress;$count--) {
					&print_kml_wpt($count);
					$waypointcount++;
			}
			print "$waypointcount waypoints dumped to kml file in reverse\n";
		}
	}

	&print_kml_trk_head;
	my $lastoutput;
	my $outputlines = 0;
	if (!defined $REVERSE) {
		for (my $count=$start;$count<=$stop;$count++) {
			if (($count == $start) or ($count == $stop) or (&dist_squared($lastoutput,$count) >= $mindistSquared)) {
				&print_kml_trk($count);
				$lastoutput=$count;
				$outputlines++;
			}	
		}
	}
	else {
		for (my $count=$stop;$count>=$start;$count--) {
			if (($count == $start) or ($count == $stop) or (&dist_squared($lastoutput,$count) >= $mindistSquared)) {
				&print_kml_trk($count);
				$lastoutput=$count;
				$outputlines++;
			}	
		}
	}
	&print_kml_trk_tail;
	
	&print_kml_tail;
	close $OUTFILE;
	my $diff = $csv_seconds[$stop]-$csv_seconds[$start];
	print "$inputlines read in, $outputlines written out\n";
	print "length of trace is $diff seconds\n\n";
}

sub print_kml_head{
	print $OUTFILE "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n";
	print $OUTFILE "<kml xmlns=\"http://earth.google.com/kml/2.0\">\n";
	print $OUTFILE "<Document>\n";
	print $OUTFILE "  <name>My Trips</name>\n";
	print $OUTFILE "  <open>1</open>\n";
	print $OUTFILE "  <Folder>\n";
	print $OUTFILE "    <name>My Places</name>\n";
	print $OUTFILE "    <open>0</open>\n";
}


sub print_kml_wpt{
	my $index = shift;
#	print $OUTFILE "<wpt lat=\"$csv_lat[$buttonpress[$index]]\" lon=\"$csv_long[$buttonpress[$index]]\">\n";
#	print $OUTFILE "<name>waypoint $index</name>\n";
#	print $OUTFILE "</wpt>\n";
}


sub print_kml_trk_head{
#close waypoint folder first
	print $OUTFILE "  </Folder>\n";
	print $OUTFILE "  <Folder>\n";
	print $OUTFILE "    <name>My Path</name>\n";
	print $OUTFILE "    <open>0</open>\n";
	print $OUTFILE "    <Placemark>\n";
	print $OUTFILE "      <name>Path-$incsvfile</name>\n";
	print $OUTFILE "      <Style>\n";
	print $OUTFILE "        <LineStyle>\n";
	print $OUTFILE "          <color>ff00ffff</color>\n";
	print $OUTFILE "          <width>3.0</width>\n";
	print $OUTFILE "        </LineStyle>\n";
	print $OUTFILE "      </Style>\n";
	print $OUTFILE "      <LineString>\n";
	print $OUTFILE "        <tessellate>1</tessellate>\n";
	print $OUTFILE "        <coordinates>\n";

}

sub print_kml_trk{
	my $index = shift;
	print $OUTFILE "          $csv_long[$index],$csv_lat[$index]\n";
}

sub print_kml_trk_tail{
	print $OUTFILE "        </coordinates>\n";
	print $OUTFILE "      </LineString>\n";
	print $OUTFILE "    </Placemark>\n";
	print $OUTFILE "  </Folder>\n";
}

sub print_kml_tail{
	print $OUTFILE "  </Document>\n";
	print $OUTFILE "</kml>\n";
}


# based on code borrowed from http://www.indo.com/distance/dist.pl
# convert a string which looks like "34:45:12N" into
# degrees.  Also accepts "34.233N" etc.
sub parse_location {
    my($str) = @_;
    my($latlong);

	if (defined $str) {
# turn 34 deg 45' 12.34" into 34d45'12.34 for parse_degrees		
		$str =~ s/ //g;
		$str =~ s/deg/d/g;
		
		$str =~ /^([0-9:.\260'"d -]*)([NSEW])$/i;
	
		$latlong = &parse_degrees($1);
		if (defined $latlong) {
			$latlong *= (($2 eq "N" || $2 eq "n" || $2 eq "E" || $2 eq "e") ? 1.0 : -1.0);
			return ($latlong);
		}
	}
	
	return undef;
}

# borrowed from http://www.indo.com/distance/dist.pl
# convert a string like 34:45:12.34 or 38:40 or 34.124 or 34d45'12.34"
# or 25° 02' 30" to degrees
# (also handles a leading `-')
sub parse_degrees {
    my($str) = @_;
    my($d,$m,$s,$sign);

    # yeah, this could probably be done with one regexp.
    if ($str =~ /^\s*(-?)([\d.]+)\s*(:|d|\260)\s*([\d.]+)\s*(:|\')\s*([\d.]+)\s*\"?\s*$/) {
	$sign = ($1 eq "-") ? -1.0 : 1.0;
	$d = $2 + 0.0;
	$m = $4 + 0.0;
	$s = $6 + 0.0;
    } elsif ($str =~ /^\s*(-?)([\d.]+)\s*(:|d|\260)\s*([\d.]+)(\')?\s*$/) {
	$sign = ($1 eq "-") ? -1.0 : 1.0;
	$d = $2 + 0.0;
	$m = $4 + 0.0;
	$s = 0.0;
    } elsif ($str =~ /^\s*(-?)([\d.]+)(d|\260)?\s*$/) {
	$sign = ($1 eq "-") ? -1.0 : 1.0;
	$d = $2 + 0.0;
	$m = 0.0;
	$s = 0.0;
    } else {
	die "parse_degrees: can't parse $str\n";
    }
    return ($sign * ($d + ($m / 60.0) + ($s / 3600.0)));
}

sub parse_csv {
	return quotewords(",",0, $_[0]);
}

################### unused stuff below this line #####################

### kml data still to be implemented

#	<Style id="PhotoIconHighlight">
#		<IconStyle>
#			<scale>1.1</scale>
#			<Icon>
#				<href>http://maps.google.com/mapfiles/kml/pal4/icon38.png</href>
#			</Icon>
#		</IconStyle>
#		<LabelStyle>
#			<scale>1.1</scale>
#		</LabelStyle>
#	</Style>
#	<StyleMap id="PhotoIconPair">
#		<Pair>
#			<key>normal</key>
#			<styleUrl>#PhotoIconNormal</styleUrl>
#		</Pair>
#		<Pair>
#			<key>highlight</key>
#			<styleUrl>#PhotoIconHighlight</styleUrl>
#		</Pair>
#	</StyleMap>
#	<Style id="PhotoIconNormal">
#		<IconStyle>
#			<Icon>
#				<href>http://maps.google.com/mapfiles/kml/pal4/icon46.png</href>
#			</Icon>
#		</IconStyle>
#	</Style>

#		<Placemark id="g2id_1109">
#			<name>New signs in Sowden Lane (covered up)</name>
#			<Snippet maxLines="2"></Snippet>
#			<description><![CDATA[<a href="http://www.morgad.co.uk/gallery2/v/cycle01/08-signs.JPG.html">
#<img src="http://192.168.1.13/gallery2/main.php?g2_view=core.DownloadItem&g2_itemId=1110"></a><br>]]></description>
#			<LookAt>
#				<longitude>-3.428278</longitude>
#				<latitude>50.643947</latitude>
#				<altitude>0</altitude>
#				<range>3000</range>
#				<tilt>45</tilt>
#				<heading>0</heading>
#			</LookAt>
#			<styleUrl>#PhotoIconPair</styleUrl>
#			<Point>
#			<coordinates>-3.428278,50.643947,0</coordinates>
#			</Point>
#		</Placemark>

__END__

=head1 NAME

csv2kml - convert QstartZ csv file into GoogleEarth kml or gpx file, also geotag photos.

=head1 SYNOPSIS

csv2kml F<csvfilename> [options] 

Options: 
 --outgpxlong F<filename>
 --outgpxshort F<filename>
 --outkmllong F<filename>
 --outkmlshort F<filename>
 --hdop F<number>
 --mindist F<number>
 --dryrun
 --reverse
 --help
 --man

(* = not implemented yet)

=head1 OPTIONS

=over 8

=item --B<outgpxlong> F<filename>

The name of the output gpx file that is only HDOP filtered (default is no output)

=item --B<outgpxshort> F<filename>

The name of the output gpx file that is HDOP filtered, as well as cropped from second to last button presses, and with all points less than MINDIST metres apart removed (default is no output)

=item --B<outkmllong> F<filename>

The name of the output kml file (same data as long gpx file) (default is no output)

=item --B<outkmlshort> F<filename>

The name of the output kml file (same data as short gpx file) (default is no output)

=item --B<hdop> F<number>

Max HDOP value for filtering bad gps data (default 2.5)

=item --B<mindist> F<number>

Value in metres for filtering points too close together on 'short' output (default is 20)

=item --B<dryrun> 

Suppress all output to files, just display messages/statistics on screen

=item --B<reverse> 

Reverses order of output tracks and waypoints so it goes from the end of the csv file to the beginning

=item --B<help>

Print a brief help message and exits.

=item --B<man>

Prints the manual page and exits.

=back


=head1 DESCRIPTION

B<csv2kml> will convert a QstartZ csv file into a GoogleEarth kml file or a gpx file, the first button press
in the csv data is assumed to be simultaneous with first jpeg (when sorted by time)
the offset between these two times used to write accurate
GPS data to EXIF headers of jpegs.

Button presses and jpeg pictures are output to the kml/gpx file as waypoints,
'T' events used to generate ground tracks in kml/gpx file, the second and the last button 
presses set the endpoints of the track. Track is filtered for points less than
L<mindist|"item/mindist number"> Metres apart and all points where HDOP is greater than L<hdop|"item/hdop number"> are totally ignored.

=head1 EXAMPLES

Current example of use - it does not do much yet ...

 $ ./csv2kml newton01.csv --outgpxshort newton-exeter.gpx
 9 jpegs found.
 Ref button press was at index # 8735
 Start button press was at index # 8750
 Stop button press was at index # 10893

 attempting to output to newton-exeter.gpx, from 8750 to 10893
 3 waypoints dumped to gpx file
 12262 read in, 1418 written out
 length of trace is 2143 csv_seconds

 $ ./csv2kml newton01.csv --reverse --outgpxshort exeter-newton.gpx
 9 jpegs found.
 Ref button press was at index # 8735
 Start button press was at index # 8750
 Stop button press was at index # 10893

 attempting to output to exeter-newton.gpx, from 10893 to 8750
 3 waypoints dumped to gpx file
 12262 read in, 1418 written out
 length of trace is 2143 csv_seconds

=head1 DOWNLOAD

changelog is here L<http://www.morgad.co.uk/csv2kml-log.txt>

sourcecode is here L<http://www.morgad.co.uk/csv2kml>

(both files last updated __DATESTAMP__)

=head1 COPYRIGHT

 Copyright 2007, 2008 by David R Morgan
 released under GPL v2, use at your own risk
 patches for corrections/improvements welcome

=cut