#!/usr/bin/perl -w
# NAME: wsg842mercator.pl
# AIM: Convert a WSG84 lon,lat pair to a a spheical mercator point
# 31/10/2011 geoff mclane http://geoffair.net/mperl
# from : http://www.mail-archive.com/gdal-dev@lists.osgeo.org/msg03050.html
# /* spherical mercator for Google, VE, Yahoo etc
#  * epsg:900913 R= 6378137 
#  * x = longitude
#  * y= R*ln(tan(pi/4 + latitude/2)
#  */
# Also see : http://code.google.com/p/gepy/source/browse/trunk/tile.py
# Coordinate extent of Earth in EPSG:900913:
# [-20037508.342789244, -20037508.342789244,  20037508.342789244, 20037508.342789244]
# The constant 20037508.342789244 comes from the circumference of the Earth  in meters, 
# which is 40 thousand kilometers, the coordinate origin is in the  middle of extent.  
# In fact you can calculate the constant as: 2 * math.pi * 6378137 / 2.0  
# Polar areas with abs(latitude) bigger then 85.05112878 are clipped off.    
#What are zoom level constants (pixels/meter) for pyramid with EPSG:900913?  
# whole region is on top of pyramid (zoom=0) covered by 256x256 pixels tile,  
# every lower zoom level resolution is always divided by two  
# initialResolution = 20037508.342789244 * 2 / 256 = 156543.03392804062    

use strict;
use warnings;
use File::Basename;  # split path ($name,$dir,$ext) = fileparse($file [, qr/\.[^.]*/] )
use Cwd;
use Math::Trig;
my $perl_dir = 'C:\GTools\perl';
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 = $perl_dir."\\temp.$pgmname.txt";
open_log($outfile);

# user variables
my $VERS = "0.0.1 2011-10-31";
my $load_log = 0;
my $in_file = '';
my $verbosity = 0;
my $reverse = 0;
my $debug_on = 0;
my $def_file = 'def_file';

my $PI = 3.1415926535897932384626433832795029;
my $D2R = $PI / 180;
my $R2D = 180 / $PI;
my $RADIUS = 6378137;

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

my $u_lon = 120;
my $u_lat = 30;

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" );
    }
}

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 convert_lon_lat($$) {
    my ($lon,$lat) = @_;
    my $x = $RADIUS * $D2R * $lon;
    my $y = $RADIUS * log(tan($D2R*(45 + ($lat / 2))));
    return $x,$y;
}

sub convert_mx_my($$) {
    my ($x,$y) = @_;
    my $lon = ($x / $RADIUS) * $R2D;
    my $lat = ((atan( exp( $y / $RADIUS ) ) * $R2D) - 45) * 2;
    return $lon,$lat;
}

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 in_world_range($$) {
    my ($lon,$lat) = @_;
    my $bad = 0;
    if (($lon < -180)||($lon > 180)) {
        prt("Longitude $lon OUT OF RANGE -180 to 180 degrees!\n");
        $bad++;
    }
    if (($lat < -90)||($lat > 90)) {
        prt("Latitude $lat OUT OF RANGE -90 to 90 degrees!\n");
        $bad++;
    }
    if ($bad) {
        pgm_exit(1,"Out of range input.\n");
    }
}

sub show_extents()
{
    prt("Coordinate extent of Earth in EPSG:900913\n");
    # [-20037508.342789244, -20037508.342789244,  20037508.342789244, 20037508.342789244]
    my $mx1 = -20037508.342789244;
    my $my1 = -20037508.342789244;
    my ($lon1,$lat1) = convert_mx_my($mx1,$my1);
    prt("Mercator x = $mx1, y = $my1 == lon $lon1, lat $lat1\n");
    my $mx2 = 20037508.342789244;
    my $my2 = 20037508.342789244;
    my ($lon2,$lat2) = convert_mx_my($mx2,$my2);
    prt("Mercator x =  $mx2, y =  $my2 == lon  $lon2, lat  $lat2\n");
    if (VERB9()) {
        prt("And testing the reversal...\n");
        ($mx1,$my1) = convert_lon_lat($lon1,$lat1);
        prt("Mercator x = $mx1, y = $my1 == lon $lon1, lat $lat1\n");
        ($mx2,$my2) = convert_lon_lat($lon2,$lat2);
        prt("Mercator x =  $mx2, y =  $my2 == lon  $lon2, lat  $lat2\n");
    }
    prt("\n");
}

#########################################
### MAIN ###
parse_args(@ARGV);
##prt( "$pgmname: in [$cwd]: Hello, World...\n" );
##process_in_file($in_file);
in_world_range($u_lon,$u_lat);
show_extents() if (VERB5());
prt( "Convert lon = $u_lon, lat = $u_lat\n");
my ($mx,$my) = convert_lon_lat($u_lon,$u_lat);
prt( "Mercator x = $mx, y = $my\n");
my ($nlon,$nlat) = convert_mx_my( $mx, $my );
prt( "Back to lon = $nlon, lat = $nlat\n");

pgm_exit(0,"");
########################################
sub give_help {
    prt("$pgmname: version $VERS\n");
    prt("Usage: $pgmname [options] longitude latitude\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(" --reverse     (-r) = Reverse lon and lat.\n");
}
sub need_arg {
    my ($arg,@av) = @_;
    pgm_exit(1,"ERROR: [$arg] must have a following argument!\n") if (!@av);
}

sub parse_args {
    my (@av) = @_;
    my ($arg,$sarg,$cnt,@arr);
    $cnt = 0;
    while (@av) {
        $arg = $av[0];
        if (($arg =~ /^-/) && !($arg =~ /^-(\d|\.)+$/)) {
            $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);
                    }
                }
                prt("Verbosity = $verbosity\n") if (VERB1());
            } elsif ($sarg =~ /^r/) {
                $reverse = 1;
            } else {
                pgm_exit(1,"ERROR: Invalid argument [$arg]! Try -?\n");
            }
        } else {
            if ($cnt == 0) {
                if ($arg =~ /,/) {
                    @arr = split(",",$arg);
                    if (scalar @arr == 2) {
                        $u_lon = $arr[0];
                        $u_lat = $arr[1];
                    } else {
                        pgm_exit(1,"ERROR: Invalid argument [$arg]! Did not split into 2, lon,lat\n");
                    }
                    $cnt = 2;
                } else {
                    $u_lon = $arg;
                    $cnt++;
                }
            } elsif ($cnt == 1) {
                $u_lat = $arg;
                $cnt++;
            } else {
                pgm_exit(1,"ERROR: Invalid argument [$arg]! Already have lon = $u_lon and lat = $u_lat!\n");
            }
        }
        shift @av;
    }

    if ($cnt != 2) {
        pgm_exit(1,"ERROR: lon and lat NOT found in command!\n");
    }
    if ($reverse) {
        my $arg = $u_lon;
        $u_lon = $u_lat;
        $u_lat = $arg;
    }
}

# eof - wsg842mercator.pl
