#!/usr/bin/perl -w
# NAME: maptile.pl
# AIM: Convert a lat long zoom to a tile number
# from : http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
#
# Directory structure - like 
# OSM 'standard' style	 http://[abc].tile.openstreetmap.org/zoom/x/y.png	 0-19
#
# Zoom Coverage     Tiles
# 0    1 tile world 1
# 1    2 X 2        4
# 2    4 x 4        16
# n    2^n x 2^n    2^2n
# 12   4096 x 4096  16,777,216
# 16                  4,294,967,296 = 2^32 
# 17                 17,179,868,184
# 18                 68,719,476,736
# 19   Max. mapnik  274.877.906,944
# X goes from 0 (left edge is 180 deg W) to 2^zoom - 1 (right edge is 180 deg E
# Y goes from 0 (top edge is 85.0511 to 2^zoom - 1 (bottom edge is 85.0511 deg S) in a Mercator projection
# For the curious, the num 85.0511 is the result of arctan(sinh(Pi)). By using this bound, the entire map 
# becomes a (very large) square. See also the Osmarender bug.
#
# 17/02/2014 geoff mclane http://geoffair.net/mperl
use strict;
use warnings;
use File::Basename;  # split path ($name,$dir,$ext) = fileparse($file [, qr/\.[^.]*/] )
use Math::Trig; # qw(great_circle_distance great_circle_direction deg2rad rad2deg);
##use POSIX;
use POSIX qw/floor/;
use Cwd;
my $os = $^O;
my $perl_dir = '/home/geoff/bin';
my $PATH_SEP = '/';
my $temp_dir = '/tmp';
if ($os =~ /win/i) {
    $perl_dir = 'C:\GTools\perl';
    $temp_dir = $perl_dir;
    $PATH_SEP = "\\";
}
unshift(@INC, $perl_dir);
require 'lib_utils.pl' or die "Unable to load 'lib_utils.pl' Check paths in \@INC...\n";
# log file stuff
our ($LF);
my $pgmname = $0;
if ($pgmname =~ /(\\|\/)/) {
    my @tmpsp = split(/(\\|\/)/,$pgmname);
    $pgmname = $tmpsp[-1];
}
my $outfile = $temp_dir.$PATH_SEP."temp.$pgmname.txt";
open_log($outfile);

# user variables
my $VERS = "0.0.2 2014-01-13";
my $load_log = 0;
my $in_lat = '';
my $in_lon = '';
my $in_zoom = '';
my $verbosity = 0;
my $out_file = '';

# ### DEBUG ###
my $debug_on = 0;
my $def_lat = 37.78756; #//38.7;
my $def_lon = -122.0;
my $def_zoom = 9;  # //11;  // hard to choose then best init zoom

### program variables
my @warnings = ();
my $cwd = cwd();

sub VERB1() { return $verbosity >= 1; }
sub VERB2() { return $verbosity >= 2; }
sub VERB5() { return $verbosity >= 5; }
sub VERB9() { return $verbosity >= 9; }

sub show_warnings($) {
    my ($val) = @_;
    if (@warnings) {
        prt( "\nGot ".scalar @warnings." WARNINGS...\n" );
        foreach my $itm (@warnings) {
           prt("$itm\n");
        }
        prt("\n");
    } else {
        prt( "\nNo warnings issued.\n\n" ) if (VERB9());
    }
}

sub pgm_exit($$) {
    my ($val,$msg) = @_;
    if (length($msg)) {
        $msg .= "\n" if (!($msg =~ /\n$/));
        prt($msg);
    }
    show_warnings($val);
    close_log($outfile,$load_log);
    exit($val);
}


sub prtw($) {
   my ($tx) = shift;
   $tx =~ s/\n$//;
   prt("$tx\n");
   push(@warnings,$tx);
}

sub process_in_file($) {
    my ($inf) = @_;
    if (! open INF, "<$inf") {
        pgm_exit(1,"ERROR: Unable to open file [$inf]\n"); 
    }
    my @lines = <INF>;
    close INF;
    my $lncnt = scalar @lines;
    prt("Processing $lncnt lines, from [$inf]...\n");
    my ($line,$inc,$lnn);
    $lnn = 0;
    foreach $line (@lines) {
        chomp $line;
        $lnn++;
        if ($line =~ /\s*#\s*include\s+(.+)$/) {
            $inc = $1;
            prt("$lnn: $inc\n");
        }
    }
}

sub mlog_10 {
  my $n = shift;
  return log($n)/log(10);
}


# from Theo HTML
#// The math
#// http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Tile_servers
#	function lon2tile( lon, zoom ) {
#		return Math.floor( ( lon + 180 ) / 360 * pow( 2, zoom ) );
#	}
#	function lat2tile( lat, zoom ) {
#		return Math.floor(( 1 - Math.log( Math.tan( lat * pi / 180) + 1 / cos( lat * pi / 180)) / pi )/2 * pow(2, zoom) );
#	}
#	function tile2lon( x, z ) {
#		return ( x / pow( 2, z ) * 360 - 180 );
#	}
#	function tile2lat( y, z ) {
#		var n = pi - 2 * pi * y / pow( 2, z );
#		return 180 / pi * Math.atan( 0.5 * ( Math.exp( n ) - Math.exp( -n ) ));
#	}
sub lat2tilef($$) {
    my ($lat,$zoom) = @_;
    my $lat_rad = deg2rad($lat);
    my $n = 2 ** $zoom;
    my $ytile = floor((1.0 - mlog_10(tan($lat_rad) + (1 / cos($lat_rad))) / pi) / 2.0 * $n);
    return $ytile;
}

sub lon2tilef($$) {
    my ($lon,$zoom) = @_;
    my $n = 2 ** $zoom;
    my $xtile = floor($n * (($lon + 180) / 360));
    return $xtile;
}


# Python
# Lon./lat. to tile numbers
# import math
# def deg2num(lat_deg, lon_deg, zoom):
#  lat_rad = math.radians(lat_deg)
#  n = 2.0 ** zoom
#  xtile = int((lon_deg + 180.0) / 360.0 * n)
#  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
#  return (xtile, ytile)
sub lat2tile($$) {
    my ($lat,$zoom) = @_;
    my $lat_rad = deg2rad($lat);
    my $n = 2 ** $zoom;
    my $ytile = int((1.0 - mlog_10(tan($lat_rad) + (1 / cos($lat_rad))) / pi) / 2.0 * $n);
    return $ytile;
}

sub lon2tile($$) {
    my ($lon,$zoom) = @_;
    my $n = 2 ** $zoom;
    my $xtile = int($n * (($lon + 180) / 360));
    return $xtile;
}


sub show_mappy() {
    prt("Converting lat $in_lat, lon $in_lon, zoom = $in_zoom to mappy\n");
    my $lat_rad = deg2rad($in_lat);
    my $n = 2 ** $in_zoom;
    my $xtile = int($n * (($in_lon + 180) / 360));
    my $ytile = int((1.0 - mlog_10(tan($lat_rad) + (1 / cos($lat_rad))) / pi) / 2.0 * $n);
    my $x = lon2tile($in_lon,$in_zoom);
    my $y = lat2tile($in_lat,$in_zoom);
    prt("xtile=$xtile, ytile=$ytile ($x/$y)\n");
}

sub test() {
    my $z = 7;
    my ($x,$y,$mx,$my,$tx,$ty,$cnt,$cx,$cy);
    my %xdupes = ();
    prt("Running test...\n");
    #for ($x = 0; $x < 180; $x += 0.1) {
    #    for ($y = 0; $y < 90; $y += 0.1) {
    my $xmin = 99999999999;
    my $xmax = -99999999999;
    for ($x = -179; $x < 180; $x += 0.1) {
    #    for ($y = -89; $y < 90; $y += 0.1) {
        for ($y = -85; $y < 85; $y += 0.1) {
            $mx = lon2tile($x,$z);
            $my = lat2tile($y,$z);
            $tx = lon2tilef($x,$z);
            $ty = lat2tilef($y,$z);
            #if (($mx != $tx)||($my != $ty)) {
            #    prt("lon/lat/zoom $x/$y/$z mx=$mx my=$my tx=$tx ty=$ty\n");
            #}
            if ($mx != $tx) {
                prt("lon/lat/zoom $x/$y/$z DIFF mx/tx=$mx/$tx\n");
            }
            if ($my != $ty) {
                prt("lon/lat/zoom $x/$y/$z DIFF my/ty=$my/$ty\n");
            }
            ###prt("lon/lat/zoom $x/$y/$z mx=$mx my=$my\n");
            if (!defined $xdupes{$mx}) {
                $xdupes{$mx} = 1;
                $cnt++;
                $cx = sprintf("%.2f",$x);
                $cy = sprintf("%.2f",$y);
                prt("lon/lat/zoom $cx/$cy/$z mx=$mx my=$my tx=$tx ty=$ty\n");

            }
            $xmax = $mx if ($mx > $xmax);
            $xmin = $mx if ($mx < $xmin);
        }
    }
    # at zoom level 19 that is 3591 directories...
    prt("For zoom level $z, that is $cnt x directories... x min $xmin, max $xmax\n");
    $load_log = 1;
    pgm_exit(0,"");

}
sub test2() {
    my $val = sqrt(274877906944);
    prt("sqrt(274877906944) = $val x $val\n");

    my ($z);
    for ($z = 0; $z < 20; $z++) {
        $val = 2 ** $z;
        prt("$z = $val x $val\n");
    }

    exit(1);
}
test();

#########################################
### MAIN ###
parse_args(@ARGV);
#process_in_file($in_lat);
show_mappy();
pgm_exit(0,"");
########################################

sub need_arg {
    my ($arg,@av) = @_;
    pgm_exit(1,"ERROR: [$arg] must have a following argument!\n") if (!@av);
}

sub is_a_double($) {
    my $val = shift;
    return 1 if ($val =~ /^-?\d+\.?\d*$/); # { print "is a real number\n" }
    return 1 if ($val =~ /^\d+$/); # { print "is a whole number\n" }
    return 1 if ($val =~ /^-?\d+$/); # { print "is an integer\n" }
    return 1 if ($val =~ /^[+-]?\d+$/); # { print "is a +/- integer\n" }
    return 1 if ($val =~ /^-?(?:\d+(?:\.\d*)?|\.\d+)$/); # { print "is a decimal number\n" }
    return 1 if ($val =~ /^([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))?$/); # { print "a C float\n" }
    return 0 if ($val =~ /\D/); # { print "has nondigits\n" }
    return 0;
}

sub parse_args {
    my (@av) = @_;
    my ($arg,$sarg);
    my $verb = VERB2();
    my $cnt = 0;
    while (@av) {
        $arg = $av[0];
        if (($arg =~ /^-/)&&(!is_a_double($arg))) {
            $sarg = substr($arg,1);
            $sarg = substr($sarg,1) while ($sarg =~ /^-/);
            if (($sarg =~ /^h/i)||($sarg eq '?')) {
                give_help();
                pgm_exit(0,"Help exit(0)");
            } elsif ($sarg =~ /^v/) {
                if ($sarg =~ /^v.*(\d+)$/) {
                    $verbosity = $1;
                } else {
                    while ($sarg =~ /^v/) {
                        $verbosity++;
                        $sarg = substr($sarg,1);
                    }
                }
                $verb = VERB2();
                prt("Verbosity = $verbosity\n") if ($verb);
            } elsif ($sarg =~ /^l/) {
                if ($sarg =~ /^ll/) {
                    $load_log = 2;
                } else {
                    $load_log = 1;
                }
                prt("Set to load log at end. ($load_log)\n") if ($verb);
            } elsif ($sarg =~ /^o/) {
                need_arg(@av);
                shift @av;
                $sarg = $av[0];
                $out_file = $sarg;
                prt("Set out file to [$out_file].\n") if ($verb);
            } else {
                pgm_exit(1,"ERROR: Invalid argument [$arg]! Try -?\n");
            }
        } else {
            if ($cnt == 0) {
                $in_lat = $arg;
                prt("Set latitude input to [$in_lat]\n") if ($verb);
            } elsif ($cnt == 1) {
                $in_lon = $arg;
                prt("Set longitude input to [$in_lon]\n") if ($verb);
            } elsif ($cnt == 2) {
                $in_zoom = $arg;
                prt("Set zomm input to [$in_zoom]\n") if ($verb);
            } else {
                pgm_exit(1,"What is this? $arg! Already have lat $in_lat, lon $in_lon, zoom = $in_zoom\n");
            }
            $cnt++;
        }
        shift @av;
    }

    if ($debug_on) {
        prtw("WARNING: DEBUG is ON!\n");
        if (length($in_lat) ==  0) {
            $in_lat = $def_lat;
            $in_lon = $def_lon;
            $in_zoom = $def_zoom;
            prt("Set DEFAULTS lat $in_lat, lon $in_lon, zoom = $in_zoom\n");
        }
    }
    if (length($in_lat) ==  0) {
        pgm_exit(1,"ERROR: No latitude found in command!\n");
    }
    if (length($in_lon) ==  0) {
        pgm_exit(1,"ERROR: No longitude found in command!\n");
    }
    if (length($in_zoom) ==  0) {
        pgm_exit(1,"ERROR: No zoom found in command!\n");
    }
}

sub give_help {
    prt("$pgmname: version $VERS\n");
    prt("Usage: $pgmname [options] latitude longitude zoom\n");
    prt("Options:\n");
    prt(" --help  (-h or -?) = This help, and exit 0.\n");
    prt(" --verb[n]     (-v) = Bump [or set] verbosity. def=$verbosity\n");
    prt(" --load        (-l) = Load LOG at end. ($outfile)\n");
    prt(" --out <file>  (-o) = Write output to this file.\n");
}

# eof - template.pl