Определение расстояния между географическими точками в MySQL

При разработке современного сайта часто возникает необходимость реализовать функционал вывода близлежащих географических точек. Самым оптимальным способом решения этой задачи является перекладывание работы по реализации определения точек на плечи MySQL. Если конкретней, то нам будут нужны возможности пространственных расширений MySQL (до версии 5.0.16 эти расширения были доступны только для MyISAM, более поздние версии MySQL поддерживают работу пространственных расширений с InnoDB, NDB, BDB и ARCHIVE).

Расстояние между точками будет вычисляться по формуле Хаверсина ^[1]. Формула позволяет получать расстояние между точками с очень низкой погрешностью (величина погрешности прямо пропорциональна расстоянию между точками, и не превышает 10-20 километров при вычислении очень больших расстояний, например между штаб-квартирой Google в Калифорнии (37.422045, -122.084347) и оперным театром в Сиднее, Австралия (-33.856553, 151.214696)).

Для корректного вычисления близлежащих точек с расстоянием, понадобится создать следующие хранимые функции и процедуры:

функция geodist

geodist() определяет расстояние между точками по их координатам.

-- число 6367 - это радиус Земли в километрах, для использования решения в других единицах измерения, достаточно сконвертировать радиус Земли в нужную единицу измерения
DELIMITER $$
DROP FUNCTION IF EXISTS geodist $$
CREATE FUNCTION geodist (
  src_lat DECIMAL(9,6), src_lon DECIMAL(9,6),
  dst_lat DECIMAL(9,6), dst_lon DECIMAL(9,6)
) RETURNS DECIMAL(6,2) DETERMINISTIC
BEGIN
 SET @dist := 6367 * 2 * ASIN(SQRT(
      POWER(SIN((src_lat - ABS(dst_lat)) * PI()/180 / 2), 2) +
      COS(src_lat * PI()/180) *
      COS(ABS(dst_lat) * PI()/180) *
      POWER(SIN((src_lon - dst_lon) * PI()/180 / 2), 2)
    ));
 RETURN @dist;
END $$
DELIMITER ;

функция geodist_pt

geodist_pt() является оберткой для geodist(), и работает с координатами точек в виде объекта типа POINT.

DELIMITER $$
DROP FUNCTION IF EXISTS geodist_pt $$
CREATE FUNCTION geodist_pt (src POINT, dst POINT) 
RETURNS DECIMAL(6,2) DETERMINISTIC
BEGIN
  RETURN geodist(X(src), Y(src), X(dst), Y(dst));
END $$
DELIMITER ;

процедура geobox_pt

Вычисляет координаты верхнего левого и нижнего правого угла области поиска с помощью процедуры geobox(), затем конвертирует полученные координаты в объекты типа POINT.

-- pt -> центральная точка области поиска
-- dist -> расстояние от центра в километрах
-- top_lft -> верхний левый угол области поиска (объект типа POINT)
-- bot_rgt -> нижний правый угол области поиска (объект типа POINT)
DELIMITER $$
DROP PROCEDURE IF EXISTS geobox_pt $$
CREATE PROCEDURE geobox_pt (
    IN pt POINT, IN dist DECIMAL(6,2),
    OUT top_lft POINT, OUT bot_rgt POINT
) DETERMINISTIC
BEGIN

  CALL geobox(X(pt), Y(pt), dist, @lat_top, @lon_lft, @lat_bot, @lon_rgt);
  SET top_lft := POINT(@lat_top, @lon_lft);
  SET bot_rgt := POINT(@lat_bot, @lon_rgt);
END $$
DELIMITER ;

процедура geobox

Вычисляет координаты области поиска

-- src_lat, src_lon -> центральная точка области поиска
-- dist -> расстояние от центра в километрах
-- lat_top, lon_lft -> координаты верхнего левого угла области поиска
-- lat_bot, lon_rgt -> координаты нижнего правого угла области поиска
DELIMITER $$
DROP PROCEDURE IF EXISTS geobox $$
CREATE PROCEDURE geobox (
  IN src_lat DECIMAL(9,6), IN src_lon DECIMAL(9,6), IN dist DECIMAL(6,2),
  OUT lat_top DECIMAL(9,6), OUT lon_lft DECIMAL(9,6),
  OUT lat_bot DECIMAL(9,6), OUT lon_rgt DECIMAL(9,6)
) DETERMINISTIC
BEGIN
  SET lat_top := src_lat - (dist / 69);
  SET lon_lft := src_lon - (dist / ABS(COS(RADIANS(src_lat)) * 69));
  SET lat_bot := src_lat + (dist / 69);
  SET lon_rgt := src_lon + (dist / ABS(COS(RADIANS(src_lat)) * 69));
END $$
DELIMITER ;

Пример использования:

Допустим, что есть таблица, в которой хранятся географические данные. Таблица имеет такую структуру:

CREATE TABLE geo (
         id INT,
         name VARCHAR(100),
         x DECIMAL(9,6),
         y DECIMAL (9,6)
);

тогда, для получения всех населенных пунктов, расположенных на расстоянии 200 километров от исходной точки понадобится сделать следующее:

-- в переменную src сохраняем координаты нужной точки в виде объекта POINT (системная функция Point() создает объект типа POINT из двух строковых координат)
SELECT @src := Point(x,y) FROM geo WHERE name = 'Москва';
-- формируем "область поиска" с заданным радиусом
CALL geobox_pt(@src, 200.0, @top_lft, @bot_rgt);
-- достаем данные 
SELECT g.name, geodist(X(@src), Y(@src), x, y) AS dist
FROM geo g
WHERE x BETWEEN X(@top_lft) AND X(@bot_rgt)
AND y BETWEEN Y(@top_lft) AND Y(@bot_rgt)
HAVING dist < 200.0
ORDER BY dist desc;

Автор: xoyk

Источник ^[2]