[Keith Vetter] -- for another weekend project I needed to get some USGS topographical maps images. The USGS and Microsoft have combined to create 
TerraServer (see http://terraserver.homeadvisor.msn.com) which is a vast, 
public domain, database of the USGS topo maps and aerial images in many 
different resolutions. The only complication in using TerraServer is that it 
works in UTM coordinates as opposed to latitude and longitude. Last
weekend I finally got writing the script that does the conversion into UTM
and grabs the images.

As a side project, I also wrote this little utility to share with others 
the core technology of getting maps off of TerraServer. This is just a 
very basic front-end GUI: you specify latitude and longitude (or UTM if 
you know it), the resolution and type of image and hit the "Get Map" 
button and it will display that image.

Notes about the images: each image is 200x200 pixels; the resolution
is in meters per pixel; not all resolutions are available for all image
types.

For grins, try looking at latitude 32.15112, longitude -110.83032,
resolution 4m/p and as an image.
  

======
#!/bin/sh
# The next line is executed by /bin/sh, but not tcl \
    exec wish $0 ${1+"$@"}

##+##########################################################################
#
# TkTopoMap
#
# Simple front end to Microsoft's terraserver which servers up topo
# maps and images. Has to convert from latitude and longitude into
# Universal Transverse Mercator coordinates. The conversion is performed
# using the WGS-84 datum.
# by Keith Vetter
#
# Revisions:
# KPV Aug 20, 2002 - initial revision
#
##+##########################################################################
#############################################################################

package require Tk
if {[catch {package require Img}]} {
    wm withdraw .
    tk_messageBox -icon error -title TkTopoMap \
        -message "TkTopoMap\n\nError: this program requires the Img package."
    exit
}
package require http 2.0

proc Init {} {
    global env state 

    set state(latitude) ""                      ;# Value we'll map
    set state(longitude) ""                     ;# Value we'll map
    set state(utm) ""                           ;# utm of lat/long
    set state(vscale) "16 m/p"                  ;# Resolution of image
    set state(vtheme) Topo                      ;# Type of image
    set state(fetching) 0

    # Terraserver chunks images so we must adjust utm by this value
    array set state {
        mult,10  200    mult,11  400   mult,12  800   mult,13    1600
        mult,14  3200   mult,15  6400  mult,16  12800 mult,17  25600
        mult,18  51200  mult,19  102400
    }

    # Map resolutions available from the server
    set state(resolutions) [list "1 m/p" "2 m/p" "4 m/p" \
                                "8 m/p" "16 m/p" "32 m/p" \
                                "64 m/p" "128 m/p" \
                                "256 m/p" "512 m/p"]
    set s 9
    foreach r $state(resolutions) {
        set state(res,$r) [incr s]              ;# String => scale
        set state(i,res,$s) $r                  ;# Scale => string
    }
    set state(themes) [list Topo Image]         ;# Type of images available
    set state(theme,Image) 1
    set state(theme,Topo)  2

    set state(url) "http://terraserver-usa.com/tile.ashx"
    append state(url) "?T=\$TTHEME&S=\$SSCALE&X=\$XX&Y=\$YY&Z=\$ZZONE&W=0"
    
    if {[info exists env(temp)]} {              ;# We put image into here
        set state(tempfile) [file join $env(temp) get_topo.temp]
    } elseif {[info exists env(tmp)]} {
        set state(tempfile) [file join $env(tmp) get_topo.temp]
    } elseif {[file isdirectory /tmp]} {
        set state(tempfile) [file join /tmp get_topo.temp]
    } else {
        set state(tempfile) get_topo.temp
    }
    trace variable state w DoTrace
}
##+##########################################################################
# 
# Display
# 
# Sets up our simple dialog
# 
proc Display {} {
    global state

    wm title . TkTopoMap
    label .l
    
    frame .f0
    frame .f1a -relief sunken -bd 2
    frame .f1b -relief sunken -bd 2
    frame .f2 -relief sunken -bd 2
    frame .f2.a ; frame .f2.b
    frame .f3 -relief sunken -bd 2
    pack .f0 -side top -fill x -ipady 3 -ipadx 3
    pack .f3 .f2 -side bottom -fill x -ipady 3 -ipadx 3
    pack .f2.a .f2.b -side left -fill both -expand 1
    pack .f1a .f1b -side left -fill both -ipady 3 -ipadx 3
    
    label .title -text TkTopoMap \
        -font [eval font create [font actual [.l cget -font]] -size 18]
    label .title2 -text "by Keith Vetter" \
        -font [eval font create [font actual [.l cget -font]] -size 10]
    button .b -text "About TkTopoMap" -command about -font [.title2 cget -font]
    pack .title -in .f0 -side top
    pack .b -in .f0 -expand 1 -side left -pady 5
    

    label .latlong -text "Latitude/Longitude" \
        -font [eval font create [font actual [.l cget -font]] -size 10] 

    label .lat -text Latitude
    label .long -text Longitude
    entry .elat -width 15 -justify center -textvariable state(latitude)
    entry .elong -width 15 -justify center -textvariable state(longitude)
    label .spacer -text ""
    button .go1a -text " GetMap " -command GoLatLong
    grid .latlong - -in .f1a
    grid .lat .elat -sticky ew -in .f1a
    grid .long .elong -sticky ew -in .f1a
    grid .spacer - -sticky ew -in .f1a
    grid .go1a - -in .f1a -pady 10

    label .utm -text "UTM" \
        -font [eval font create [font actual [.l cget -font]] -size 10] 
    label .northing -text "Northing"
    entry .enorthing -width 15 -justify center -textvariable state(northing)
    label .easting -text Easting
    entry .eeasting -width 15 -justify center -textvariable state(easting)
    label .zone -text Zone
    entry .ezone -width 15 -justify center -textvariable state(zone)
    button .go1b -text " GetMap " -command GoUTM

    grid .utm - -in .f1b
    grid .northing .enorthing -sticky ew -in .f1b
    grid .easting .eeasting -sticky ew -in .f1b
    grid .zone .ezone -sticky ew -in .f1b
    grid .go1b - -in .f1b -pady 10
    
    label .sc  -text "Resolution"
    eval tk_optionMenu .esc state(vscale) $state(resolutions)
    .esc config -width 8
    label .theme -text Theme
    eval tk_optionMenu .etheme state(vtheme) $state(themes)
    .etheme config -width 8
    pack .sc .esc -side left -in .f2.a
    pack .theme .etheme -side left -in .f2.b
    
    label .terra -text "TerraServer:" -justify right
    label .eterra -textvariable state(terra) -width 30 -relief sunken
    canvas .c -bd 1 -relief ridge -width 200 -height 200
    .c create text 100 100 -anchor c -text "no image" -tag what
    image create photo map
    .c create image 0 0 -anchor nw -image map
    
    grid .terra .eterra -sticky ew -in .f3
    grid .c - -in .f3
}
##+##########################################################################
# 
# DoTrace
# 
# Traces the state variable so we know when to enable certain buttons.
# 
proc DoTrace {arg1 arg2 op} {
    global state

    if {$state(fetching)} {                     ;# Always off while fetching
        .go1a config -state disabled
        .go1b config -state disabled
        return
    }

    set new disabled                            ;# Check the lat/long button
    if {$state(latitude) != "" && $state(longitude) != ""} {
        set new normal
    }
    .go1a config -state $new
    
    set new disabled                            ;# Check the utm button
    if {$state(northing) != "" && $state(easting) != "" && $state(zone) != ""} {
        set new normal
    }
    .go1b config -state $new
}
##+##########################################################################
# 
# GoUTM
# 
# Handles getting map from UTM coordinates.
# 
proc GoUTM {} {
    global state
    set emsg ""
    
    if {! [string is double $state(northing)]} {
        set emsg "Can't figure out northing: '$state(northing)'"
    } elseif {! [string is double $state(easting)]} {
        set emsg "Can't figure out easting: '$state(easting)'"
    } elseif {! [string is double $state(zone)]} {
        set emsg "Can't figure out zone: '$state(zone)'"
    }
    if {$emsg != ""} {
        tk_messageBox -icon error -message $emsg
        return
    }
    set state(latitude) ""
    set state(longitude) ""
    GetMap
}
##+##########################################################################
# 
# GoLatLong
# 
# Handles getting map from latitude and longitude.
# 
proc GoLatLong {} {
    global state

    # Extract the lat/long from the entry boxes
    set state(latitude) [string trim $state(latitude)]
    set lat $state(latitude)
    regsub -nocase {N *(.*)} $lat {\1} lat
    regsub -nocase {S *(.*)} $lat {-\1} lat
    if {[regexp {^[\d.]+$} $lat]} {
        ;
    } elseif {[regexp {^([\d.]+)\s+([\d.]+)\s+([\d.]+)$} $lat => a b c]} {
        set lat [expr {$a + $b/60.0 + $c/60.0/60.0}]
    } else {
        set emsg "Can't figure out $state(latitude)"
        tk_messageBox -icon error -message $emsg
        return
    }
    
    set state(longitude) [string trim $state(longitude)]
    set long $state(longitude)
    regsub -nocase {W *(.*)} $long {-\1} long
    regsub -nocase {E *(.*)} $long {\1} long
    if {[regexp {^-?[\d.]+$} $long]} {
        ;
    } elseif {[regexp {^(-?)([\d.]+)\s+([\d.]+)\s+([\d.]+)$} $long \
                   => west a b c]} {
        set long [expr {$a + $b/60.0 + $c/60.0/60.0}]
        if {$west != ""} { set long [expr {-1 * $long}]}
    } else {
        set emsg "Can't figure out $state(longitude)"
        tk_messageBox -icon error -message $emsg
        return
    }

    # Convert first to UTM coordinates
    set state(utm) [ll2utm $lat $long]
    foreach {state(northing) state(easting) state(zone)} $state(utm) break

    GetMap
}
##+##########################################################################
# 
# GetMap
# 
# Fetches and displays the map for this utm from the terraserver. We
# fetch the image straight into a temp file for efficiency.
# 
proc GetMap {} {
    global state

    # Adjust scale for this image's theme
    set state(scale) $state(res,$state(vscale))
    set state(theme) $state(theme,$state(vtheme))
    if {$state(theme) == 2 && $state(scale) == 10} { set state(scale) 11 }
    if {$state(theme) == 1 && $state(scale) >= 17} { set state(scale) 16 }
    set state(vscale) $state(i,res,$state(scale))
    
    catch {image delete map}
    .c itemconfig what -text "fetching image"
    set state(fetching) 1
    
    # Chunk utm for
    set mult $state(mult,$state(scale))
    set XX [expr {int($state(easting)  / $mult)}]
    set YY [expr {int($state(northing) / $mult)}]
    set SSCALE $state(scale)
    set TTHEME $state(theme)
    set ZZONE $state(zone)
    set state(xurl) [subst $state(url)]
    set state(terra) "X=$XX  Y=$YY  Z=$ZZONE  MULT=$mult"
    after 1 GetMap2
}
##+##########################################################################
# 
# GetMap2
# 
# This gets the actual map. It's a separate procedure so we can call it
# as an after task.
# 
proc GetMap2 {} {    
    global state
    ;# Fetch the image from terraserver into our tempfile
    set f [open $state(tempfile) w]
    set token [::http::geturl $state(xurl) -channel $f]
    ::http::wait $token
    close $f

    ;# Display the image
    image create photo map -file $state(tempfile)
    file delete $state(tempfile)
    set state(fetching) 0
}
##+##########################################################################
# 
# ll2utm
# 
# Convert latitude and longitude into Universal Transverse Mercator (UTM)
# coordinates. Lots of fun math which I got off the web.
# 
proc ll2utm {latitude longitude} {
    set PI [expr {atan(1) * 4}]
    set K0 0.9996

    # WGS-84
    set er 6378137                              ;# EquatorialRadius
    set es2 0.00669438                          ;# EccentricitySquared
    set es4 [expr {$es2 * $es2}]
    set es6 [expr {$es2 * $es2 * $es2}]

    # Must be in the range -180 <= long < 180
    while {$longitude < -180} { set longitude [expr {$longitude + 360}]}
    while {$longitude >= 180}  { set longitude [expr {$longitude - 360}]}

    # Now convert
    set lat_rad [expr {$latitude * $PI / 180.0}]
    set long_rad [expr {$longitude * $PI / 180.0}]
    
    set zone [expr {int(($longitude + 180) / 6) + 1}]
    if {$latitude >= 56.0 && $latitude < 64.0 &&
        $longitude >= 3.0 && $longitude < 12.0} {
        $zone = 32
    }
    if { $latitude >= 72.0 && $latitude < 84.0 } {
        if { $longitude >= 0.0  && $longitude <  9.0 } {$zone = 31;}
        if { $longitude >= 9.0  && $longitude < 21.0 } {$zone = 33;}
        if { $longitude >= 21.0 && $longitude < 33.0 } {$zone = 35;}
        if { $longitude >= 33.0 && $longitude < 42.0 } {$zone = 37;}
    }
    # +3 puts origin in middle of zone    
    set long_origin [expr {( $zone - 1 ) * 6 - 180 + 3}]
    set long_origin_rad [expr {$long_origin * $PI / 180.0}]
    set eccPrimeSquared [expr {$es2 / ( 1.0 - $es2 )}]
    set N [expr {$er / sqrt( 1.0 - $es2 * sin( $lat_rad ) * sin( $lat_rad ) )}]
    set T [expr {tan( $lat_rad ) * tan( $lat_rad )}]
    set C [expr {$eccPrimeSquared * cos( $lat_rad ) * cos( $lat_rad )}]
    set A [expr {cos( $lat_rad ) * ( $long_rad - $long_origin_rad )}]
    set M [expr { $er * ( \
       (1.0 - $es2 / 4 - 3 * $es4 / 64 - 5 * $es6 / 256) * $lat_rad \
       - (3 * $es2 / 8 + 3 * $es4 / 32 + 45 * $es6 / 1024) * sin(2 * $lat_rad) \
       + (15 * $es4 / 256 + 45 * $es6 / 1024 )             * sin(4 * $lat_rad) \
       - (35 * $es6 / 3072 )                               * sin(6 * $lat_rad) \
                              )}]
    set easting [expr {$K0 * $N * ( $A + ( 1 - $T + $C ) * $A * $A * $A / 6 \
       + ( 5 - 18 * $T + $T * $T + 72 * $C - 58 * $eccPrimeSquared ) * \
       $A * $A * $A * $A * $A / 120 ) + 500000.0}]
    set northing [expr {$K0 * ( $M + $N * tan( $lat_rad ) * \
       ( $A * $A / 2 + ( 5 - $T + 9 * $C + 4 * $C * $C ) * \
       $A * $A * $A * $A / 24 + ( 61 - 58 * $T + $T * $T + \
       600 * $C - 330 * $eccPrimeSquared ) * \
       $A * $A * $A * $A * $A * $A / 720 ) )}]

    if {$latitude < 0} {  ;# 1e7 meter offset for southern hemisphere  
        set northing [expr {$northing + 10000000.0}]
    }
    
    set northing [expr {int($northing)}]
    set easting [expr {int($easting)}]
    return [list $northing $easting $zone]
}
##+##########################################################################
# 
# about
# 
# Your basic about box.
# 
proc about {} {
    set msg ""
    append msg "TkTopoMap\nby Keith Vetter\n\n"
    append msg "TkTopoMap is a simple front end for TerraServer, Microsoft's\n"
    append msg "database of USGS aerial images and topographic maps.\n"
    append msg "This data is free to everyone.\n\n"

    append msg "Fetching images from TerraServer is straightforward except\n"
    append msg "for the conversion from latitude and longitude into\n"
    append msg "Universal Transverse Mercator (UTM) coordinates. The\n"
    append msg "conversion is done using the WGS-84 datum.\n\n"
    append msg "For more details see http://terraserver.homeadvisor.msn.com."

    tk_messageBox -message $msg -title TkTopoMap
    
}
################################################################
################################################################
################################################################

Init
Display

set state(latitude) 40.069
set state(longitude) -82.517 
set state(latitude) "37 48 9"
set state(longitude) "-122 13 29"

if {0} {                                        ;# Beale Airforce base
    set state(latitude) 32.15112
    set state(longitude) -110.83032
    set state(vscale) "4 m/p"   
    set state(vtheme) Image
}
======