- Main Page
- Packages
- Classes
- Files

- File List
- File Members

# C:/Documents and Settings/student/My Documents/Visual Studio 2008/projects/system_controller/UtmConvert/UtmConvert/ConvertLatLonUtm.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         //references:
    00008         //Converting UTM to Latitude and Longitude (Or Vice Versa)
    00009         //Steven Dutch, Natural and Applied Sciences, University of Wisconsin - Green Bay
    00010         //http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM
    00011
    00012     public class ConvertLatLonUtm {
    00013         Datum _datum;
    00014         public Datum Datum {
    00015             get { return _datum; }
    00016             set { _datum = value; }
    00017         }
    00018
    00019         double _easting = 491887;
    00020         public double Easting {
    00021             get { return _easting; }
    00022         }
    00023
    00024         public double x {
    00025             get { return _easting; }
    00026         }
    00027
    00028         double _northing = 4876938;
    00029         public double Northing {
    00030             get { return _northing; }
    00031         }
    00032
    00033         public double y {
    00034             get { return _northing; }
    00035         }
    00036
    00037         string _zone = "";
    00038         public string Zone {
    00039             get { return _zone; }
    00040         }
    00041
    00042         int _centralMeridian = 0;
    00043         public int CentralMeridian {
    00044             get { return _centralMeridian; }
    00045         }
    00046
    00047         double _lat;
    00048         public double Latitude {
    00049             get { return _lat; }
    00050         }
    00051
    00052         double _lon;
    00053         public double Longitude {
    00054             get { return _lon; }
    00055         }
    00056
    00057         public ConvertLatLonUtm() {
    00058             _datum = new Datum();
    00059
    00060         }
    00061
    00062         const double k0 = 0.9996; //scale along long0
    00063         double e = 0.08; //approximately. This is the eccentricity of the earth's elliptical cross-section.
    00064         double eTic = 0.007; //The quantity e' only occurs in even powers so it need only be calculated as e'^2.
    00065         public void convertLatLonToUtm(double lat, double lon) {
    00066             calcCM(lon);
    00067             e = Math.Sqrt(1.0 - (Math.Pow(_datum.b / _datum.a, 2)));
    00068             eTic = Math.Pow(e, 2) / (1.0 - Math.Pow(e, 2));
    00069             double n = (_datum.a - _datum.b) / (_datum.a + _datum.b);
    00070             //This is the radius of curvature of the earth in the meridian plane.
    00071             double rho = _datum.a * (1.0 - Math.Pow(e, 2)) / Math.Pow(1.0
    00072                 - Math.Pow(e * Math.Sin(lat), 2), 3.0 / 2.0);
    00073             //This is the radius of curvature of the earth perpendicular to the meridian plane. It is also the
    00074             //distance from the point in question to the polar axis, measured perpendicular to the earth's surface.
    00075             double nu = _datum.a / Math.Pow(1.0 - Math.Pow(e * Math.Sin(lat), 2), 1.0 / 2.0);
    00076             //in radians (This differs from the treatment in the Army reference)
    00077             double p = lon - long0;
    00078             //=a*(1-n+(5*n*n/4)*(1-n)+(81*n^4/64)*(1-n))
    00079             double A0 = _datum.a * (1.0 - n + (5.0 * Math.Pow(n,2) / 4.0) * (1.0 - n)
    00080                 + (81.0 * Math.Pow(n, 4) / 64.0) * (1.0 - n));
    00081             //=(3*a*n/2)*(1-n-(7*n*n/8)*(1-n)+55*n^4/64)
    00082             double B0 = (3.0 * _datum.a * n / 2.0) * (1.0 - n - (7.0 * Math.Pow(n,2) / 8.0)
    00083                 * (1.0 - n) + 55.0 * Math.Pow(n, 4) / 64.0);
    00084             //=(15*a*n*n/16)*(1-n+(3*n*n/4)*(1-n))
    00085             double C0 = (15.0 * _datum.a * Math.Pow(n,2) / 16.0)
    00086                 * (1.0 - n + (3.0 * Math.Pow(n,2) / 4.0) * (1.0 - n));
    00087             //=(35*a*n^3/48)*(1-n+11*n*n/16)
    00088             double D0 = (35.0 * _datum.a * Math.Pow(n, 3) / 48.0) * (1.0 - n + 11.0 * Math.Pow(n, 2) / 16.0);
    00089             //=(315*a*n^4/51)*(1-n)
    00090             double E0 = (315.0 * _datum.a * Math.Pow(n, 4) / 51.0) * (1.0 - n);
    00091             //NOAA Meridional Arc
    00092             //=A0*lat-B0*SIN(2*lat)+C0*SIN(4*lat)-D0*SIN(6*lat)+E0*SIN(8*lat)
    00093             double S = A0 * lat - B0 * Math.Sin(2.0 * lat) + C0 * Math.Sin(4.0 * lat)
    00094                 + E0 * Math.Sin(8.0 * lat);
    00095
    00097             //double M = _datum.a * ((1.0 - Math.Pow(e, 2) / 4.0 - 3.0 * Math.Pow(e, 4)
    00098             //    / 64.0 - 5.0 * Math.Pow(e, 6) / 256.0) * lat
    00099             //    - (3.0 * Math.Pow(e, 2) / 8.0 + 3.0 * Math.Pow(e, 4)
    00100             //    / 32.0 + 45.0 * Math.Pow(e, 6) / 1024) * Math.Sin(2.0 * lat)
    00101             //    + (15.0 * Math.Pow(e, 4) / 256.0 + 45.0 * Math.Pow(e, 6) / 1024.0)
    00102             //    * Math.Sin(4 * lat)
    00103             //    - (35.0 * Math.Pow(e, 6) / 3072.0) * Math.Sin(6.0 * lat));
    00104
    00105             //=S*k0
    00106             double K1 = S * k0;
    00107             //=nu*SIN(lat)*COS(lat)*k0/2
    00108             double K2 = k0 * nu * Math.Sin(lat) * Math.Cos(lat) / 2.0;
    00109             //=((nu*SIN(lat)*COS(lat)^3)/24)*(5-TAN(lat)^2+9*e1sq*COS(lat)^2+4*e1sq^2*COS(lat)^4)*k0
    00110             double K3 = (k0 * nu * Math.Sin(lat) * Math.Pow(Math.Cos(lat), 3) / 24.0)
    00111                 * (5.0 - Math.Pow(Math.Tan(lat), 2) + 9.0 * eTic
    00112                 * Math.Pow(Math.Cos(lat), 2)
    00113                 + 4.0 * Math.Pow(eTic, 2) * Math.Pow(Math.Cos(lat), 4));
    00114             _northing = K1 + K2 * Math.Pow(p, 2) + K3 * Math.Pow(p, 4);
    00115             //=nu*COS(lat)*k0
    00116             double K4 = nu * Math.Cos(lat) * k0;
    00117             //=(COS(lat))^3*(nu/6)*(1-TAN(lat)^2+e1sq*COS(lat)^2)*k0
    00118             double K5 = Math.Pow(Math.Cos(lat), 3) * (nu / 6.0)
    00119                 * (1 - Math.Pow(Math.Tan(lat), 2)
    00120                 + eTic * Math.Pow(Math.Cos(lat), 2)) * k0;
    00121             _easting = 500000.0 + (K4 * p + K5 * Math.Pow(p, 3));
    00122             char c = 'N';
    00123             //implicit casting is suppost to be going on here.
    00124             if (lat > 0)
    00125                 c = (char)(((int)((int)(lat * (180.0 / Math.PI))) / 8) + 78);
    00126             else
    00127                 c = (char)(((int)((int)(lat * (180.0 / Math.PI))) / 8) + 77);
    00128             if (c >= 'O') c++; //the letters I and O are not used because they can
    00129             if (c == 'I') c--; //be confused with numbers.
    00130             _zone = ((int)(((lon * (180.0 / Math.PI)) + 180) / 6) + 1).ToString() + c;
    00131         }
    00132
    00133         double long0 = -123; //central meridian of zone (radians) // maybe should be 123
    00134         private void calcCM(double lon) {
    00135             if (lon < 0)
    00136                 _centralMeridian = (int)(((int)(Math.Abs(lon)
    00137                     * (180.0 / Math.PI))) / 6) * -6 - 3;
    00138             else
    00139                 _centralMeridian = (int)(((int)(Math.Abs(lon)
    00140                     * (180.0 / Math.PI))) / 6) * 6 + 3;
    00141             long0 = (double)_centralMeridian * (Math.PI / 180.0);
    00142         }
    00143
    00144
    00145         public void convertUtmToLatLon(double x, double y, string zone) {
    00146             e = Math.Sqrt(1.0 - (Math.Pow(_datum.b / _datum.a, 2)));
    00147             eTic = Math.Pow(e, 2) / (1.0 - Math.Pow(e, 2));
    00148             if(int.Parse(zone.Substring(0, zone.Length - 1)) > 30)
    00149                 calcCM(ConvertDegRad.getRadians((((int.Parse(zone
    00150                     .Substring(0, zone.Length - 1))- 30) * 6) - 3).ToString(), 'E'));
    00151             else
    00152                 calcCM(ConvertDegRad.getRadians(((180 - (int.Parse(zone.Substring
    00153                     (0, zone.Length - 1)) * 6)) + 3).ToString(), 'W'));
    00154             x = x - 500000.0; //subtract 500,000 from conventional UTM coordinate
    00155             double M = y / k0;
    00156             double mu = M / (_datum.a * (1 - Math.Pow(e, 2) / 4.0 - 3.0
    00157                 * Math.Pow(e, 4) / 64.0 - 5.0 * Math.Pow(e, 6) / 256.0));
    00158             double e1 = (1.0 - Math.Pow(1.0 - Math.Pow(e, 2), 0.5))
    00159                 / (1.0 + Math.Pow(1.0 - Math.Pow(e, 2), 0.5));
    00160             //= (3 * e1 / 2 - 27 e1^3 / 32 ..)
    00161             double J1 = 3.0 * e1 / 2.0 - 27.0 * Math.Pow(e1, 3) / 32.0;
    00162             //= (21 * e1^2 / 16 - 55 * e1^4 / 32 ..)
    00163             double J2 = 21.0 * Math.Pow(e1, 2) / 16.0 - 55.0 * Math.Pow(e1, 4) / 32.0;
    00164             //(151 * e1^3 / 96 ..)
    00165             double J3 = 151.0 * Math.Pow(e1, 3) / 96.0;
    00166             //(1097 * e1^4 / 512 ..)
    00167             double J4 = 1097.0 * Math.Pow(e1, 4) / 512.0;
    00168             //footprint latitude =mu+J1*SIN(2*mu)+J2*SIN(4*mu)+J3*SIN(6*mu)+J4*SIN(8*mu)
    00169             double fp = mu + J1 * Math.Sin(2.0 * mu) + J2 * Math.Sin(4.0 * mu)
    00170                 + J3 * Math.Sin(6.0 * mu) + J4 * Math.Sin(8.0 * mu);
    00171             //e'^2 * cos^2(fp)
    00172             double C1 = Math.Pow(eTic, 2) * Math.Pow(Math.Cos(fp), 2);
    00173             //tan^2(fp)
    00174             double T1 = Math.Pow(Math.Tan(fp), 2);
    00175             //a(1 - e^2) / (1 - e^2 * sin^2(fp))^(3 / 2)
    00176             double R1 = _datum.a * (1.0 - Math.Pow(e, 2))
    00177                 / Math.Pow(1.0 - Math.Pow(e, 2) * Math.Pow(Math.Sin(fp), 2), 3.0 / 2.0);
    00178             //a / (1 - e^2 * sin^2(fp))^(1/2)
    00179             double N1 = _datum.a / Math.Pow(1.0 - Math.Pow(e, 2) * Math.Pow(Math.Sin(fp), 2), 1.0 / 2.0);
    00180             //x / (N1 * k0)
    00181             double D = x / (N1 * k0);
    00182             //N1 * tan(fp) / R1
    00183             double Q1 = N1 * Math.Tan(fp) / R1;
    00184             //(D^2 / 2)
    00185             double Q2 = Math.Pow(D, 2) / 2.0;
    00186             //(5 + 3 * T1 + 10 * C1 - 4 * C1^2 - 9 * e'^2) * D^4 / 24
    00187             double Q3 = (5.0 + 3.0 * T1 + 10.0 * C1 - 4.0 * Math.Pow(C1,2)
    00188                 - 9.0 * Math.Pow(eTic,2)) * Math.Pow(D,4) / 24.0;
    00189             //(61 + 90 * T1 + 298 * C1 + 45 * T1^2 - 3 * C1^2 - 252 * e'^2) * D^6 / 720
    00190             double Q4 = (61.0 + 90.0 * T1 + 289.0 * C1 + 45.0 * Math.Pow(T1, 2)
    00191                 - 3.0 * Math.Pow(C1, 2) - 252.0 * Math.Pow(eTic, 2)) * Math.Pow(D, 6) / 720.0;
    00192             _lat = fp - Q1 * (Q2 - Q3 + Q4);
    00193             double Q5 = D;
    00194             //(1 + 2 * T1 + C1)D^3 / 6
    00195             double Q6 = (1.0 + 2.0 * T1 + C1) * Math.Pow(D,3) / 6.0;
    00196             //(5 - 2 * C1 + 28 * T1 - 3 * C1^2 + 8 * e'^2 + 24 * T1^2) * D^5 / 120
    00197             double Q7 = (5.0 - 2.0 * C1 + 28.0 * T1 - 3.0 * Math.Pow(C1,2) + 8.0 * Math.Pow(eTic,2)
    00198                 + 24.0 * Math.Pow(T1,2)) * Math.Pow(D,5) / 120.0;
    00199             //long0 + (Q5 - Q6 + Q7) / cos(fp)
    00200             _lon = long0 + (Q5 - Q6 + Q7) / Math.Cos(fp);
    00201         }
    00202
    00203     }
    00204 }