UtmConvert: C:/Documents and Settings/student/My Documents/Visual Studio 2008/projects/UtmConvert/UtmConvert/ConvertLatLonUtm.cs Source File

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