C:/Documents and Settings/student/My Documents/Visual Studio 2008/projects/UtmConvert/UtmConvert/LatLongToUtm.cs

Go to the documentation of this file.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace UtmConvert {
    class LatLongToUtm {
        /*
Symbols

    * lat = latitude of point
    * long = longitude of point
    * long0 = central meridian of zone
    * k0  = scale along long0 = 0.9996. Even though it's a constant, we retain it as a separate symbol to keep the numerical coefficients simpler, also to allow for systems that might use a different Mercator projection.
    * e = SQRT(1-b^2/a^2) = .08 approximately. This is the eccentricity of the earth's elliptical cross-section.
    * e'^2 = (ea/b)^2 = e^2/(1-e^2) = .007 approximately. The quantity e' only occurs in even powers so it need only be calculated as e'2.
    * n = (a-b)/(a+b)
    * rho = a(1-e^2)/(1-e^2 sin^2(lat))^(3/2). This is the radius of curvature of the earth in the meridian plane.
    * nu = a/(1-e^2 sin^2(lat))^(1/2). This is the radius of curvature of the earth perpendicular to the meridian plane. It is also the distance from the point in question to the polar axis, measured perpendicular to the earth's surface.
    * p = (long-long0) in radians (This differs from the treatment in the Army reference)

Calculate the Meridional Arc

S is the meridional arc through the point in question (the distance along the earth's surface from the equator). All angles are in radians.

    * S = A'lat - B'sin(2lat) + C'sin(4lat) - D'sin(6lat) + E'sin(8lat), where lat is in radians and 
    * A' = a[1 - n + (5/4)(n^2 - n^3) + (81/64)(n^4 - n^5) ...]
    * B' = (3 tan/2)[1 - n + (7/8)(n^2 - n^3) + (55/64)(n^4 - n^5) ...]
    * C' = (15 tan^2/16)[1 - n + (3/4)(n^2 - n^3)  ...]
    * D' = (35 tan^3/48)[1 - n + (11/16)(n^2 - n^3)  ...]
    * E' = (315 tan^4/512)[1 - n  ...]

The USGS gives this form, which may be more appealing to some. (They use M where the Army uses S)

    * M = a[(1 - e^2/4 - 3e^4/64 - 5e^6/256 ....)lat
       - (3e^2/8 + 3e^4/32 + 45e^6/1024...)sin(2lat) 
      + (15e^4/256 + 45e^6/1024 + ....)sin(4lat)
       - (35e^6/3072 + ....) sin(6lat) + ....)] where lat is in radians

This is the hard part. Calculating the arc length of an ellipse involves functions called elliptic integrals, which don't reduce to neat closed formulas. So they have to be represented as series.
Converting Latitude and Longitude to UTM

All angles are in radians.

y = northing = K1 + K2p^2 + K3p^4, where

    * K1 = Sk0,
    * K2 = k0 nu sin(lat)cos(lat)/2 = k0 nu sin(2 lat)/4
    * K3 = [k0 nu sin(lat)cos^3(lat)/24][(5 - tan^2(lat) + 9e'^2 cos^2(lat) + 4e'^4 cos4^(lat)]

x = easting = K4p + K5p^3, where

    * K4 = k0 nu cos(lat)
    * K5 = (k0 nu cos^3(lat)/6)[1 - tan^2(lat) + e'^2 cos^2(lat)]

Easting x is relative to the central meridian. For conventional UTM easting add 500,000 meters to x.

What the Formulas Mean

The hard part, allowing for the oblateness of the Earth, is taken care of in calculating S (or M). So K1 is simply the arc length along the central meridian of the zone corrected by the scale factor. Remember, the scale is a hair less than 1 in the middle of the zone, and a hair more on the outside.
All the higher K terms involve nu, the local radius of curvature (roughly equal to the radius of the earth or roughly 6,400,000 m), trig functions, and powers of e'2 ( = .007 ). So basically they are never much larger than nu. Actually the maximum value of K2 is about nu/4 (1,600,000), K3 is about nu/24 (267,000) and K5 is about nu/6 (1,070,000). Expanding the expressions will show that the tangent terms don't affect anything.
If we were just to stop with the K2 term in the northing, we'd have a quadratic in p. In other words, we'd approximate the parallel of latitude as a parabola. The real curve is more complex. It will be more like a hyperbola equatorward of about 45 degrees and an ellipse poleward, at least within the narrow confines of a UTM zone. (At any given latitude we're cutting the cone of latitude vectors with an inclined plane,
so the resulting intersection will be a conic section. Since the projection cylinder has a curvature, the exact curve is not a conic but the difference across a six-degree UTM zone is pretty small.)  Hence the need for higher order terms. Now p will never be more than +/-3 degrees = .05 radians, so p2 is always less than .0025 (1/400) and p4 is always less than .00000625 (1/160000). Using a spreadsheet, it's easy to see how the individual terms vary with latitude. K2p2 never exceeds 4400 and K3p4 is at most a bit over 3. That is, the curvature of a parallel of latitude across a UTM zone is at most a little less than 4.5 km and the maximum departure from a parabola is at most a few meters.
K4 is what we'd calculate for easting in a simple-minded way, just by calculating arc distance along the parallel of latutude. But, as we get farther from the central meridian, the meridians curve inward, so our actual easting will be less than K4. That's what K5 does. Since p is never more than +/-3 degrees = .05 radians, p3 is always less than .000125 (1/8000). The maximum value of K5p3 is about 150 meters.  

References

http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm Steven Dutch, Natural and Applied Sciences, University of Wisconsin - Green Bay
Snyder, J. P., 1987; Map Projections - A Working Manual. U.S. Geological Survey Professional Paper 1395, 383 p. If you are at all serious about maps you need this book.
Army, Department of, 1973; Universal Transverse Mercator Grid, U. S. Army Technical Manual TM 5-241-8, 64 p.
NIMA Technical Report 8350.2, "Department of Defense World Geodetic System 1984, Its Definition and Relationships with Local Geodetic Systems," Second Edition, 1 September 1991 and its supplements. The report is available from the NIMA Combat Support Center and its stock number is DMATR83502WGS84. Non-DoD requesters may obtain the report as a public sale item from the U.S. Geological Survey, Box 25286, Denver Federal Center, Denver, Colorado 80225 or by phone at 1-800-USA-MAPS.   
         */
        Datum _datum;
        public Datum Datum {
            get { return _datum; }
            set { _datum = value; }
        }

        double long0 = -123; //central meridian of zone (radians) // maybe should be 123
        const double k0 = 0.9996; //scale along long0 
        double e = 0.08; //approximately. This is the eccentricity of the earth's elliptical cross-section.
        double eTicSq = 0.007; //The quantity e' only occurs in even powers so it need only be calculated as e'^2.
        double n = 0;
        double rho = 0; //This is the radius of curvature of the earth in the meridian plane.
        double nu = 0; //This is the radius of curvature of the earth perpendicular to the meridian plane. It is also the distance from the point in question to the polar axis, measured perpendicular to the earth's surface.
        double p = 0; //in radians (This differs from the treatment in the Army reference)
        //=a*(1-n+(5*n*n/4)*(1-n)+(81*n^4/64)*(1-n))
        double A0 = 0;
        //=(3*a*n/2)*(1-n-(7*n*n/8)*(1-n)+55*n^4/64)
        double B0 = 0;
        //=(15*a*n*n/16)*(1-n+(3*n*n/4)*(1-n))
        double C0 = 0;
        //=(35*a*n^3/48)*(1-n+11*n*n/16)
        double D0 = 0;
        //=(315*a*n^4/51)*(1-n)
        double E0 = 0;
        //NOAA Meridional Arc
        //=A0*lat-B0*SIN(2*lat)+C0*SIN(4*lat)-D0*SIN(6*lat)+E0*SIN(8*lat)
        double S = 0;
        //=S*k0
        double K1 = 0;
        //=nu*SIN(lat)*COS(lat)*k0/2
        double K2 = 0;
        //=((nu*SIN(lat)*COS(lat)^3)/24)*(5-TAN(lat)^2+9*e1sq*COS(lat)^2+4*e1sq^2*COS(lat)^4)*k0
        double K3 = 0;
        //=nu*COS(lat)*k0
        double K4 = 0;
        //=(COS(lat))^3*(nu/6)*(1-TAN(lat)^2+e1sq*COS(lat)^2)*k0
        double K5 = 0;
        //army Meridional Arc
        //double M = 0;

        double _easting = 491887;
        public double Easting {
            get { return _easting; }
        }
        public double x {
            get { return _easting; }
        }
        double _northing = 4876938;
        public double Northing {
            get { return _northing; }
        }
        public double y {
            get { return _northing; }
        }

        string _zone = "";
        public string Zone {
            get { return _zone; }
        }

        int _centralMeridian = 0;
        public int CentralMeridian { 
            get { return _centralMeridian; }
        }


        public LatLongToUtm() {
            _datum = new Datum();
            
        }



        public void getUtm(double lat, double lon) {
            calcCM(lon);
            e = Math.Sqrt(1.0 - (Math.Pow(_datum.b / _datum.a, 2)));
            eTicSq = Math.Pow(e, 2) / (1.0 - Math.Pow(e, 2));
            n = (_datum.a - _datum.b) / (_datum.a + _datum.b);
            rho = _datum.a * (1.0 - Math.Pow(e, 2)) / Math.Pow(1.0
                - Math.Pow(e * Math.Sin(lat), 2), 3.0 / 2.0);
            nu = _datum.a / Math.Pow(1.0 - Math.Pow(e * Math.Sin(lat), 2), 1.0 / 2.0);
            p = lon - long0;
            A0 = _datum.a * (1.0 - n + (5.0 * Math.Pow(n,2) / 4.0) * (1.0 - n) 
                + (81.0 * Math.Pow(n, 4) / 64.0) * (1.0 - n));
            B0 = (3.0 * _datum.a * n / 2.0) * (1.0 - n - (7.0 * Math.Pow(n,2) / 8.0) 
                * (1.0 - n) + 55.0 * Math.Pow(n, 4) / 64.0);
            C0 = (15.0 * _datum.a * Math.Pow(n,2) / 16.0) 
                * (1.0 - n + (3.0 * Math.Pow(n,2) / 4.0) * (1.0 - n));
            D0 = (35.0 * _datum.a * Math.Pow(n, 3) / 48.0) * (1.0 - n + 11.0 * Math.Pow(n, 2) / 16.0);
            E0 = (315.0 * _datum.a * Math.Pow(n, 4) / 51.0) * (1.0 - n);
            
            S = A0 * lat - B0 * Math.Sin(2.0 * lat) + C0 * Math.Sin(4.0 * lat)
                + E0 * Math.Sin(8.0 * lat);

            //M = _datum.a * ((1.0 - Math.Pow(e, 2) / 4.0 - 3.0 * Math.Pow(e, 4)
            //    / 64.0 - 5.0 * Math.Pow(e, 6) / 256.0) * lat
            //    - (3.0 * Math.Pow(e, 2) / 8.0 + 3.0 * Math.Pow(e, 4)
            //    / 32.0 + 45.0 * Math.Pow(e, 6) / 1024) * Math.Sin(2.0 * lat)
            //    + (15.0 * Math.Pow(e, 4) / 256.0 + 45.0 * Math.Pow(e, 6) / 1024.0)
            //    * Math.Sin(4 * lat)
            //    - (35.0 * Math.Pow(e, 6) / 3072.0) * Math.Sin(6.0 * lat));
            
            K1 = S * k0;
            K2 = k0 * nu * Math.Sin(lat) * Math.Cos(lat) / 2.0;
            K3 = (k0 * nu * Math.Sin(lat) * Math.Pow(Math.Cos(lat), 3) / 24.0)
                * (5.0 - Math.Pow(Math.Tan(lat), 2) + 9.0 * eTicSq
                * Math.Pow(Math.Cos(lat), 2)
                + 4.0 * Math.Pow(eTicSq, 2) * Math.Pow(Math.Cos(lat), 4));
            _northing = K1 + K2 * Math.Pow(p, 2) + K3 * Math.Pow(p, 4);
            K4 = nu * Math.Cos(lat) * k0;
            K5 = Math.Pow(Math.Cos(lat), 3) * (nu / 6.0)
                * (1 - Math.Pow(Math.Tan(lat), 2)
                + eTicSq * Math.Pow(Math.Cos(lat), 2)) * k0;            
            _easting = 500000 + (K4 * p + K5 * Math.Pow(p, 3));

            char c = 'N';
            if (lat > 0)
                c = (char)(((int)((int)(lat * (180.0 / Math.PI))) / 8) + 78);
            else
                c = (char)(((int)((int)(lat * (180.0 / Math.PI))) / 8) + 77);
            if (c >= 'O') c++;

            if (c == 'I') c--;
            _zone = ((int)(((lon * (180.0 / Math.PI)) + 180) / 6) + 1).ToString() + c;
        }

        private void calcCM(double lon) {
            if (lon < 0)
                _centralMeridian = (int)(((int)(Math.Abs(lon) * (180.0 / Math.PI))) / 6) * -6 - 3;
            else
                _centralMeridian = (int)(((int)(Math.Abs(lon) * (180.0 / Math.PI))) / 6) * 6 + 3;
            long0 = (double)_centralMeridian * (Math.PI / 180.0);
        }

    }
}

---