Convert Latitude/Longitude to UTM

Converts Latitude and Longitude to UTM coordinates. The first step in using IGOR as a GPS utility.
Based on cpp code found on the internet by Chuck Gantz

#pragma rtGlobals=1     // Use modern global access method.
CONSTANT  DEG2RAD = 0.01745329;
CONSTANT  RAD2DEG = 57.29577951;
Structure GPSpos
    variable UTMNorthing
    variable UTMEasting
    string UTMzone
Endstructure

Function LLtoUTM( ReferenceEllipsoid,   Lat,  Long, gps)
//converts lat/long to UTM coords.  Equations from USGS Bulletin 1532
//East Longitudes are positive, West longitudes are negative.
//North latitudes are positive, South latitudes are negative
//Lat and Long are in decimal degrees
    //Written by Chuck Gantz- chuck.gantz@globalstar.com
    variable ReferenceEllipsoid, Lat, Long
    struct GPSpos &gps
   
    variable UTMEasting,  UTMNorthing
    string UTMZone
   
    make/o/d/n=(23) EquatorialRadius, SquareOfEccentricity
     make/o/n=23/T EllipsoidName
     EquatorialRadius={6377563,     6378160, 6377397, 6377484, 6378206, 6378249, 6377276, 6378166, 6378150, 6378160, 6378137, 6378200, 6378270, 6378388, 6378245, 6377340, 6377304, 6378155, 6378160, 6378165, 6378145, 6378135, 6378137}
    SquareOfEccentricity={.00667054,.006694542,.006674372,.006674372,.006768658,.006803511,.006637847,.006693422,.006693422,.006694605,.00669438,.006693422,.00672267,.00672267,.006693422,.00667054,.006637847,.006693422,.006694542,.006693422,.006694542,.006694318,.00669438}
    EllipsoidName={"Airy","Australian National","Bessel 1841","Bessel 1841 (Nambia) ","Clarke 1866","Clarke 1880","Everest","Fischer 1960 (Mercury) ","Fischer 1968","GRS 1967", "GRS 1980", "Helmert 1906","Hough",  "International", "Krassovsky", "Modified Airy", "Modified Everest",    "Modified Fischer 1960", "South American 1969","WGS 60","WGS 66","WGS-72", "WGS-84"}


    variable a = EquatorialRadius[ReferenceEllipsoid]
    variable eccSquared = SquareOfEccentricity[ReferenceEllipsoid];
    variable k0 = 0.9996;

    variable LongOrigin;
    variable eccPrimeSquared;
    variable N, T, C, AA, M;
   
//Make sure the longitude is between -180.00 .. 179.9
    variable LongTemp = (Long+180)-floor((Long+180)/360)*360-180; // -180.00 .. 179.9;

    variable LatRad = Lat*DEG2RAD;
    variable LongRad = LongTemp*DEG2RAD;
    variable LongOriginRad;
    variable    ZoneNumber;

    ZoneNumber = floor((LongTemp + 180)/6) + 1;
 
    if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 )
        ZoneNumber = 32;
    endif
  // Special zones for Svalbard
    if( Lat >= 72.0 && Lat < 84.0 )
      if  ( LongTemp >= 0.0  && LongTemp <  9.0 )
        ZoneNumber = 31;
      elseif( LongTemp >= 9.0  && LongTemp < 21.0 )
        ZoneNumber = 33;
      elseif(LongTemp >= 21.0 && LongTemp < 33.0 )
        ZoneNumber = 35;
      elseif(LongTemp >= 33.0 && LongTemp < 42.0 )
        ZoneNumber = 37;
      endif
     endif
    LongOrigin = (ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
    LongOriginRad = LongOrigin * DEG2RAD;

    //compute the UTM Zone from the latitude and longitude
   
    eccPrimeSquared = (eccSquared)/(1-eccSquared);

    N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
    T = tan(LatRad)*tan(LatRad);
    C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
    AA = cos(LatRad)*(LongRad-LongOriginRad);
   
    M = (1  - eccSquared/4      - 3*eccSquared*eccSquared/64    - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
    M -=  (3*eccSquared/8   + 3*eccSquared*eccSquared/32    + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
    M += (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
    M -=  (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad)
    M *= a
   
    UTMEasting = (k0*N*(AA+(1-T+C)*AA*AA*AA/6+ (5-18*T+T*T+72*C-58*eccPrimeSquared)*AA*AA*AA*AA*AA/120)+ 500000.0);

    UTMNorthing = (k0*(M+N*tan(LatRad)*(AA*AA/2+(5-T+9*C+4*C*C)*AA*AA*AA*AA/24+ (61-58*T+T*T+600*C-330*eccPrimeSquared)*AA*AA*AA*AA*AA*AA/720)));
    if(Lat < 0)
        UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
    endif

    gps.UTMEasting = UTMEasting
    gps.UTMNorthing = UTMNorthing
    gps.UTMZone = num2str(ZoneNumber)+ UTMLetterDesignator(Lat)
    print gps
End

Function/t UTMLetterDesignator( Lat)
variable Lat
//This routine determines the correct UTM letter designator for the given latitude
//returns "Z" if latitude is outside the UTM limits of 84N to 80S
    //Written by Chuck Gantz- chuck.gantz@globalstar.com
    string LetterDesignator;

    if((84 >= Lat) && (Lat >= 72))
        LetterDesignator = "X";
    elseif((72 > Lat) && (Lat >= 64))
         LetterDesignator = "W";
    elseif((64 > Lat) && (Lat >= 56))
        LetterDesignator = "V";
    elseif((56 > Lat) && (Lat >= 48))
        LetterDesignator = "U";
    elseif((48 > Lat) && (Lat >= 40))
     LetterDesignator = "T";
    elseif((40 > Lat) && (Lat >= 32))
     LetterDesignator = "S";
    elseif((32 > Lat) && (Lat >= 24))
     LetterDesignator = "R";
    elseif((24 > Lat) && (Lat >= 16))
     LetterDesignator = "Q";
    elseif((16 > Lat) && (Lat >= 8))
    LetterDesignator = "P";
    elseif(( 8 > Lat) && (Lat >= 0))
    LetterDesignator = "N";
    elseif(( 0 > Lat) && (Lat >= -8))
    LetterDesignator = "M";
    elseif((-8> Lat) && (Lat >= -16))
    LetterDesignator = "L";
    elseif((-16 > Lat) && (Lat >= -24))
     LetterDesignator = "K";
    elseif((-24 > Lat) && (Lat >= -32))
     LetterDesignator = "J";
    elseif((-32 > Lat) && (Lat >= -40))
    LetterDesignator = "H";
    elseif((-40 > Lat) && (Lat >= -48))
    LetterDesignator = "G";
    elseif((-48 > Lat) && (Lat >= -56))
     LetterDesignator = "F";
    elseif((-56 > Lat) && (Lat >= -64))
    LetterDesignator = "E";
    elseif((-64 > Lat) && (Lat >= -72))
     LetterDesignator = "D";
    elseif((-72 > Lat) && (Lat >= -80))
    LetterDesignator = "C";
    else
         LetterDesignator = "Z"; //This is here as an error flag to show that the Latitude is outside the UTM limits
    endif
    return LetterDesignator;
End
Hi,
how about using the more compact
Function/S UTMLetterDesignator(Lat)
variable Lat
    string LetterDesignator
    if( -80 <= Lat && Lat <= 84 )
        LetterDesignator = "CDEFGHJKLMNPQRSTUVWXX"[(Lat + 80) / 8]
    else
        LetterDesignator = "Z"; //This is here as an error flag to show that the Latitude is outside the UTM limits
    endif
    return LetterDesignator;
End

and why do you use Function/T? Is this an undocumented flag?

Wolfgang Harneit
Hi Wolfgang,
I just copied it from C code to IGOR, I wasn't looking to make it fancy, I just did as little as possible to make it work ;-). You are right, it should be Function/S.
Your UTM designator function is a lot more elegant.

Andy
I just came across this site. It's nice to know my code is still being used after all these years.

I agree Wolfgang's code is more elegant. The reason I used the long dragged out code was just to make it clear what I was doing. By writing out
if((84 >= Lat) && (Lat >= 72))
LetterDesignator = "X";
elseif((72 > Lat) && (Lat >= 64))
LetterDesignator = "W";
etc., it's obvious what the LetterDesignators mean. Since compactness and speed weren't important, the pedogical advantage seemed the way to go.

Chuck
The function parameters need fixing, too:

Function LLtoUTM( ReferenceEllipsoid,   Lat,  Long,)
    variable ReferenceEllipsoid, Lat, Long
    struct GPSpos gps


should be:

Function LLtoUTM( ReferenceEllipsoid,   Lat,  Long, gps)
    variable ReferenceEllipsoid, Lat, Long
    struct GPSpos &gps  // output: need the & here to return structure

--Jim Prouty
Software Engineer, WaveMetrics, Inc.

Forum

Support

Gallery

Igor Pro 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More