欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

通过两点经纬度和范围距离计算范围四个顶点的经纬度

程序员文章站 2022-07-05 18:05:00
...

   参考博文:https://blog.csdn.net/zengmingen/article/details/68490497

private static final double mile = 1609.344;//米-》英里换算
//通过起点、终点的经纬度求方位角
double azimuth = gps2d(originLat, originLong, destinationLat, destinationLong);
double range = 50;//范围 单位:英里mile
double dist = (range+1)*mile;//由于计算顶点经纬度的方法存在误差,所以距离需要+1
//求范围四个顶点的经纬度
Map<String,Object> origin1 = computerThatLonLat(originLong, originLat, azimuth+90, dist);
Map<String,Object> origin2 = computerThatLonLat(originLong, originLat, azimuth-90, dist);
Map<String,Object> destination1 = computerThatLonLat(destinationLong, destinationLat, azimuth+90, dist);
Map<String,Object> destination2 = computerThatLonLat(destinationLong, destinationLat, azimuth-90, dist);
    /* 
     * 大地坐标系资料WGS-84 长半径a=6378137 短半径b=6356752.3142 扁率f=1/298.2572236 
     */  
    /** 长半径a=6378137 */  
    private static double a = 6378137;  
    /** 短半径b=6356752.3142 */  
    private static double b = 6356752.3142;  
    /** 扁率f=1/298.2572236 */  
    private static double f = 1 / 298.2572236;  
  
    /** 
     * 经纬度换成弧度 (°)  
     * @param d 度 
     * @return 弧度 
     */  
    private static double rad(double d) {  
        return d * Math.PI / 180.0;  
    } 
  
    /** 
     * 弧度换成度  
     * @param x 弧度 
     * @return 度 
     */  
    private static double deg(double x) {  
        return x * 180 / Math.PI;  
    } 
    /** 
     * 通过一个点的经纬度、距离、方位角,计算另一点经纬度
     * @param lon 经度 
     * @param lat 纬度 
     * @param lonlat 已知点经纬度 
     * @param brng 方位角 
     * @param dist 距离(米) 
     */  
    public static Map<String,Object> computerThatLonLat(double lon, double lat, double brng, double dist) {
    	Map<String,Object> map = new HashMap<String,Object>();
    	double alpha1 = rad(brng);  
        double sinAlpha1 = Math.sin(alpha1);  
        double cosAlpha1 = Math.cos(alpha1);    
        double tanU1 = (1 - f) * Math.tan(rad(lat));  
        double cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1));  
        double sinU1 = tanU1 * cosU1;  
        double sigma1 = Math.atan2(tanU1, cosAlpha1);  
        double sinAlpha = cosU1 * sinAlpha1;  
        double cosSqAlpha = 1 - sinAlpha * sinAlpha;  
        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 cos2SigmaM=0;  
        double sinSigma=0;  
        double cosSigma=0;  
        double sigma = dist / (b * A), sigmaP = 2 * Math.PI;  
        while (Math.abs(sigma - sigmaP) > 1e-12) {  
            cos2SigmaM = Math.cos(2 * sigma1 + sigma);  
            sinSigma = Math.sin(sigma);  
            cosSigma = Math.cos(sigma);  
            double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)  
                    - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));  
            sigmaP = sigma;  
            sigma = dist / (b * A) + deltaSigma;  
        }  
  
        double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;  
        double lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,  
                (1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));  
        double lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);  
        double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));  
        double L = lambda - (1 - C) * f * sinAlpha  
                * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));  
  
        double revAz = Math.atan2(sinAlpha, -tmp); // final bearing  
        double new_lng = lon+deg(L);
        double new_lat = deg(lat2);
        map.put("longitude", new_lng);
        map.put("latitude", new_lat);
        return map;
    }    
    
    /**
     * 经纬度方位角计算
     * @param lat_a 起点纬度
     * @param lng_a 起点经度
     * @param lat_b 终点纬度
     * @param lng_b 终点经度
     * @return 方位角度数
     */
    public static double gps2d(double lat_a, double lng_a, double lat_b, double lng_b) {
        double d = 0;
        lat_a=lat_a*Math.PI/180;
        lng_a=lng_a*Math.PI/180;
        lat_b=lat_b*Math.PI/180;
        lng_b=lng_b*Math.PI/180;     
        d=Math.sin(lat_a)*Math.sin(lat_b)+Math.cos(lat_a)*Math.cos(lat_b)*Math.cos(lng_b-lng_a);
        d=Math.sqrt(1-d*d);
        d=Math.cos(lat_b)*Math.sin(lng_b-lng_a)/d;
        d=Math.asin(d)*180/Math.PI;
        return d;
     }