Ищем широту и долготу с помощью shell-скриптов

http://www.linuxjournal.com/article/10589

Work the Shell - Exploring Lat/Lon with Shell Scripts
November 1st, 2009 by Dave Taylor in HOWTOs

С развитием систем позиционирования на земной поверхности, связанными с мобильными устройствами, наиболее важным аспектом стало наличие удобного метода регистрации положения. Стандартный метод - показать координаты, то есть удаление от экватора и нулевого меридиана (географические широта и долгота). Стандартную нотацию понимают такие сервисы, как Google Maps, Yahoo Maps, MapQuest и подобные.

С точки зрения использования информации в shell-скриптах, нас будут интересовать как географические координаты, так и возможность, например, вычислить расстояние между двумя точками на земной поверхности.

Это может показаться довольно сложным делом, до тех пор, пока вы не узнаете о довольном простом API, предоставляемом Yahoo Maps. Этот API позволяет сформировать специальный URL, содержащий адрес объекта, передать его на сервер, взамен получить объект XML, содержащий значения широты и долготы нужного места.

Например, вы хотите ознакомиться с таким адресом: 1600 Pennsylvania Avenue, Washington, DC. Знаю, вы уже видели это место на картинках. Каковы его координаты?


$ u='http://api.maps.yahoo.com/ajax/geocode'
$ a='?appid=onestep&qt=1&id=m&qs=1600+pennsylvania+ave+washington+dc'
$ curl "$u$a"
YGeoCode.getMap({"GeoID"      : "m",
                 "GeoAddress" : "1600 pennsylvania ave washington dc",
                 "GeoPoint"   : {"Lat" : 38.89859,
                                 "Lon" : -77.035971},
                 "GeoMID"     : false,
                 "success"    : 1} ,1);
<!-- xm6.maps.re3.yahoo.com uncompressed/chunked
     Tue Aug  4 12:16:51 PDT 2009 -->

Обратите внимание на то, что в действительности вывод координат возвращается в виде двух строк. В дальнейших примерах этот вывод будет переформатирован для пущего удобства.

В данном объекте видим координаты: широта Latitude = 38.89859, долгота Longitude = -77.035971. Поместим эти данные в нужной форме (“38.89859,-77.035971”) в Google Maps и увидим картинку места (рис. 1).


Рисунок 1. Белый Дом

Это был адрес Белого Дома.

Начнём конструировать простой скрипт, принимающий адрес в формате улица-дом, город, штат и возвращающий координаты.

Первая часть проста: заберём ввод в командной строке и переформатируем его в URL-подобный синтаксис. Затем добавим к нему адрес Yahoo API и заберем полученный вывод:

Скрипт whereis.sh:


#!/bin/sh

url='http://api.maps.yahoo.com/ajax/geocode'
args='?appid=onestep&qt=1&id=m&qs='
converter="$url$args"

addr="$(echo $* | sed 's/ /+/g')"
curl -s "$converter$addr"
exit 0

Протестируем с помощью произвольных адресов:


$ sh whereis.sh 2001 Blake Street, Denver, CO
YGeoCode.getMap({"GeoID"      : "m",
                 "GeoAddress" : "2001 Blake Street, Denver, CO",
                 "GeoPoint"   : {"Lat" : 39.754386,
                                 "Lon" : -104.994261},
                 "GeoMID"     : false,
                 "success"    : 1}, 1);
<!-- x1.maps.sp1.yahoo.com uncompressed/chunked
     Tue Aug  4 12:37:44 PDT 2009 -->

Далее потребуется очистить вывод. Выделим широту и долготу, всё остальное отфильтруем и отбросим. Это можно проделать с помощью различных инструментов типа awk или Perl, но поскольку я учитель, буду использовать cut :).

Для этого подсчитаем двойные кавычки в блоках вывода. Итого: двенадцать кавычек до координат и пятнадцать кавычек после долготы. Учтём это дело:


$ sh whereis.sh 2001 Blake Street, Denver, CO | cut -d\" -f13-15
:39.754386,"Lon":-104.994261},

В целом хорошо, чтобы улучшить результат, будем использовать не диапазон полей, а отдельные поля:


$ sh whereis.sh 2001 Blake Street, Denver, CO | cut -d\" -f13,15
:39.754386,":-104.994261},

Это уже почти на 99% то, что нужно. Зачистим шум. Для этого вернёмся в сам скрипт и кое что поправим:


curl -s "$converter$addr" | \
    cut -d\" -f13,15 | \
    sed 's/[^0-9\.\,\-]//g'

Протестируем:


$ sh whereis.sh 2001 Blake Street, Denver, CO
39.754386,104.994261,

Почти. То есть очень, очень близко. Запятая в конце мешает. Хммм....

Добавим вторую подстановку в sed:


sed 's/[^0-9\.\,\-]//g;s/,$//'

Ну вот, получили то, что собирались получить в начале конструирования. Проверим на других адресах:


$ sh whereis.sh 1313 S. Disneyland Drive, Anaheim CA
33.814413,-117.924424

Ага, парковка в Диснейленде, Калифорния.

Теперь расстояние между двумя точками.

Это будет потруднее. Так как координаты угловые, потребуется некоторая математика. Я нашёл некий 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;

Тут радиус Земли R и он равен 6371 километр. Поскольку Земля не вполне шар, а некий сфероид, ошибки неизбежны. Однако, посмотрим, что же всё-таки получится.

Для вычислений в shell обычно используют bc, и для наших целей, видимо, его хватит за глаза. Несмотря на его некоторый, скажем, GNU-синтаксис.

Например, как задать число пи с помощью bc (это пример из документации bc - прим. перев.):


pi=$(echo "scale=10; 4*a(1)" | bc -l)

Поскольку bc работает с радианами, а не градусами, нам еще потребуется конвертация одного в другое. Обратим внимание также на то, что нам потребуется работать сразу с двумя адресами, четырьмя значениями:


$ sh farapart.sh \
  "1600 pennsylvania ave, washington dc" \
  "1313 s. disneyland drive, anaheim, ca"

Lat/long for 1600 pennsylvania ave, washington dc

= 38.89859, -77.035971

Lat/long for 1313 s. disneyland drive, anaheim, ca

= 33.814413, -117.924424

Конвертируем градусы в радианы.

Способ ясен:


Radians = degrees * ( pi / 180 )

Пи, естественно, равно 3.1415926535897932384.

Заданные значения координат:

41.878658, -87.640404

А тут значения, конвертированные в радианы:

0.7309204767, -1.529613605

Для примера вычислим одно значение с помощью bc:


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

Всё здорово, но только для вычисления дистанции потребуется задействовать функцию atan2(), которая не является частью bc.

Похоже, задачка немного слишком трудна для bc. После некоторого размышления и битья головой об стену в течение пары часов, я решил, что маленькая программка на С решает проблему.

Программа получает две пары значений координат в градусах и вычисляет расстояние между точками в милях.


Listing 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 миль. Это на 10% короче, чем даёт Гугль Мапс, но возможно это оттого, что нет прямого хода по дорогам?

Теперь у нас в наличии все детали. Соединим их вместе.

Синтетический скрипт получился довольно коротким за счет того, что в нем вызывается программа на С для вычисления дистанции.


#!/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

Его можно сделать еще короче, если научить С-программу принимать пары координат 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 Maps дает это расстояние как расстояние между двумя городами, и оно составляет по их алгоритму 3015 миль.

Придется отлаживать математику.

Видимо что-то не так в расчете. До сих пор мы вводили координаты и любовались на картинку, где это место расположено. Естественно, такой подход не дает гарантий, что в автоматическом режиме мы получим именно нужные координаты.

Отладку данной игрушки мы оставляем читателям. Присылайте ваши решения и лучшие из них мы опубликуем в следующем месяце.

Dave Taylor:
has been involved with UNIX since he first logged in to the on-line network in 1980. That means that, yes, he's coming up to the 30-year mark now. You can find him just about everywhere on-line, but start here: www.DaveTaylorOnline.com. In addition to all his other projects, Dave is now a film critic. You can read his reviews at www.DaveOnFilm.com.

Назад