Universal Transverse Mercator

Keith Vetter 2004-04-15 : Universal Transverse Mercator (UTM) is another coordinate system for the earth that has some nice advantages over the more familiar latitude and longitude especially at large scales. The TerraServer online database of USGS topo maps and aerial photographs is based on UTM coordinates (see TkTopomap).

The basic idea behind UTM is that the world is divided into 60 north-south zones. Within each zone, a coordinate is measured in meters north from the equator--northing--and meters east from a central meridian--easting. For gory details see http://mac.usgs.gov/mac/isb/pubs/factsheets/fs07701.html or [broken www.nps.gov/prwi/readutm.htm]: try [L1 ].

Below are two routines to convert from UTM into lat/lon and vice versa. It uses the WGS-84 datum but a table is provided if you want to use other ones. As usual a simple demo is also provided (warning, absolutely no error checking is done.)

##+############################################################# #############
#
# ll2utm -- Converts 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)}]
    if {$latitude > 84.0 || $latitude < -80.0} {
        set letter "Z"
    } else {
        set l [expr {int(($latitude + 80) / 8.0)}]
        set letter [string index "CDEFGHJKLMNPQRSTUVWXX" $l]
    }
    
    return [list $northing $easting $zone $letter]
}
##+##########################################################################
#
# utm2ll -- Converts Universal Transverse Mercator (UTM) into
# latitude and longitude coordinates. More fun math which I got off
# the web.
#
proc utm2ll {northing easting zone {letter S}} {
    set PI [expr {atan(1) * 4}]
    set K0 0.9996
    
    # WGS-84
    set er 6378137                              ;# EquatorialRadius
    set es2 0.00669438                          ;# EccentricitySquared
    set es2x [expr {1.0 - $es2}]
    
    set x [expr {$easting - 500000.0}]
    set northernHemisphere [expr {$letter >= "N"}]
    set y [expr {$northing - ($northernHemisphere ? 0.0 : 10000000.0)}]
    set long_origin [expr {($zone - 1) * 6 - 180 + 3}] ;# +3 puts in middle
    set ep2 [expr {$es2 / $es2x}]
    set e1 [expr {(1.0 - sqrt($es2x)) / (1.0 + sqrt($es2x))}]
    set M [expr {$y / $K0}]
    set mu [expr {$M / ($er * (1.0 - $es2 /4.0 - 3 * $es2 * $es2 /64.0
                               - 5 * $es2 * $es2 * $es2 /256.0))}]
    set phi [expr {$mu + (3 * $e1 / 2 - 27 * $e1 * $e1 * $e1 / 32 ) * sin(2*$mu)
               + (21 * $e1 * $e1 / 16 - 55 * $e1 * $e1 * $e1 * $e1 / 32) *
                   sin(4*$mu ) + (151 * $e1 * $e1 * $e1 / 96 ) * sin(6*$mu)}]
    set N1 [expr {$er / sqrt(1.0 - $es2 * sin($phi) * sin($phi))}]
    set T1 [expr {tan($phi) * tan($phi)}]
    set C1 [expr {$ep2 * cos($phi) * cos($phi)}]
    set R1 [expr {$er * $es2x / pow(1.0 - $es2 * sin($phi) * sin($phi), 1.5)}]
    set D [expr {$x / ($N1 * $K0)}]
    set latitude [expr {$phi - ($N1 * tan($phi) / $R1) * ($D * $D / 2
           - (5 + 3 * $T1 + 10 * $C1 - 4 * $C1 * $C1 - 9 * $ep2) *
           $D * $D * $D * $D / 24
           + (61 + 90 * $T1 + 298 * $C1 + 45 * $T1 * $T1
           - 252 * $ep2 - 3 * $C1 * $C1 ) * $D * $D * $D * $D * $D * $D / 720)}]
    set latitude [expr {$latitude * 180.0 / $PI}]
    set longitude [expr {($D - (1 + 2 * $T1 + $C1) * $D * $D * $D / 6
           + (5 - 2 * $C1 + 28 * $T1 - 3 * $C1 * $C1
           + 8 * $ep2 + 24 * $T1 * $T1) * $D * $D * $D * $D * $D / 120)
           / cos($phi)}]
    set longitude [expr {$long_origin + $longitude * 180.0 / $PI}]
 
    return [list $latitude $longitude]
}
 
## Other datums which could be used here instead:
##
## Name                         EquatorialRadius        EccentricitySquared
## ----                         ----------------        -------------------
## Airy                         6377563                 0.00667054
## Australian National          6378160                 0.006694542
## Bessel 1841                  6377397                 0.006674372
## Bessel 1841 (Nambia)         6377484                 0.006674372
## Clarke 1866                  6378206                 0.006768658
## Clarke 1880                  6378249                 0.006803511
## Everest                      6377276                 0.006637847
## Fischer 1960 (Mercury)       6378166                 0.006693422
## Fischer 1968                 6378150                 0.006693422
## GRS 1967                     6378160                 0.006694605
## GRS 1980                     6378137                 0.00669438
## Helmert 1906                 6378200                 0.006693422
## Hough                        6378270                 0.00672267
## International                6378388                 0.00672267
## Krassovsky                   6378245                 0.006693422
## Modified Airy                6377340                 0.00667054
## Modified Everest             6377304                 0.006637847
## Modified Fischer 1960        6378155                 0.006693422
## South American 1969          6378160                 0.006694542
## WGS 60                       6378165                 0.006693422
## WGS 66                       6378145                 0.006694542
## WGS-72                       6378135                 0.006694318
## WGS-84                       6378137                 0.00669438
 
 
################################################################
#
# DEMO code
#
package require Tk
 
proc DoDisplay {} {
 
    labelframe .ll -text "Lat/Lon" -padx 5
    label .ll.lat -text Latitude:
    entry .ll.elat -textvariable C(lat)
    label .ll.lon -text Longitude:
    entry .ll.elon -textvariable C(lon)
    grid .ll.lat .ll.elat -sticky w
    grid .ll.lon .ll.elon -sticky w
    grid rowconfigure .ll 1000 -weight 1
 
    labelframe .utm -text UTM -padx 5
    label .utm.north -text Northing:
    entry .utm.enorth -textvariable C(north)
    label .utm.east -text Easting:
    entry .utm.eeast -textvariable C(east)
    label .utm.zone -text Zone:
    entry .utm.ezone -textvariable C(zone)
    label .utm.letter -text Letter:
    entry .utm.eletter -textvariable C(letter)
    
    grid .utm.north .utm.enorth -sticky w
    grid .utm.east .utm.eeast -sticky w
    grid .utm.zone .utm.ezone -sticky w
    grid .utm.letter .utm.eletter -sticky w
 
    frame .2
    button .2.utm -text "==>" -command [list Convert 1]
    button .2.ll -text "<==" -command [list Convert 0]
    grid .2.utm -pady 5
    grid .2.ll
 
    grid .ll .2 .utm -sticky ns -pady {5 10} -padx 5
    grid configure .2 -padx 20
    
    image create photo ::img::blank -width 1 -height 1
    button .about -image ::img::blank -highlightthickness 0 -command \
        [list tk_messageBox -message "UTM Converter\nby Keith Vetter\nMarch, 2004"]
    place .about -in . -relx 1 -rely 1 -anchor se
    focus .ll.elat
}
proc lat2dec {value} {
    set cnt [scan "$value 0 0 0" "%g %g %g" l1 l2 l3]
    set dec [expr {abs($l1) + $l2 / 60.0 + $l3 / 3600.0}]
    if {$l1 < 0} {set dec [expr {-1 * $dec}]}
    return $dec
}    
proc Convert {toUTM} {
    global C
    
    if {$toUTM} {
        foreach var {north east zone letter} { set C($var) "" }
        foreach var {lat lon} {
            set C($var) [string trim $C($var)]
            if {$C($var) eq ""} return
            set $var [lat2dec $C($var)]
        }
        set utm [ll2utm $lat $lon]
        foreach {C(north) C(east) C(zone) C(letter)} $utm break
    } else {
        foreach var {lat lon} { set C($var) "" } 
        foreach var {north east zone letter} { 
            set C($var) [string trim $C($var)]
            if {$C($var) eq "" && $var ne "letter"} return
        }
        set ll [utm2ll $C(north) $C(east) $C(zone) $C(letter)]
        foreach {C(lat) C(lon)} $ll break
    }
}
DoDisplay
set C(lat) "40 4 4.5"
set C(lon) [lat2dec "-82 31 12.6"]

PWE: The first line of ll2utm proc

 if {$longitude > 0} {set longitude [expr {-1 * $longitude}]} 

makes this code correct for only the western hemisphere. Probably more correct is to remove the first line and enter western hemisphere longitudes as negative numbers (as seems to be the normal convention). fixed


PD: These formulae appear to be the same as in Snyder http://onlinepubs.er.usgs.gov/djvu/PP/PP_1395.pdf - see formulae 8-9 to 8-25. This is very slow. A better solution would involve an extension to use the proj.4 library http://trac.osgeo.org/proj/wiki/ProjAPI


Screenshots


uniquename 2013aug18

The screenshot above is on an 'external' site (flickr) and is likely to go dead. Also the image is too small to read the text. Here is another screenshot of the UTM GUI --- 'locally stored', on the disk drives of the server of this wiki --- a larger image.

vetter_UniversalTransverseMercator_wiki13044_screenshot_616x142.jpg

When the GUI first came up, the entry fields on the right were blank. A click on the right-arrow caused those entry fields to be filled.

That little dimple on the lower right corner of the GUI is like the 'About' button on many of Vetter's other GUI's. If you click on it, a popup shows the text 'UTM Converter, by Keith Vetter, 2004'. You can see the definition of that button in the code above.


Related topics: geodesy