GPS 纠偏，点距

 1     internal class EvilTransform
 2     {
 3         private const double pi = 3.14159265358979324;
 4 
 5         //
 6         // Krasovsky 1940
 7         //
 8         // a = 6378245.0, 1/f = 298.3
 9         // b = a * (1 - f)
10         // ee = (a^2 - b^2) / a^2;
11         private const double a = 6378245.0;
12         private const double ee = 0.00669342162296594323;
13 
14         //
15         // World Geodetic System ==> Mars Geodetic System
16         public static void Transform(double wgLat, double wgLon, out double mgLat, out double mgLon)
17         {
18             if (OutOfChina(wgLat, wgLon))
19             {
20                 mgLat = wgLat;
21                 mgLon = wgLon;
22                 return;
23             }
24             double dLat = TransformLat(wgLon - 105.0, wgLat - 35.0);
25             double dLon = TransformLon(wgLon - 105.0, wgLat - 35.0);
26             double radLat = wgLat/180.0*pi;
27             double magic = Math.Sin(radLat);
28             magic = 1 - ee*magic*magic;
29             double sqrtMagic = Math.Sqrt(magic);
30             dLat = (dLat*180.0)/((a*(1 - ee))/(magic*sqrtMagic)*pi);
31             dLon = (dLon*180.0)/(a/sqrtMagic*Math.Cos(radLat)*pi);
32             mgLat = wgLat + dLat;
33             mgLon = wgLon + dLon;
34         }
35 
36         private static bool OutOfChina(double lat, double lon)
37         {
38             if (lon < 72.004 || lon > 137.8347)
39                 return true;
40             if (lat < 0.8293 || lat > 55.8271)
41                 return true;
42             return false;
43         }
44 
45         private static double TransformLat(double x, double y)
46         {
47             double ret = -100.0 + 2.0*x + 3.0*y + 0.2*y*y + 0.1*x*y + 0.2*Math.Sqrt(Math.Abs(x));
48             ret += (20.0*Math.Sin(6.0*x*pi) + 20.0*Math.Sin(2.0*x*pi))*2.0/3.0;
49             ret += (20.0*Math.Sin(y*pi) + 40.0*Math.Sin(y/3.0*pi))*2.0/3.0;
50             ret += (160.0*Math.Sin(y/12.0*pi) + 320*Math.Sin(y*pi/30.0))*2.0/3.0;
51             return ret;
52         }
53 
54         private static double TransformLon(double x, double y)
55         {
56             double ret = 300.0 + x + 2.0*y + 0.1*x*x + 0.1*x*y + 0.1*Math.Sqrt(Math.Abs(x));
57             ret += (20.0*Math.Sin(6.0*x*pi) + 20.0*Math.Sin(2.0*x*pi))*2.0/3.0;
58             ret += (20.0*Math.Sin(x*pi) + 40.0*Math.Sin(x/3.0*pi))*2.0/3.0;
59             ret += (150.0*Math.Sin(x/12.0*pi) + 300.0*Math.Sin(x/30.0*pi))*2.0/3.0;
60             return ret;
61         }
62     }

 1 public struct EarthPoint
 2     {
 3         public const double Ea = 6378137; // 赤道半径 WGS84标准参考椭球中的地球长半径(单位:m)  
 4         public const double Eb = 6356725; // 极半径   
 5         public readonly double Ec;
 6         public readonly double Ed;
 7         public readonly double Jd;
 8         public readonly double Latidute;
 9         public readonly double Longitude;
10         public readonly double Wd;
11 
12         public EarthPoint(double longitude, double latidute)
13         {
14             Longitude = longitude;
15             Latidute = latidute;
16             Jd = Longitude*Math.PI/180; //转换成角度
17             Wd = Latidute*Math.PI/180; //转换成角度
18             Ec = Eb + (Ea - Eb)*(90 - Latidute)/90;
19             Ed = Ec*Math.Cos(Wd);
20         }
21 
22         public double Distance(EarthPoint point)
23         {
24             double dx = (point.Jd - Jd)*Ed;
25             double dy = (point.Wd - Wd)*Ec;
26             return Math.Sqrt(dx*dx + dy*dy);
27         }
28 
29         public static double GetDistance(
30             double latidute1,
31             double longitude1,
32             double latidute2,
33             double longitude2)
34         {
35             var p1 = new EarthPoint(longitude1, latidute1);
36             var p2 = new EarthPoint(longitude2, latidute2);
37             return p1.Distance(p2);
38         }
39     }