经纬度点距离的那点儿事

前段时间整地址解析和道路解析一直会遇到经纬度点到点的距离以及一个点到一个道路的距离的问题,现在把相关的工具代码及用法贴出来做个备忘:

1.点到线段的距离,其中PointD只是一个坐标点的结构体而已

 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         }

2.点到点的距离

 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         }

以上已经说了坐标系中,点到点的距离以及点到线的距离,那么问题来了,经纬度转坐标系怎么转呢。。。。

以下是经纬度转高斯坐标,用法就是GeoToGauss(经度,纬度,0,0,ref x,ref y,MLP);其中x,y是要转换出来的变量,MLP这个是比如跟线段比较的话,那么MLP=(lon1+lon2)/2;

 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         }

以上代码经过验证,误差不大,其实直接用经纬度,不转坐标系也没啥大问题。

注:以上仅个人技术见解,如果有更好的方法或者对以上方式有修改意见的,希望大家提出来,一起验证修改

原文地址:https://www.cnblogs.com/djzny/p/4111401.html