Shell 技巧 - 计算两个经纬度点之间的距离

作者:Dave Taylor

上个月,我在这专栏的结尾留下了一个脚本,它可以返回两个地址的经纬度值,最终目的是让脚本计算这两个点之间的距离。 例如

$ farapart.sh "union station, denver co" \
      "union station, chicago il"
Calculating lat/long for union station, denver co
= 39.75288, -105.000473
And calculating lat/long for union station, chicago il
= 41.878658, -87.640404

实际上,计算距离的公式相当复杂。 这是我上个月展示的数学运算的 JavaScript 实现

var R    = 6371;       // kilometers
var dLat = (lat2-lat1);
var dLon = (lon2-lon1);
var a    = Math.sin(dLat/2) * Math.sin(dLat/2) +
           Math.cos(lat1.toRad()) * Math.cos(lat2.toRad()) *
           Math.sin(dLon/2) * Math.sin(dLon/2);
var c    = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d    = R * c;

不用说,将此转换为 shell 脚本会有点棘手,但由于我们必须使用比 Bash 内置功能更复杂的数学工具,这也意味着我们有许多选项可供选择,包括 Perl、awk 和 bc。 实际上,我们也可以编写一个快速的 C 程序来解决给定四个变量的这个方程,但真的,当我可以使其变得复杂时,为什么要使其变得容易呢? 如果我想要容易的,我会拿出一些 Perl,对吧? 上个月,我承诺了一些 bc,所以让我们看看我们是否可以让那个老旧的应用程序完成繁重的工作。

度到弧度

数学上的第一步是将我们从地图系统获得的经纬度值从度转换为弧度。 这结果证明是很简单的

弧度 = 度 * ( π / 180 )

π 当然是 3.1415926535897932384。

给定如下值

41.878658, -87.640404

那么,这些值的弧度等价物是

0.7309204767, -1.529613605

为了熟悉 bc,这里有一个简单的命令行方式来计算其中一个值

echo "scale=8; -87.640404 * ( 3.14159265 / 180)" | bc

这一切都很好,但事实证明,我探索的用于计算两点之间距离的不同方程需要 atan2() 函数,而该函数不是 bc 的一部分。

与其用头撞老式的墙壁直到头破血流,我不如认输并承认这对于 shell 脚本和 bc 来说可能是一个有点太复杂的数学问题。

Dave 认输了!

在花费了比我想承认的还要多的时间试图让它在 bc 中正常工作后,我将“认输”并暂时切换到另一种编程语言。 我将跳到 C 中写几行代码,并快速写出一个简单的程序,该程序在给定以度为单位的两个经纬度对的情况下,计算它们之间以英里为单位的距离(清单 1)。

清单 1. C 距离程序

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define EARTH_RADIUS       (6371.0072 * 0.6214)
#define TORADS(degrees)    (degrees * (M_PI / 180))

main(int argc, char **argv)
{
   double lat1, long1, lat2, long2;
   double dLat, dLong, a, c, d;

   lat1  = TORADS(atof(argv[1]));
   long1 = TORADS(atof(argv[2]));
   lat2  = TORADS(atof(argv[3]));
   long2 = TORADS(atof(argv[4]));

   dLat  = lat2 - lat1;
   dLong = long2 - long1;

   a = sin(dLat/2) * sin(dLat/2) +
       cos(lat1) * cos(lat2) * sin(dLong/2) * sin(dLong/2);
   c = 2 * atan2(sqrt(a), sqrt(1-a));

   printf("%g\n", EARTH_RADIUS * c);
}

它有效吗? 让我们看看

$ distance 39.75288 -105.000473 41.878658 -87.640404
917.984

这似乎是合理的。 这两点之间的大圆距离是 917 英里。 当然,Google 地图显示距离大约远 10%,但这也许是因为没有通过公路的直线距离路线?

当然,这个公式也存在误差,因为地球不是一个完美的球体,而是一个扁球体,其直径因测量位置而异。 但就我们的目的而言,让我们忽略这个问题。 您可以在 Google 上搜索以了解 Vincenty 公式之类的内容,但这超出了这个荒谬的题外话的范围。

现在我们拥有了我们需要的所有部分:位置到经纬度以及两个经纬度点之间的距离。 让我们把它们放在一起。

将它们整合在一起

为了使一切都正常工作,我实际上对原始脚本进行了修改和删减,使其更加简洁,当然,还调用了清单 1 中所示的 C “distance” 程序。 [清单 1 也可在我们的 FTP 站点 ftp://ftp.linuxjournal.com/pub/lj/listings/issue188/10606.tgz 上找到。] 准备好了吗? 它出奇的短

#!/bin/sh
converter="http://api.maps.yahoo.com/ajax/
↪geocode?appid=onestep&qt=1&id=m&qs="

tmpfile="/tmp/bc.script.$$"

# Get lat/long for point 1
addr="$(echo $1 | sed 's/ /+/g')"
values="$(curl -s $converter$addr | \
  cut -d\" -f13,15 | \
  sed 's/[^0-9\.\,\-]//g;s/,$//')"

lat1=$(echo $values | cut -d, -f1)
long1=$(echo $values | cut -d, -f2)

# Get lat/long for point 2
addr="$(echo $2 | sed 's/ /+/g')"
values="$(curl -s $converter$addr | \
  cut -d\" -f13,15 | \
  sed 's/[^0-9\.\,\-]//g;s/,$//')"

lat2=$(echo $values | cut -d, -f1)
long2=$(echo $values | cut -d, -f2)

# Now we have the lat/long for both points, let's
# figure out the distance between them...
dist=$(./distance $lat1 $long1 $lat2 $long2)
echo "$1 to $2 is $dist miles"
exit 0

如果我们调整 C 程序以接受 x,y 位置对,脚本会更短,但我会将这个留给您。 相反,让我们做一些测试

$ farapart.sh \
      "union station, denver, co" \
      "union station, chicago, il"
union station, denver, co to
    union station, chicago, il is 917.984 miles

现在,再来一些更模糊的例子如何

$ farapart.sh "long beach, ca" "boston, ma"
long beach, ca to boston, ma is 2597.53 miles

哎呀,该死,这似乎太短了。 让我们看看 Yahoo 地图报告的这两个城市之间的距离是多少。 果然,它报告说行程应该是 3,015 英里,而不是 2,597 英里。

调试数学公式

在某个地方存在一个错误,导致我们得到的结果不佳。 我猜想 C 程序中存在某种严重的舍入误差(因为我们可以通过实验验证我们获得的经纬度信息是有效的,只需将其插入地图应用程序并查看它将我们放置在哪里)。

但是,我已经在这个例子上耗尽了精力。 事实证明,它比我预期的要棘手得多,我将其作为练习留给您,亲爱的读者,看看您是否可以找出 C 程序中哪里出了问题,并将您的修复报告给我们。 我们将在下个月发布其中最好的! 同时,在下个月的专栏中,我将回到更关于 shell 而更少关于数学的内容。 我的意思是,见鬼,当我攻读计算机科学学位时,我不喜欢数学,所以我现在为什么要玩数学呢?

Dave Taylor 自 1980 年首次登录在线网络以来就一直参与 UNIX。 这意味着,是的,他现在即将迎来 30 周年。 您几乎可以在任何在线地方找到他,但请从这里开始:www.DaveTaylorOnline.com。 除了他的所有其他项目外,Dave 现在还是一名影评人。 您可以在 www.DaveOnFilm.com 阅读他的评论。

加载 Disqus 评论