使用awk查找公里中两个经纬度位置之间的距离

Question

我想从中间站计算两个站（前后）之间的距离。每个站都有一个 GPS 位置。例如：

P1      20.2    70
P2      21      70.3
P3      21.5    70.4
P4      22      71
P5      22.75   71.6
P6      23      72
P7      23.2    72.4
P8      24      73.3
P9      24.5    74
P10     25.1    74.3

这里第 1 列中的每个 P 是特定站点的站点名称，并分别在第 2 列和第 3 列中具有纬度和经度位置。

我想计算每个站点与其前面和后面的距离，即

For P2 what is the distance between P1 and P3 in KM?
For P3 What is the distance between P2 and P4 in KM?
For P4 what is the distance between P3 and P5 in KM?
........
........
For P9 what is the distance between P8 and P10 in KM?

所需的输出如下所示：

P2      149.62
P3      134.27
P4      190.60
P5      155.56
P6      100.97
P7      180.41
P8      226.77
P9      163.53

我使用以下公式计算了距离，但我不确定在基于 GPS 的位置中使用它是否正确。

d = sqrt(pow(lat2-lat1, 2) + pow(lon2-lon1, 2))

我正在为上述内容编写awk 代码，但遇到了问题。 Karafa（见答案）根据我的要求帮助编写了awk 脚本。这是

awk ' {k[NR] = $1
         a[NR] = $2
         b[NR] = $3}
    END {for(i=2; i<NR; i++)
           print k[i], (sqrt((a[i+1]-a[i-1])^2 + (b[i+1]-b[i-1])^2))*110}' file.txt

它工作得非常好，但是，在得到 Dawg 的评论（见 cmets）后，我意识到 Haversine 公式将是正确的方法，我从这里 https://rosettacode.org/wiki/Haversine_formula#AWK 处理。

# syntax: GAWK -f HAVERSINE_FORMULA.AWK
# converted from Python
BEGIN {
    distance(36.12,-86.67,33.94,-118.40) # BNA to LAX
    exit(0)
}
function distance(lat1,lon1,lat2,lon2,  a,c,dlat,dlon) {
    dlat = radians(lat2-lat1)
    dlon = radians(lon2-lon1)
    lat1 = radians(lat1)
    lat2 = radians(lat2)
    a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
    c = 2 * atan2(sqrt(a),sqrt(1-a))
    printf("distance: %.4f km\n",6372.8 * c)
}
function radians(degree) { # degrees to radians
    return degree * (3.1415926 / 180.)
} 

现在我想用上面的 Haversine 公式更新解决方案。我这样做如下，但遇到很多语法错误：

awk ' {k[NR] = $1
             a[NR] = $2
             b[NR] = $3}
        END {for(i=2; i<NR; i++)
    function distance(a[i-1],b[i-1],a[i+1],b[i+1],  a,c,dlat,dlon) {
        dlat = radians(a[i+1]-a[i-1])
        dlon = radians(b[i+1]-b[i-1])
        lat1 = radians(a[i-1])
        lat2 = radians(a[i+1])
        a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
        c = 2 * atan2(sqrt(a),sqrt(1-a))
        #printf("distance: %.4f km\n",6372.8 * c)  #I have disabled it
    }
    function radians(degree) { # degrees to radians
        return degree * (3.1415926 / 180.)
    }
               print k[i], 6372.8 * c}' file.txt

Answer 1

这将以最小的变化解决您的awk 问题，但对于距离，您使用的是欧几里得范数，这在球体上无效，但对于短距离来说是很好的近似值。

$ awk ' {k[NR] = $1
         a[NR] = $2
         b[NR] = $3}
    END {for(i=2; i<NR; i++)
           print k[i], sqrt((a[i+1]-a[i-1])^2 + (b[i+1]-b[i-1])^2)}' file

P2 1.36015
P3 1.22066
P4 1.73277
P5 1.41421
P6 0.917878
P7 1.64012
P8 2.06155
P9 1.48661

不确定您使用的是哪种缩放方式。
以 P2 为例

dist^2 = (21.5-20.2)^2 + (70.4-70)^2 = 1.3^2 + 0.4^2 = 1.69 + 0.16 = 1.85 = (1.36)^2

要使用基于度数的球面距离，只需插入正确的距离函数

$ awk 'function distance(lat1,lon1,lat2,lon2,  a,c,dlat,dlon) {
         dlat = radians(lat2-lat1)
         dlon = radians(lon2-lon1)
         lat1 = radians(lat1)
         lat2 = radians(lat2)
         a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
         c = 2 * atan2(sqrt(a),sqrt(1-a))
         return 6372.8 * c}
      function radians(degree) {return degree * (3.1415926 / 180.)}  
        {k[NR] = $1
         a[NR] = $2
         b[NR] = $3}
    END {for(i=2; i<NR; i++)
           print k[i], distance(a[i+1],b[i+1],a[i-1],b[i-1])}' file


P2 150.453
P3 132.736
P4 186.056
P5 151.428
P6 96.0023
P7 173.071
P8 217.711
P9 158.759