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

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

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


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

функция geodist

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

-- число 6371 - это радиус Земли в километрах, для использования решения в других единицах измерения, достаточно сконвертировать радиус Земли в нужную единицу измерения
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 := 6371 * 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;
Ads
AdBlock has stolen the banner, but banners are not teeth — they will be back

More

Comments 41

    +5
    А где тесты скорости? Думаю, при нескольких миллионах записей скорость будет не ахти. К слову, в своё время с проблемой скорости справился используя предварительную выборку, типа:

    latitude - 2
    longitude - 2
    
    ...  WHERE latitude > 1 AND latitude < 3 AND longitude > 1 AND longitude < 3
    

    А потом уже искал более точными алгоритмами, если это требовалось.
      0
      К слову, в статье это есть
      WHERE x BETWEEN X(@top_lft) AND X(@bot_rgt)
      AND y BETWEEN Y(@top_lft) AND Y(@bot_rgt)
      
        0
        Да, но там есть и
        SELECT g.name, geodist(X(@src), Y(@src), x, y) AS dist
        

        т.е. функция вызовется всех записей в таблице, и только потом уже выполнится условие.
          +1
          Проверил, на MySQL нет разницы в производительности между

          select a
          from t
          where a > ? and a < ? and a*a < ?
          

          и
          select a, a*a as s
          from t
          where a > ? and a < ?
          having s < ?
          


          Планировщик показывает одинаковую картину, при наличии индекса по «a»
            0
            На каких объёмах тестировали? Да и a*a — не сильно сложный пример.
              +1
              На малых. Цель понять принципиальную возможность ограничения области поиска с помощью индекса. Такая возможность подтверждается показаниями explain в обоих случаях.

              a*a — достаточный пример неиндексируемого выражения, чтобы проверить описанное выше.

              Скорость вычисления расстояния с помощью хранимой процедуры — это отдельный вопрос, не имеющий отношения к вашему замечанию о «предварительной выборке».
                0
                В комментарии ниже я приводил пример с сложной формулой вычисления. Она не является хранимой процедурой, но суть, в данном случае, та же.
                –1
                т.е. утверждение «функция вызовется всех записей в таблице» неверное
        0
        Удобное решение, особенно для генерации отчетов или промежуточных данных для нечастых операций. Можно спокойно применять в desktop-приложениях и в мобильных приложениях (если, конечно, БД поддерживает функции). На его основе можно также сделать и вычисление длины более сложных путей.

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

          Вопрос к знатокам: а оптимальнее ли перекладывать это все на БД? Я остановился на варианте делаем выборку слабым алгоритмом точек 200 и их уже вне БД обрабатываем более точными алгоритмами, выбирая оптимальные точки.
            +2
            а не проще в этом случае перейти на PostgreSQL или MongoDB?
            у них обоих есть расширения для индексации 2d и 2dSphere координат, выборка в итоге достаточна шустрая, без лишних телодвижений
              –2
              С PostgreSQL не сталкивался, а в MongoDB это хорошо в случае, когда выборка только по координатам. Если есть какие-то ещё условия выборки — работает не так уж и быстро.
                0
                Можете привести конкретный пример и объем данных?
                  0
                  База вида

                  | lat | lng | param1 | param2 | param3 |

                  Нужно найти места в районе 10 км + формула для подсчёта веса параметров, например (param1 + param2)/param3 (формула более сложная, параметров около 5, точно уже не помню).
                  Объём бд — 3+ миллиона записей и постоянно пополняется. Хранить вес не возможно, так как параметры очень часто изменяются.

                  Пробовали перенести на монго — он не осиливал подсчёт по формуле за разумное время. На мускуле после оптимизаций выходило до пол-секунды(зависимо от объёма данных после грубо уточняющей выборки).
                    +1
                    Немного непонятно определение «Хранить вес не возможно, так как параметры очень часто изменяются.». На сколько часто? Раз в секунду? Или же они зависят от какого-то третьего значения, которого нет в этой таблице, который динамически меняется?
                      0
                      Зависимо от нагрузки — от нескольких раз в минуту, до раза в час. И да, они зависят от внешних параметров(текущее время и дата).
                      0
                      Извините, читал по-дагонали… PostGIS сюда не подойдет?
                      Насколько я помню в PostgreSQL много очень процедур для вычисления пространственных данных, насчет MySQL не скажу.
                  0
                  Поддержу, Postgresql и Postgis — сильно проще. Поиск ближайших из коробки, кстати, не так давно релизнули.
                  В России этими проблемами занимается Олег Бартунов, вот его презентация на эту тему.
                +4
                Формула называется по-русски «формула гаверсинусов».
                  +7
                  > формуле Хаверсина
                  Оу, с заглавной буквы — вИликий учОный Хаверсин (-:
                  Haversine == haversed sine / half the version, en.wikipedia.org/wiki/Versine
                  The versine or versed sine, versin(θ), is a trigonometric function equal to 1 − cos(θ) and 2sin2(½θ). It appeared in some of the earliest trigonometric tables and was once widespread, but it is now little-used. There are several related functions, most notably the haversine, half the versine, known in the haversine formula of navigation.
                    +1
                    > half the version
                    half the versine
                    (damn you autocorrect).
                    0
                    Я так понимаю, число 6367 в километрах взяли как среднее значение радиуса земли между 6356.78 к полюсам и 6378.14 к экватору?
                      0
                      Вставлю свои три копейки. Описанное выше будет работать пока расстояние будет довольно небольшим, то есть до тех пор покуда ваш радиус не вылетит за пределы -180/180 или -90/90. А дальше уже начинаются проблемы. По крайней мере я так и не отыскал нормального решения к сожалению. Может быть кто-то подскажет в какую сторону смотреть?
                        0
                        Простите, а на какой планете у вас угловые координаты вылетают за пределы [-180, 180], [-90, 90]?
                          0
                          планета Земля. Если есть желание сделать выборку в большом радиусе (5000 км) и вы находитесь скажем на Аляске, а другой человек где-то в Сибири(стоит отметить что радиус должен быть достаточно большим что бы вышеназванное было истиной), то тот человек что находится в Сибири обнаружен не будет(если искать относительно человека в Аляске), ибо вы упретесь в -180. И что дальше делать мне в общем-то непонятно. К сожалению, как я уже сказал выше пока что решения не нашел и буду рад предложениям. Понимаю, что случай в принципе довольно искуственный, но все же имеет право на существование.

                          P.S. Можно конечно брутфорсить(считать длинны в реалтайме), но при количестве координат, выходящим за пределы десятков тысяч позиций это станет неприлично долгим.
                        +3
                        Когда в свое время у нас стояла такая же задача, мы решили ее путем хранения в базе кроме широты и долготы еще дополнительно координат в метрах переведенных из градусов по проекции меркатор.
                        Таким образом мы рассчитывали метрическое окно попадания в коде, а запрос к базе уже был просто на попадание в диапазон.

                        P.S. Формула
                          +1
                          Это же вроде будет работать только для представления Земли как идеальной сферы, а не, скажем, для WGS84?
                            +1
                            Вроде есть для этого? ST_Distance_Sphere более быстрая и ST_Distance_Spheroid более точная
                              0
                              Вник в суть. И подумал, может если есть возможно перенести эти расчеты на клиентскую сторону?

                              var R = 6371; // km
                              var dLat = (lat2-lat1).toRad();
                              var dLon = (lon2-lon1).toRad();
                              var lat1 = lat1.toRad();
                              var lat2 = lat2.toRad();
                              
                              var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
                                      Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2); 
                              var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 
                              var d = R * c;
                              
                                0
                                наткнулся еще на класс написанный на js (если надо напишу в лс)
                                  0
                                  Если задача заключается лишь в том что бы подсчитать расстояние между точками, то тут сложности нет. Можно начать отсюда, а дальше по линкам внутри википедии.
                                  Если задача с довольно большой точностью выбрать достаточно большой «квадрат» с нужной стороной(сторона превышает несколько сотен или тысяч километров) — то вот тут задача становится ВЕСЬМА нетривиальной.
                                +2
                                Between по 2 ключам никогда индексируется. Если уж Вы используйте spatial — функции, там есть замечательная MBRContains и SPATIAL KEY.
                                И кстати да, например для вывода на карту быстрее делать выборку с MBRContains с LIMIT и ORDER BY created DESC, а потом уже проверять расстояния. Тогда действительно работает быстро.
                                  0
                                  По-моему, у вас ошибка в формуле, вот здесь:
                                  src_lat - ABS(dst_lat)
                                  

                                  Зачем здесь брать модуль? Если вы возьмете две точки на одном меридиане, но с противоположными широтами, то ваша формула будет считать их за одну точку и вернет ноль в качестве расстояния.

                                  Кроме того, вот тут модуль тоже не нужен:
                                  COS(ABS(dst_lat) * PI()/180)
                                  

                                  поскольку косинус — четная функция, т.е. cos(x) = cos(-x).
                                    0
                                    для большей производительности, я бы реализовал в UDF
                                    поточнее и побыстрее будет

                                    static const double DEG_TO_RAD = 0.017453292519943295769236907684886;
                                    ///  Earth's quatratic mean radius for WGS-84
                                    static const double EARTH_RADIUS_IN_METERS = 6372797.560856;
                                     
                                    /**  Computes the arc, in radian, between two WGS-84 positions.
                                      *
                                      * The result is equal to <code>Distance(from,to)/EARTH_RADIUS_IN_METERS
                                      *     2*asin(sqrt(h(d/EARTH_RADIUS_IN_METERS )))
                                      *
                                      *    - d is the distance in meters between 'from' and 'to' positions.
                                      *    - h is the haversine function: <code>h(x)=sin²(x/2)
                                      *
                                      * The haversine formula gives:
                                      *    h(d/R) = h(from.lat-to.lat)+h(from.lon-to.lon)+cos(from.lat)*cos(to.lat)
                                      *
                                      *  http://en.wikipedia.org/wiki/Law_of_haversines
                                      */

                                    double ArcInRadians(const Position& from, const Position& to) {
                                        double latitudeArc  = (from.lat - to.lat) * DEG_TO_RAD;
                                        double longitudeArc = (from.lon - to.lon) * DEG_TO_RAD;
                                        double latitudeH = sin(latitudeArc * 0.5);
                                        latitudeH *= latitudeH;
                                        double lontitudeH = sin(longitudeArc * 0.5);
                                        lontitudeH *= lontitudeH;
                                        double tmp = cos(from.lat*DEG_TO_RAD) * cos(to.lat*DEG_TO_RAD);
                                        return 2.0 * asin(sqrt(latitudeH + tmp*lontitudeH));
                                    }
                                     
                                    /**  Computes the distance, in meters, between two WGS-84 positions.
                                      *
                                      * The result is equal to <code>EARTH_RADIUS_IN_METERS*ArcInRadians(from,to)</code>
                                      *
                                      *  ArcInRadians
                                      */

                                    double DistanceInMeters(const Position& from, const Position& to) {
                                        return EARTH_RADIUS_IN_METERS*ArcInRadians(from, to);
                                    }
                                    • UFO just landed and posted this here
                                        0
                                        а каков процент ошибок? Ведь это все еще вероятностный поиск…
                                        • UFO just landed and posted this here
                                        0
                                        Кстати хорошим решением может быть Sphinx.
                                        Произвести индексацию таблицы с координатами по id, lat, lon
                                        и искать уже в сфинксе по встроеному в него geodist'у.
                                        Вполне шустро на 1млн. записей и результаты поиска сразу в метрах от точки поиска.
                                          0
                                          Вот тут кажется ошибочка
                                          ...
                                          DROP PROCEDURE IF EXISTS geobox $$
                                          CREATE PROCEDURE geobox (
                                          ...
                                            SET lat_top := src_lat - (dist / 69);
                                          ...
                                            SET lat_bot := src_lat + (dist / 69);
                                          


                                          Latitude увеличивается на север, т.е. TOP > BOTTOM, а у вас наоборот. Или я запутался в чем-то?
                                            0
                                            И вот ещё ошибка:
                                            Написано:
                                            -- dist -> расстояние от центра в километрах
                                            А по факту работает с милями. Ибо 69 — это количество миль в 1 градусе longtitude. Для километров там должна быть цифра 111.2
                                              0

                                              В 2016-м году не изобретайте велосипед, но используйте PostgreSQL 9.5+ с PostGIS 2.2+, там это всё есть и работает шустро (если создать индекс).


                                              Вот, например, запрос из одного моего тестового задания, ищущий 3 ближайших точки из таблицы (я тестировал на 100 000) и его план запроса:


                                              SELECT
                                                ST_X(vehicles.position::geometry) AS longitude,
                                                ST_Y(vehicles.position::geometry) AS latitude,
                                                vehicles.position <-> 'POINT(37.058515 55.923175)'::geography AS distance_sphere,
                                                ST_Distance(points.position, 'POINT(37.058515 55.923175)'::geography) AS distance_spheroid
                                              FROM vehicles
                                              WHERE available = true
                                              ORDER BY vehicles.position <-> 'POINT(37.058515 55.923175)'::geography ASC
                                              LIMIT 3;

                                              Limit  (cost=0.28..1.66 rows=3 width=56) (actual time=15.580..15.658 rows=3 loops=1)
                                                ->  Index Scan using vehicles_position_idx on vehicles  (cost=0.28..22812.28 rows=49500 width=56) (actual time=15.578..15.655 rows=3 loops=1)
                                                      Order By: ("position" <-> '0101000020E6100000A8A9656B7D8742400EBE30992AF64B40'::geography)
                                                      Filter: available
                                                      Rows Removed by Filter: 12
                                              Planning time: 0.215 ms
                                              Execution time: 15.751 ms

                                              Я как раз использовал формулу гаверсинусов для расчёта расстояния и получаемое значение очень близко к значению, возвращаемому оператором geography <-> geography, который же как раз умеет ускоряться за счёт GiST-индексов. Однако функция ST_Distance считает дистанцию более точно, используя не сферу, а более приближенный к форме земли сфероид.


                                              Если интересно, то можете посмотреть:



                                              А так же можно поставить себе PostGIS, например, из следующего Docker-образа: https://hub.docker.com/r/mdillon/postgis/ и поиграться запросами следующего вида:


                                              SELECT
                                                points.position <-> 'POINT(37.058515 55.923175)'::geography                  AS "a <-> b",
                                                ST_Distance(points.position, 'POINT(37.058515 55.923175)'::geography)        AS "ST_Distance(a, b)",
                                                ST_Distance(points.position, 'POINT(37.058515 55.923175)'::geography, true)  AS "ST_Distance(a, b,  true)",
                                                ST_Distance(points.position, 'POINT(37.058515 55.923175)'::geography, false) AS "ST_Distance(a, b, false)"
                                              FROM (
                                                SELECT
                                                  ST_SetSRID(ST_MakePoint(37.61778 + (n*random() - 5000.00)/50000.00, 55.75583 + (n*random() - 5000.00)/50000.00),4326) AS position
                                                FROM generate_series(1,3) As n
                                              ) points;

                                              Only users with full accounts can post comments. Log in, please.