Shell 技巧 - 计算两个经纬度点之间的距离
上个月,我在这专栏的结尾留下了一个脚本,它可以返回两个地址的经纬度值,最终目的是让脚本计算这两个点之间的距离。 例如
$ 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 来说可能是一个有点太复杂的数学问题。
在花费了比我想承认的还要多的时间试图让它在 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 英里。
Dave Taylor 自 1980 年首次登录在线网络以来就一直参与 UNIX。 这意味着,是的,他现在即将迎来 30 周年。 您几乎可以在任何在线地方找到他,但请从这里开始:www.DaveTaylorOnline.com。 除了他的所有其他项目外,Dave 现在还是一名影评人。 您可以在 www.DaveOnFilm.com 阅读他的评论。