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

GIS算法总结(求经纬度,求距离,求方位角)

程序员文章站 2022-06-11 16:10:36
...

公司项目中用到的一些GIS算法,记录一下。所有经纬度都是基于WGS-84标准,为Google Earth上的经纬度。
1.已知两个点的经纬度,求它们的直线距离

python代码实现

from math import *
def getDistance(latA, lonA, latB, lonB):
    ra = 6378140  # radius of equator: meter
    rb = 6356755  # radius of polar: meter
    flatten = (ra - rb) / ra  # Partial rate of the earth
    # change angle to radians
    radLatA = radians(latA)
    radLonA = radians(lonA)
    radLatB = radians(latB)
    radLonB = radians(lonB)
 
    pA = atan(rb / ra * tan(radLatA))
    pB = atan(rb / ra * tan(radLatB))
    x = acos(sin(pA) * sin(pB) + cos(pA) * cos(pB) * cos(radLonA - radLonB))
    c1 = (sin(x) - x) * (sin(pA) + sin(pB))**2 / cos(x / 2)**2
    c2 = (sin(x) + x) * (sin(pA) - sin(pB))**2 / sin(x / 2)**2
    dr = flatten / 8 * (c1 - c2)
    distance = ra * (x + dr)
    return distance
    
distance = getDistance(30.25482684590,120.01375138640,30.25582684590,120.01375138640)
print(distance, '米')

2.已知两个点的经纬度,求它们的方位角

python代码实现

def getDegree(latA, lonA, latB, lonB):
    """
    Args:
        point p1(latA, lonA)
        point p2(latB, lonB)
    Returns:
        bearing between the two GPS points,
        default: the basis of heading direction is north
    """
    radLatA = radians(latA)
    radLonA = radians(lonA)
    radLatB = radians(latB)
    radLonB = radians(lonB)
    dLon = radLonB - radLonA
    y = sin(dLon) * cos(radLatB)
    x = cos(radLatA) * sin(radLatB) - sin(radLatA) * cos(radLatB) * cos(dLon)
    brng = degrees(atan2(y, x))
    brng = (brng + 360) % 360
    return brng
    
brng = getDegree(30.25482684590,120.01375138640,30.25582684590,120.01375138640)
print(brng)

3.已知一个点的经纬度,方位角,距离,求另一点的经纬度

python代码实现

from math import *


def rad(d):
    return d * pi / 180.0


def deg(x):
    return x * 180 / pi


def getLonAndLat(lat, lng, brng, dist):
    a = 6378137
    b = 6356752.3142
    f = 1.0 / 298.257223563

    lon1 = lng * 1
    lat1 = lat * 1
    s = dist
    alpha1 = rad(brng)
    sinAlpha1 = sin(alpha1)
    cosAlpha1 = cos(alpha1)

    tanU1 = (1 - f) * tan(rad(lat1))
    cosU1 = 1 / sqrt((1 + tanU1 * tanU1))
    sinU1 = tanU1 * cosU1
    sigma1 = atan2(tanU1, cosAlpha1)
    sinAlpha = cosU1 * sinAlpha1
    cosSqAlpha = 1 - sinAlpha * sinAlpha
    uSq = cosSqAlpha * (a * a - b * b) / (b * b)
    A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))

    sigma = s / (b * A)
    sigmaP = 2 * pi
    while abs(sigma - sigmaP) > 1e-12:
        cos2SigmaM = cos(2 * sigma1 + sigma)
        sinSigma = sin(sigma)
        cosSigma = cos(sigma)
        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 = s / (b * A) + deltaSigma

    tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1
    lat2 = atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
                 (1 - f) * sqrt(sinAlpha * sinAlpha + tmp * tmp))
    lambdA = atan2(sinSigma * sinAlpha1, cosU1*cosSigma - sinU1 * sinSigma * cosAlpha1)
    C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
    L = lambdA - (1-C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM+C * cosSigma * (-1+2 * cos2SigmaM * cos2SigmaM)))
    revAz = atan2(sinAlpha, -tmp)
    lngLatObj = {'latitude': deg(lat2), 'longitude': lon1 + deg(L), }
    return lngLatObj


print(getLonAndLat(30.2548268459,120.0137513864,0,110.85670061685299))

JavaScript代码实现

            /**
             * 度换成弧度
             * @param  {Float} d  度
             * @return {[Float}   弧度
             */
            function rad (d){
                return d * Math.PI / 180.0;
            };
            /**
             * 弧度换成度
             * @param  {Float} x 弧度
             * @return {Float}   度
             */
            function deg (x){
                return x*180/Math.PI;
            };
            /**
             * 
             * @param {*} lng 经度 113.3960698
             * @param {*} lat 纬度 22.941386
             * @param {*} brng 方位角 45   ---- 正北方:000°或360°  正东方:090° 正南方:180°  正西方:270°
             * @param {*} dist 距离 9000
             * 
            */
            function getLonAndLat (lat,lng,brng,dist){
                //大地坐标系资料WGS-84 长半径a=6378137 短半径b=6356752.3142 扁率f=1/298.2572236
                var a=6378137; 
                var b=6356752.3142; 
                var f=1/298.257223563;

                var lon1 = lng*1;
                var lat1 = lat*1;
                var s = dist;
                var alpha1 = rad(brng);
                var sinAlpha1 = Math.sin(alpha1);
                var cosAlpha1 = Math.cos(alpha1);

                var tanU1 = (1-f) * Math.tan(rad(lat1));
                var cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1*cosU1;
                var sigma1 = Math.atan2(tanU1, cosAlpha1);
                var sinAlpha = cosU1 * sinAlpha1;
                var cosSqAlpha = 1 - sinAlpha*sinAlpha;
                var uSq = cosSqAlpha * (a*a - b*b) / (b*b);
                var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
                var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));

                var sigma = s / (b*A), sigmaP = 2*Math.PI;
                while (Math.abs(sigma-sigmaP) > 1e-12) {
                    var cos2SigmaM = Math.cos(2*sigma1 + sigma);
                    var sinSigma = Math.sin(sigma);
                    var cosSigma = Math.cos(sigma);
                    var 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 = s / (b*A) + deltaSigma;
                }

                var tmp = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1;
                var lat2 = Math.atan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1,
                    (1-f)*Math.sqrt(sinAlpha*sinAlpha + tmp*tmp));
                var lambda = Math.atan2(sinSigma*sinAlpha1, cosU1*cosSigma - sinU1*sinSigma*cosAlpha1);
                var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
                var L = lambda - (1-C) * f * sinAlpha *
                    (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));

                var revAz = Math.atan2(sinAlpha, -tmp);  // final bearing
                
                var lngLatObj = {latitude:deg(lat2),longitude:lon1+deg(L),}
                return lngLatObj;
            }
            //调用
            console.log(getLonAndLat(30.2548268459,120.0137513864,0,110.85670061685299));

相关标签: GIS

上一篇: 对CIS的认识

下一篇: arcgis文档地址