【191】Java8在大比例尺小范围地图上，根据wgs84坐标系的经纬度计算两个点之间的方向和距离

场景

本文代码在大比例迟、小范围的地图上测试过。这些地图一般是县、区、镇、街道等范围的，其测试效果较好。由于地图范围较小，可以把经纬度近似看作直线。

问题分析

方向一共分东、南、西、北、东北、西北、西南、东南共八个方向。一周是 360 度，360 度除以 8 等于 45 度。以输入的第一个点为原点，绕此点一周，每个方向占45度。如果第二个点和第一个点的线段落在对应的角度范围内，就是对应的方向。

上图中表示了方向和角度的关系。第一个点是A，第二个点是B，如图所示，B在A的东北方向。

我以第一个点A为原点，点A所在的纬度线为x轴，点A所在的经度线为y轴，可以把地图划分成四个象限。在每个象限内，第二个点B和第一个点A之间的线段与y轴的角度决定了第二个点的方向。

请看下图例子：

上图中以第三象限为例，标记出了对应的角度和方向的关系。

怎么求得角度的？

设第一个点是A，第二个点是B。以A为原点构建平面直角坐标系（如上图）。过B作一条与y轴垂直的直线，交点是C。显然ABC组成了一个直角三角形。利用经纬度可获取AB长度和AC长度，利用三角函数可计算角度。

代码实现

package zhangchao;

public class GeoUtils {
    /**
     * 输入两个点的wgs84坐标系的经纬度，计算距离，单位是米。
     * @param longitude1 第一个点的经度
     * @param latitude1  第一个点的纬度
     * @param longitude2 第二个点的经度
     * @param latitude2  第二个点的纬度
     * @return 两个点的距离，单位是米。
     */
    public static double calWgs84Distance(double longitude1, double latitude1,
                                          double longitude2, double latitude2) {
        double a = 6378137, b = 6356752.3142, f = 1 / 298.257223563;
        double L = Math.toRadians(longitude2 - longitude1);
        double U1 = Math.atan((1 - f) * Math.tan(Math.toRadians(latitude1)));
        double U2 = Math.atan((1 - f) * Math.tan(Math.toRadians(latitude2)));
        double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
        double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
        double lambda = L, lambdaP, iterLimit = 100;
        double cosSqAlpha;
        double sinSigma;
        double cos2SigmaM;
        double sigma;
        double sinLambda;
        double cosLambda;
        double cosSigma;
        do {
            sinLambda = Math.sin(lambda);
            cosLambda = Math.cos(lambda);
            sinSigma = Math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda));
            if(sinSigma == 0)
                return 0;
            cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
            sigma = Math.atan2(sinSigma, cosSigma);
            double sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
            cosSqAlpha = 1 - sinAlpha * sinAlpha;
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
//            if(isNaN(cos2SigmaM))
//                cos2SigmaM = 0;
            double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
            lambdaP = lambda;
            lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
        } while (Math.abs(lambda-lambdaP) > (1e-12) && --iterLimit>0);
        if(iterLimit == 0) {
            return -1;
        }
        double uSq = cosSqAlpha * (a * a - b * b) / (b * b);
        double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
        double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
        double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
        double s = b * A * (sigma - deltaSigma);
        double fwdAz = Math.atan2(cosU2 * sinLambda, cosU1 * sinU2 - sinU1 * cosU2 * cosLambda);
        double revAz = Math.atan2(cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda);
        return s;
    }


    /**
     * 按照wgs84坐标系输入两个点的经纬度，计算第二个点在第一个点的什么方向。
     * 方向用 东、南、西、北、东北、西北、西南、东南
     *
     * @param longitude1 第一个点的经度
     * @param latitude1  第一个点的纬度
     * @param longitude2 第二个点的经度
     * @param latitude2  第二个点的纬度
     * @return 返回 东、南、西、北、东北、西北、西南、东南
     */
    public static String calNSEW(double longitude1, double latitude1,
                          double longitude2, double latitude2) {
        if (longitude1 == longitude2 && latitude2 > latitude1) {
            return "北";
        }
        if (longitude1 == longitude2 && latitude2 < latitude1) {
            return "南";
        }
        if (latitude1 == latitude2 && longitude2 > longitude1) {
            return "东";
        }
        if (latitude1 == latitude2 && longitude2 < longitude1) {
            return "西";
        }
        // 计算两个点之间的距离
        double dis12 = calWgs84Distance(longitude1, latitude1, longitude2, latitude2);
        // 以第一个点为原点，第一个点所在的纬度线为x轴，第一个点所在的经度线为y轴，可以把地图划分成四个象限。
        // 第一象限
        if (longitude2 > longitude1 && latitude2 > latitude1) {
            double longitude3 = longitude1;
            double latitude3  = latitude2;
            double dis13 = calWgs84Distance(longitude1, latitude1, longitude3, latitude3);
            double radians = Math.acos(dis13 / dis12);
            double angle =Math.toDegrees(radians);
            if (angle < 22.5) {
                return "北";
            } else if (angle >= 22.5 && angle <= 67.5) {
                return "东北";
            } else {
                return "东";
            }
        }
        // 第二象限
        if (longitude2 < longitude1 && latitude2 > latitude1) {
            double longitude3 = longitude1;
            double latitude3  = latitude2;
            double dis13 = calWgs84Distance(longitude1, latitude1, longitude3, latitude3);
            double radians = Math.acos(dis13 / dis12);
            double angle =Math.toDegrees(radians);
            if (angle < 22.5) {
                return "北";
            } else if (angle >= 22.5 && angle <= 67.5) {
                return "西北";
            } else {
                return "西";
            }
        }
        // 第三象限
        if (longitude2 < longitude1  && latitude2 < latitude1) {
            double longitude3 = longitude1;
            double latitude3  = latitude2;
            double dis13 = calWgs84Distance(longitude1, latitude1, longitude3, latitude3);
            double radians = Math.acos(dis13 / dis12);
            double angle =Math.toDegrees(radians);
            if (angle < 22.5) {
                return "南";
            } else if (angle >= 22.5 && angle <= 67.5) {
                return "西南";
            } else {
                return "西";
            }
        }
        // 第四象限
        if (longitude2 > longitude1  && latitude2 < latitude1) {
            double longitude3 = longitude1;
            double latitude3  = latitude2;
            double dis13 = calWgs84Distance(longitude1, latitude1, longitude3, latitude3);
            double radians = Math.acos(dis13 / dis12);
            double angle =Math.toDegrees(radians);
            if (angle < 22.5) {
                return "南";
            } else if (angle >= 22.5 && angle <= 67.5) {
                return "东南";
            } else {
                return "东";
            }
        }
        return null;
    }


	public static void main(String[] args) {
        double lon1 = 120.084856;
        double lat1 = 35.891211;
        double lon2 = 120.444538;
        double lat2 = 36.390704;
        String nsew = calNSEW(lon1, lat1, lon2, lat2);
        double distance = calWgs84Distance(lon1, lat1, lon2, lat2);
        System.out.println(nsew);
        System.out.println(distance);
    }
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161