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
}