UtmConvert: C:/Documents and Settings/student/My Documents/Visual Studio 2008/projects/system_controller/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         //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 }