Generated: Sun Apr 15 11:46:52 2012 from wsg842mercator.pl 2011/11/04 8 KB.
#!/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