经纬度点距离的那点儿事

 1         /// <summary>
 2         /// 点到线段的距离公式(利用平行四边形的面积算法)
 3         /// </summary>
 4         /// <param name="P">目标点</param>
 5         /// <param name="A">线段端点A</param>
 6         /// <param name="B">线段端点B</param>
 7         /// <returns></returns>
 8         public double DistancePointToSegment(PointD P, PointD A, PointD B)
 9         {
10             //计算点到线段(a,b)的距离  
11             //double l = 0.0;
12             //double s = 0.0;
13             //l = DistancePointToPoint(A, B);
14 
15             //s = ((A.Y - P.Y) * (B.X - A.X) - (A.X - P.X) * (B.Y - A.Y)) / (l * l);
16 
17             //return (Math.Abs(s * l));
18 
19             double a, b, c;
20             a = DistancePointToPoint(B, P);
21             if (a <= 0.00001)
22                 return 0.0;
23 
24             b = DistancePointToPoint(A, P);
25             if (b < 0.00001)
26                 return 0.0;
27 
28             c = DistancePointToPoint(A, B);
29             if (c <= 0.00001)
30                 return a;
31 
32             if (a * a >= b * b + c * c)
33                 return b;
34 
35             if (b * b >= a * a + c * c)
36                 return a;
37 
38             double l = (a + b + c) / 2;
39             double s = Math.Sqrt(Math.Abs(l * (l - a) * (l - b) * (l - c)));
40 
41             return 2 * s / c;
42         }

 1         /// <summary>
 2         /// 点到点的距离
 3         /// </summary>
 4         /// <param name="ptA"></param>
 5         /// <param name="ptB"></param>
 6         /// <returns></returns>
 7         private double DistancePointToPoint(PointD ptA, PointD ptB)
 8         {
 9             return Math.Sqrt(Math.Pow(ptA.X - ptB.X, 2) + Math.Pow(ptA.Y - ptB.Y, 2));
10         }

 1         private const double B_E2 = 0.0067385254147;
 2         private const double B_C = 6399698.90178271;
 3 
 4 
 5         /// <summary>
 6         /// 将地理坐标转换成绝对的高斯坐标
 7         /// </summary>
 8         /// <param name="jd">经度(度度格式)</param>
 9         /// <param name="wd">经度(度度格式)</param>
10         /// <param name="dh">带号</param>
11         /// <param name="dh_width">带宽(度)</param>
12         /// <param name="y">y坐标(米)</param>
13         /// <param name="x">x坐标(米)</param>
14         /// <param name="lp">中央子午线(度)</param>
15         public void GeoToGauss(double jd, double wd, short dh, short dh_width,ref  double y, ref double x, double lp)
16         {
17             double t;     //  t=tgB
18             double L;     //  中央经线的经度
19             double l0;    //  经差
20             double jd_hd, wd_hd;  //  将jd、wd转换成以弧度为单位
21             double et2;    //  et2 = (e‘ ** 2) * (cosB ** 2)
22             double N;     //  N = C / sqrt(1 + et2)
23             double X;     //  克拉索夫斯基椭球中子午弧长
24             double m;     //  m = cosB * PI/180 * l0 
25             double tsin, tcos;   //  sinB,cosB
26 
27             jd = jd * 3600;
28             wd = wd * 3600;
29 
30 
31             jd_hd = jd / 3600.0 * Math.PI / 180.0;    // 将以秒为单位的经度转换成弧度
32             wd_hd = wd / 3600.0 * Math.PI / 180.0;    // 将以秒为单位的纬度转换成弧度
33 
34             //jd_hd = jd  * Math.PI / 180.0;    // 将以度为单位的经度转换成弧度
35             //wd_hd = wd  * Math.PI / 180.0;    // 将以度为单位的纬度转换成弧度
36 
37             // 如果不设中央经线（缺省参数: -1000）,则计算中央经线,
38             // 否则,使用传入的中央经线,不再使用带号和带宽参数
39             //L = (DH - 0.5) * DH_width ;      // 计算中央经线的经度
40             if (lp == -1000)
41             {
42                 L = (dh - 0.5) * dh_width;      // 计算中央经线的经度
43             }
44             else
45             {
46                 L = lp;
47             }
48 
49             l0 = jd / 3600.0 - L;       // 计算经差
50             tsin =Math.Sin(wd_hd);        // 计算sinB
51             tcos = Math.Cos(wd_hd);        // 计算cosB
52             // 计算克拉索夫斯基椭球中子午弧长X
53             X = 111134.8611 / 3600.0 * wd - (32005.7799 * tsin + 133.9238 * Math.Pow(tsin, 3)
54                 + 0.6976 * Math.Pow(tsin, 5) + 0.0039 * Math.Pow(tsin, 7)) * tcos;
55             et2 = B_E2 * Math.Pow(tcos, 2);      //  et2 = (e‘ ** 2) * (cosB ** 2)
56             N = B_C / Math.Sqrt(1 + et2);      //  N = C / sqrt(1 + et2)
57             t = Math.Tan(wd_hd);         //  t=tgB
58             m = Math.PI / 180 * l0 * tcos;       //  m = cosB * PI/180 * l0 
59             x = X + N * t * (0.5 * Math.Pow(m, 2)
60                 + (5.0 - Math.Pow(t, 2) + 9.0 * et2 + 4 * Math.Pow(et2, 2)) * Math.Pow(m, 4) / 24.0
61                 + (61.0 - 58.0 * Math.Pow(t, 2) + Math.Pow(t, 4)) * Math.Pow(m, 6) / 720.0);
62             y = N * (m + (1.0 - Math.Pow(t, 2) + et2) * Math.Pow(m, 3) / 6.0
63                 + (5.0 - 18.0 * Math.Pow(t, 2) + Math.Pow(t, 4) + 14.0 * et2
64                 - 58.0 * et2 * Math.Pow(t, 2)) * Math.Pow(m, 5) / 120.0);
65         }