- Main Page
- Packages
- Classes
- Files

- File List
- File Members

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

Go to the documentation of this file.

    00001 ﻿using System;
    00002 using System.Collections.Generic;
    00003 using System.Linq;
    00004 using System.Text;
    00005
    00006 namespace UtmConvert {
    00007     class LatLongToUtm {
    00008         /*
    00009 Symbols
    00010
    00011     * lat = latitude of point
    00012     * long = longitude of point
    00013     * long0 = central meridian of zone
    00014     * 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.
    00015     * e = SQRT(1-b^2/a^2) = .08 approximately. This is the eccentricity of the earth's elliptical cross-section.
    00016     * 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.
    00017     * n = (a-b)/(a+b)
    00018     * 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.
    00019     * 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.
    00020     * p = (long-long0) in radians (This differs from the treatment in the Army reference)
    00021
    00022 Calculate the Meridional Arc
    00023
    00024 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.
    00025
    00026     * S = A'lat - B'sin(2lat) + C'sin(4lat) - D'sin(6lat) + E'sin(8lat), where lat is in radians and
    00027     * A' = a[1 - n + (5/4)(n^2 - n^3) + (81/64)(n^4 - n^5) ...]
    00028     * B' = (3 tan/2)[1 - n + (7/8)(n^2 - n^3) + (55/64)(n^4 - n^5) ...]
    00029     * C' = (15 tan^2/16)[1 - n + (3/4)(n^2 - n^3)  ...]
    00030     * D' = (35 tan^3/48)[1 - n + (11/16)(n^2 - n^3)  ...]
    00031     * E' = (315 tan^4/512)[1 - n  ...]
    00032
    00033 The USGS gives this form, which may be more appealing to some. (They use M where the Army uses S)
    00034
    00035     * M = a[(1 - e^2/4 - 3e^4/64 - 5e^6/256 ....)lat
    00036        - (3e^2/8 + 3e^4/32 + 45e^6/1024...)sin(2lat)
    00037       + (15e^4/256 + 45e^6/1024 + ....)sin(4lat)
    00038        - (35e^6/3072 + ....) sin(6lat) + ....)] where lat is in radians
    00039
    00040 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.
    00041 Converting Latitude and Longitude to UTM
    00042
    00043 All angles are in radians.
    00044
    00045 y = northing = K1 + K2p^2 + K3p^4, where
    00046
    00047     * K1 = Sk0,
    00048     * K2 = k0 nu sin(lat)cos(lat)/2 = k0 nu sin(2 lat)/4
    00049     * K3 = [k0 nu sin(lat)cos^3(lat)/24][(5 - tan^2(lat) + 9e'^2 cos^2(lat) + 4e'^4 cos4^(lat)]
    00050
    00051 x = easting = K4p + K5p^3, where
    00052
    00053     * K4 = k0 nu cos(lat)
    00054     * K5 = (k0 nu cos^3(lat)/6)[1 - tan^2(lat) + e'^2 cos^2(lat)]
    00055
    00056 Easting x is relative to the central meridian. For conventional UTM easting add 500,000 meters to x.
    00057
    00058 What the Formulas Mean
    00059
    00060 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.
    00061 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.
    00062 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,
    00063 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.
    00064 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.
    00065
    00066 References
    00067
    00068 http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm Steven Dutch, Natural and Applied Sciences, University of Wisconsin - Green Bay
    00069 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.
    00070 Army, Department of, 1973; Universal Transverse Mercator Grid, U. S. Army Technical Manual TM 5-241-8, 64 p.
    00071 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.
    00072          */
    00073         Datum _datum;
    00074         public Datum Datum {
    00075             get { return _datum; }
    00076             set { _datum = value; }
    00077         }
    00078
    00079         double long0 = -123; //central meridian of zone (radians) // maybe should be 123
    00080         const double k0 = 0.9996; //scale along long0
    00081         double e = 0.08; //approximately. This is the eccentricity of the earth's elliptical cross-section.
    00082         double eTicSq = 0.007; //The quantity e' only occurs in even powers so it need only be calculated as e'^2.
    00083         double n = 0;
    00084         double rho = 0; //This is the radius of curvature of the earth in the meridian plane.
    00085         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.
    00086         double p = 0; //in radians (This differs from the treatment in the Army reference)
    00087         //=a*(1-n+(5*n*n/4)*(1-n)+(81*n^4/64)*(1-n))
    00088         double A0 = 0;
    00089         //=(3*a*n/2)*(1-n-(7*n*n/8)*(1-n)+55*n^4/64)
    00090         double B0 = 0;
    00091         //=(15*a*n*n/16)*(1-n+(3*n*n/4)*(1-n))
    00092         double C0 = 0;
    00093         //=(35*a*n^3/48)*(1-n+11*n*n/16)
    00094         double D0 = 0;
    00095         //=(315*a*n^4/51)*(1-n)
    00096         double E0 = 0;
    00097         //NOAA Meridional Arc
    00098         //=A0*lat-B0*SIN(2*lat)+C0*SIN(4*lat)-D0*SIN(6*lat)+E0*SIN(8*lat)
    00099         double S = 0;
    00100         //=S*k0
    00101         double K1 = 0;
    00102         //=nu*SIN(lat)*COS(lat)*k0/2
    00103         double K2 = 0;
    00104         //=((nu*SIN(lat)*COS(lat)^3)/24)*(5-TAN(lat)^2+9*e1sq*COS(lat)^2+4*e1sq^2*COS(lat)^4)*k0
    00105         double K3 = 0;
    00106         //=nu*COS(lat)*k0
    00107         double K4 = 0;
    00108         //=(COS(lat))^3*(nu/6)*(1-TAN(lat)^2+e1sq*COS(lat)^2)*k0
    00109         double K5 = 0;
    00110         //army Meridional Arc
    00111         //double M = 0;
    00112
    00113         double _easting = 491887;
    00114         public double Easting {
    00115             get { return _easting; }
    00116         }
    00117         public double x {
    00118             get { return _easting; }
    00119         }
    00120         double _northing = 4876938;
    00121         public double Northing {
    00122             get { return _northing; }
    00123         }
    00124         public double y {
    00125             get { return _northing; }
    00126         }
    00127
    00128         string _zone = "";
    00129         public string Zone {
    00130             get { return _zone; }
    00131         }
    00132
    00133         int _centralMeridian = 0;
    00134         public int CentralMeridian {
    00135             get { return _centralMeridian; }
    00136         }
    00137
    00138
    00139         public LatLongToUtm() {
    00140             _datum = new Datum();
    00141
    00142         }
    00143
    00144
    00145
    00146         public void getUtm(double lat, double lon) {
    00147             calcCM(lon);
    00148             e = Math.Sqrt(1.0 - (Math.Pow(_datum.b / _datum.a, 2)));
    00149             eTicSq = Math.Pow(e, 2) / (1.0 - Math.Pow(e, 2));
    00150             n = (_datum.a - _datum.b) / (_datum.a + _datum.b);
    00151             rho = _datum.a * (1.0 - Math.Pow(e, 2)) / Math.Pow(1.0
    00152                 - Math.Pow(e * Math.Sin(lat), 2), 3.0 / 2.0);
    00153             nu = _datum.a / Math.Pow(1.0 - Math.Pow(e * Math.Sin(lat), 2), 1.0 / 2.0);
    00154             p = lon - long0;
    00155             A0 = _datum.a * (1.0 - n + (5.0 * Math.Pow(n,2) / 4.0) * (1.0 - n)
    00156                 + (81.0 * Math.Pow(n, 4) / 64.0) * (1.0 - n));
    00157             B0 = (3.0 * _datum.a * n / 2.0) * (1.0 - n - (7.0 * Math.Pow(n,2) / 8.0)
    00158                 * (1.0 - n) + 55.0 * Math.Pow(n, 4) / 64.0);
    00159             C0 = (15.0 * _datum.a * Math.Pow(n,2) / 16.0)
    00160                 * (1.0 - n + (3.0 * Math.Pow(n,2) / 4.0) * (1.0 - n));
    00161             D0 = (35.0 * _datum.a * Math.Pow(n, 3) / 48.0) * (1.0 - n + 11.0 * Math.Pow(n, 2) / 16.0);
    00162             E0 = (315.0 * _datum.a * Math.Pow(n, 4) / 51.0) * (1.0 - n);
    00163
    00164             S = A0 * lat - B0 * Math.Sin(2.0 * lat) + C0 * Math.Sin(4.0 * lat)
    00165                 + E0 * Math.Sin(8.0 * lat);
    00166
    00168             //M = _datum.a * ((1.0 - Math.Pow(e, 2) / 4.0 - 3.0 * Math.Pow(e, 4)
    00169             //    / 64.0 - 5.0 * Math.Pow(e, 6) / 256.0) * lat
    00170             //    - (3.0 * Math.Pow(e, 2) / 8.0 + 3.0 * Math.Pow(e, 4)
    00171             //    / 32.0 + 45.0 * Math.Pow(e, 6) / 1024) * Math.Sin(2.0 * lat)
    00172             //    + (15.0 * Math.Pow(e, 4) / 256.0 + 45.0 * Math.Pow(e, 6) / 1024.0)
    00173             //    * Math.Sin(4 * lat)
    00174             //    - (35.0 * Math.Pow(e, 6) / 3072.0) * Math.Sin(6.0 * lat));
    00175
    00176             K1 = S * k0;
    00177             K2 = k0 * nu * Math.Sin(lat) * Math.Cos(lat) / 2.0;
    00178             K3 = (k0 * nu * Math.Sin(lat) * Math.Pow(Math.Cos(lat), 3) / 24.0)
    00179                 * (5.0 - Math.Pow(Math.Tan(lat), 2) + 9.0 * eTicSq
    00180                 * Math.Pow(Math.Cos(lat), 2)
    00181                 + 4.0 * Math.Pow(eTicSq, 2) * Math.Pow(Math.Cos(lat), 4));
    00182             _northing = K1 + K2 * Math.Pow(p, 2) + K3 * Math.Pow(p, 4);
    00183             K4 = nu * Math.Cos(lat) * k0;
    00184             K5 = Math.Pow(Math.Cos(lat), 3) * (nu / 6.0)
    00185                 * (1 - Math.Pow(Math.Tan(lat), 2)
    00186                 + eTicSq * Math.Pow(Math.Cos(lat), 2)) * k0;
    00187             _easting = 500000 + (K4 * p + K5 * Math.Pow(p, 3));
    00188
    00189             char c = 'N';
    00190             if (lat > 0)
    00191                 c = (char)(((int)((int)(lat * (180.0 / Math.PI))) / 8) + 78);
    00192             else
    00193                 c = (char)(((int)((int)(lat * (180.0 / Math.PI))) / 8) + 77);
    00194             if (c >= 'O') c++;
    00195
    00196             if (c == 'I') c--;
    00197             _zone = ((int)(((lon * (180.0 / Math.PI)) + 180) / 6) + 1).ToString() + c;
    00198         }
    00199
    00200         private void calcCM(double lon) {
    00201             if (lon < 0)
    00202                 _centralMeridian = (int)(((int)(Math.Abs(lon) * (180.0 / Math.PI))) / 6) * -6 - 3;
    00203             else
    00204                 _centralMeridian = (int)(((int)(Math.Abs(lon) * (180.0 / Math.PI))) / 6) * 6 + 3;
    00205             long0 = (double)_centralMeridian * (Math.PI / 180.0);
    00206         }
    00207
    00208     }
    00209 }