Индексация по широте/долготе

Проблема

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

Вы, возможно, пробовали

• INDEX(широта), INDEX(долгота) — но оптимизатор использовал только один
• INDEX(широта, долгота) — но ему все равно приходилось работать слишком тяжело
• Иногда вы приходили к полному сканированию таблицы — Ужас.

WHERE SQRT(...)< ... — Нет возможности использовать какой-либо индекс.

WHERE широта BETWEEN ... AND долгота BETWEEN... — Этот способ имеет некоторые шансы на использование таких индексов.

Цель состоит в том, чтобы просматривать только записи, «близкие» к целевой широте/долготе.

Решение — сначала принципы

Разбиение (PARTITION) в MariaDB и MySQL позволяют создавать способы создания двух кластеризованных индексов. Таким образом, если мы можем разделить (разбить) земной шар по одной оси и использовать обычную индексацию по другой, возможно, мы получим нечто, приближающееся к двумерному индексу. Этот двумерный подход значительно снижает количество обращений к диску по сравнению с одномерными подходами, тем самым ускоряя запросы «найти ближайшее».

Это работает. Не идеально, но лучше, чем альтернативы.

На чем следует производить разбиение (PARTITION)? Похоже, что широта или долгота могли бы подойти. Обратите внимание, что долгота изменяется по ширине — от 69 миль (111 км) на экваторе до 0 на полюсах. Таким образом, широта кажется лучшим выбором.

Сколько разбиений (PARTITION)? Это не очень важно. Некоторые соображения:

• 90 разбиений — 2 градуса каждое. (Я не люблю таблицы с слишком многими разбиениями; 90 кажется достаточным.)
• 50-100 — равномерно заполненные. (Это требует кода. Для 2,7 млн. названий мест 85 разбиений варьировались от 0,5 градусов до очень широких разбиений на полюсах.)
• Не должно быть более 100 разбиений, так как есть недостатки в реализации разбиения.

Как разбить (PARTITION)? Ну, MariaDB и MySQL очень привередливы. Поэтому FLOAT/DOUBLE исключаются. DECIMAL также исключается. Таким образом, мы застряли с каким-то хитрым способом. По существу, нам нужно преобразовать широту/долготу в некоторое значение типа INT и использовать разбиение (PARTITION) по диапазону.

Варианты представления

Чтобы получить тип данных, который может использоваться в разбиении (PARTITION), необходимо «масштабировать» широту и долготу. (Рассмотрим только значения *INT; другие типы данных включены для сравнения)

   Datatype           Bytes       resolution
   ------------------ -----  --------------------------------
   Deg*100 (SMALLINT)     4  1570 m    1.0 mi  Cities
   DECIMAL(4,2)/(5,2)     5  1570 m    1.0 mi  Cities
   SMALLINT scaled        4   682 m    0.4 mi  Cities
   Deg*10000 (MEDIUMINT)  6    16 m     52 ft  Houses/Businesses
   DECIMAL(6,4)/(7,4)     7    16 m     52 ft  Houses/Businesses
   MEDIUMINT scaled       6   2.7 m    8.8 ft
   FLOAT                  8   1.7 m    5.6 ft
   DECIMAL(8,6)/(9,6)     9    16cm    1/2 ft  Friends in a mall
   Deg*10000000 (INT)     8    16mm    5/8 in  Marbles
   DOUBLE                16   3.5nm     ...    Fleas on a dog

(Отсортировано по разрешению)

Что это означает…

Deg*100 (SMALLINT) — вы берете широту/долготу, умножаете на 100, округляете и сохраняете в SMALLINT. Это займет 2 байта для каждой координаты, в общей сложности 4 байта. Два объекта могут быть разделены на 1570 метров, но будут записаны как имеющие одинаковую широту и долготу.

DECIMAL(4,2) для широты и DECIMAL(5,2) для долготы займут 2 + 3 байта и не будут иметь лучшего разрешения, чем Deg*100.

Масштабированный SMALLINT — преобразовать широту в SMALLINT SIGNED, выполнив (градусы / 90 * 32767) и округление; долготу — (градусы / 180 * 32767).

FLOAT имеет 24 значимых бита; DOUBLE — 53. (Они не работают с разбиением (PARTITION), но включены для полноты. Часто люди используют DOUBLE, не понимая, насколько это избыточно и сколько места это занимает.)

Конечно, вы можете использовать DEG*1000 и другие «промежуточные» варианты, но преимуществ нет. DEG*1000 занимает столько же места, сколько DEG*10000, но имеет меньшее разрешение.

Итак, просмотрите список, чтобы узнать, какое разрешение вам нужно, а затем выберите кодирование, которое вам удобно. Однако, поскольку мы собираемся использовать широту в качестве «ключа разбиения (partition key)», она должна быть ограничена одним из типов INT. Для примера кода я буду использовать Deg*10000 (MEDIUMINT).

GCDist — вычисление «расстояния по большой окружности»

GCDist — вспомогательная функция, которая правильно вычисляет расстояние между двумя точками на земном шаре.

Код был протестирован и показал около 20 микросекунд на вызов на ПК 2011 года выпуска. Если вам нужно проверить миллион точек, это займёт 20 секунд — слишком долго для веб-приложения. Таким образом, одна из целей процедуры, которая его использует, заключается в минимизации использования этой функции. С представленным кодом функция должна вызываться всего несколько десятков или сотен раз, за исключением патологических случаев.

Конечно, вы могли бы использовать формулу Пифагора. И это работало бы для большинства приложений. Но это не потребует дополнительных усилий для вычисления по большой окружности. Кроме того, формула по большой окружности работает через полюс и через дату линию. А функция Пифагора не намного быстрее.

Для повышения эффективности GCDist понимает выбранное вами масштабирование и имеет эти данные в коде. Я выбираю «Deg*10000», поэтому функция ожидает 350000 для представления 35 градусов. Если вы выберите другое масштабирование, вам нужно будет изменить код.

GCDist() принимает 4 масштабированные значения DOUBLE — широта1, долгота1, широта2, долгота2 — и возвращает масштабированное число «градусов», представляющее расстояние.

В таблице с вариантами представления указано разрешение в 52 фута для Deg*10000 и DECIMAL(x,4). Вот как это было вычислено: для измерения диагонали между широтой/долготой (0,0) и (0.0001,0.0001) (одна «единица в последнем разряде»): GCDist(0,0,1,1) * 69.172 / 10000 * 5280 = 51.65, где

• 69.172 мили/градус широты
• 10000 единиц на градус для выбранного масштабирования
• 5280 футов/милю.

(Нет, эта функция не учитывает то, что Земля является сплюснутым сфероидом и т. д.)

Требуемая структура таблицы

Будет одна таблица (плюс нормализованные таблицы по мере необходимости). Одна таблица должна быть разбита (partitioned) и индексирована, как указано ниже.

Поля и индексы

• PARTITION BY RANGE(широта)
• широта — масштабированная широта (см. выше)
• долгота — масштабированная долгота
• PRIMARY KEY(долгота, широта, ...) — долгота должна быть первой; что-то должно быть добавлено, чтобы сделать её уникальной
• id — (необязательно) вам может понадобиться идентифицировать строки для ваших целей; AUTO_INCREMENT, если хотите
• INDEX(id) — если `id` — AUTO_INCREMENT, то этот обычный индекс (не UNIQUE, не PRIMARY KEY) необходим
• ENGINE=InnoDB — поэтому PRIMARY KEY будет «кластеризованным»
• Другие индексы — поддерживайте их на минимальном уровне (это общее правило производительности для больших таблиц)

Для большей части этого обсуждения предполагается, что широта имеет тип MEDIUMINT — масштабируется от -90 до +90, умножая на 10000. Аналогично для долготы и от -180 до +180.

PRIMARY KEY должен

• начинаться с `долгота`, так как алгоритм нуждается в «кластеризации», которую предоставит InnoDB, и
• включать `широту` где-то, так как это ключ разбиения (partition key), и
• содержать что-то, чтобы сделать ключ уникальным (долгота+широта вряд ли будет достаточно).

Процедура FindNearest будет выполнять несколько SELECT, примерно так:

      WHERE lat    BETWEEN @my_lat - @dlat
                       AND @my_lat + @dlat   -- PARTITION Pruning and bounding box
        AND lon    BETWEEN @my_lon - @dlon
                       AND @my_lon + @dlon   -- first part of PK
        AND condition                        -- filter out non-pizza parlors

Планировщик запросов

• Произведёт «обрезку» разбиений (PARTITION) на основе широты; затем
• Внутри разбиения (PARTITION), использует долготу для выполнения сканирования по диапазону; затем
• Использует «условие» для фильтрации строк, которые вам нужны, плюс повторная проверка широты. Эта конструкция приводит к очень небольшому количеству блоков диска, которые нужно прочитать, что является основной целью конструкции.

Обратите внимание, что это даже не вызывает GCDist. Это происходит на последнем проходе, когда используется ORDER BY и LIMIT.

В хранимой процедуре хранимой процедуре есть цикл. Будут выполнены как минимум два SELECT, но при правильной настройке обычно не более 6 SELECT. Из-за поиска по PRIMARY KEY каждый SELECT затрагивает только один блок таблицы, иногда больше. Подсчёт числа затронутых блоков — грубый, но эффективный способ сравнения производительности различных конструкций. Для сравнения полное сканирование таблицы, вероятно, затронет тысячи блоков. Простой индекс INDEX(широта), вероятно, приводит к затронутым сотням блоков.

Фильтр… В аргументе хранимой процедуры FindNearest содержится булево выражение («условие») для предложения WHERE. Если вам не нужен никакой фильтр, передайте «1». Для предотвращения «SQL-инъекции» не позволяйте веб-пользователям вводить произвольные выражения; вместо этого постройте «условие» из вводимых ими данных, тем самым гарантируя его безопасность.

Алгоритм

Алгоритм реализован в хранимой процедуре из-за его сложности.

• Вы предоставляете начальную ширину «квадрата» и количество элементов для поиска.
• Он строит «квадрат» вокруг вашей точки.
• Выполняется SELECT, чтобы увидеть, сколько элементов находится в квадрате.
• Цикл, удваивая ширину квадрата, пока не будет найдено достаточно элементов.
• Теперь выполняется «последний» SELECT, чтобы получить точные расстояния, отсортировать их (ORDER BY) и ограничить их желаемым количеством.
• Если пересекается полюс или линия изменения даты, используется более сложный SELECT.

Следующий раздел («Производительность») должен сделать это немного яснее, так как он продемонстрирует несколько примеров.

Производительность

Из-за всех вариаций сложно получить осмысленный эталон. Поэтому вместо этого я представлю некоторые соображения.

Каждый SELECT ограничен «квадратом», определяемым диапазоном широты и диапазоном долготы. (См. предложение WHERE, упомянутое выше, или в примере кода ниже.) Из-за того, как искривляются линии долготы, диапазон долготы «квадрата» будет больше, чем диапазон широты. Предположим, разбиение по широте составляет 3 градуса в области, где вы ищете. Это более 200 миль (более 300 км), поэтому диапазон широты, скорее всего, будет меньше, чем ширина разбиения. Тем не менее, если вы достигаете края полосы широты, квадрат может охватывать два разбиения. После обрезки разбиений до одного (иногда больше) разбиение, запрос затем ограничен диапазоном долготы. (Помните, что PRIMARY KEY начинается с `долгота`.) Если блок данных InnoDB содержит 100 строк (полезное эмпирическое правило), запрос затрагивает один (или несколько) блок. Если квадрат охватывает два (или более) разбиения, та же логика применяется к каждому разбиению.

Таким образом, сканирование квадрата потребует не более одного блока; в редких случаях — не более нескольких блоков. Количество блоков в основном не зависит от размера набора данных.

Основной случай использования этого алгоритма — когда данные значительно больше, чем помещаются в кэш (buffer_pool). Поэтому основная цель — минимизировать количество обращений к диску.

Теперь давайте рассмотрим некоторые граничные случаи и докажем, что количество блоков все же лучше (обычно), чем при использовании традиционных методов индексации.

Что если вы ищете Starbucks в густом городе? Будет десятки, может быть сотни на квадратную милю. Если вы начнете с предположения в 100 милях, SELECTs будут затрагивать множество блоков — неэффективно. В этом случае «начальное расстояние» должно быть небольшим, скажем, 2 мили. Предположим, ваше приложение хочет 10 ближайших магазинов. В этом примере вы, вероятно, найдете более 10 Starbucks в пределах 2 миль в одном блоке InnoDB в одном разбиении. Несмотря на то, что есть второй SELECT для завершения запроса, он будет затрагивать тот же блок. Итого: одно обращение к блоку — дёшево.

Предположим, вы начинаете с квадрата со стороной 5 миль. Поскольку в некоторых густонаселенных городах мира насчитывается более 200 кофеен Starbucks в радиусе 5 миль, это может означать 300 в нашем «квадрате». Это соответствует примерно 4 блокам диска и умеренному количеству ЦП для обработки 300 записей. Все еще неплохо.

Теперь представьте, что вы находитесь на океанском лайнере где-то в Тихом океане. На борту есть одна кофейня Starbucks, но вы ищете 10 ближайших. Если вы снова начнете с 2 миль, для поиска 10 объектов потребуется несколько итераций. Однако давайте пройдемся по этому процессу. Первый запрос попадет в один раздел (возможно, в два) и найдет всего один результат. Второй запрос удваивает ширину квадрата; 4 мили по-прежнему дадут один результат — тот же результат в том же блоке, который теперь кэширован, поэтому мы не будем считать это второй операцией ввода-вывода на диск. В конечном итоге квадрат будет достаточно широким, чтобы охватывать несколько разделов. Каждый дополнительный раздел будет означать один новый доступ к диску для обнаружения отсутствия объектов в квадрате. В конечном итоге квадрат достигнет Чили, Гавайев или Фиджи и найдет несколько других объектов, возможно, достаточно, чтобы остановить итерацию. Поскольку основным критерием определения количества обращений к диску является количество разделов, которые были обработаны, мы не хотим делить мир на слишком много разделов. Если, скажем, есть 40 разделов, то я только что описал случай, в котором может быть 20 обращений к диску.

Разделы в 2 градуса могут подойти для глобальной таблицы магазинов или ресторанов. Исходное расстояние в 5 миль может подойти для фильтрации кофеен Starbucks. 20 миль могут быть лучше для универмага.

Теперь давайте обсудим последний SELECT-запрос, в котором квадрат расширяется на SQRT(2), а для точного упорядочения N результатов используется формула большого круга. SQRT(2) используется в том случае, если N элементов находятся в углах «квадрата». Расширение квадрата на эту величину позволяет нам поймать любые другие объекты, которые находились сразу за пределами старого квадрата.

Прежде всего, обратите внимание, что этот последний SELECT-запрос обращается к тем же блокам(ам), к которым обращался цикл, а также, возможно, к некоторым другим блокам. Сложно предсказать, сколько дополнительных блоков может быть обработано. Вот патологический случай. Вы находитесь посреди пустыни; квадрат расширяется и расширяется. В конечном итоге он находит N объектов. Вне конечного квадрата, полученного в результате итерации, находится большой город. Теперь срабатывает последний SELECT-запрос, и он включает много объектов в этом большом городе. «Много объектов» —> много блоков —> много обращений к диску.

Обсуждение справочного кода

Вот суть хранимой процедуры FindNearest().

• Сделайте предположение о том, как близко к «мне» нужно искать.
• Определите, сколько элементов находится в «квадрате» вокруг меня после фильтрации.
• Если элементов недостаточно, повторите, удвоив ширину квадрата.
• После нахождения достаточного количества элементов или отказа от поиска, поскольку мы ищем «слишком далеко», выполните последний проход для получения всех данных, упорядоченных и ограниченных.

Обратите внимание, что цикл использует только «квадраты» диапазонов широты/долготы. Это грубо, но работает хорошо с разделением и индексацией и избегает вызова GCDist (до последнего шага). В примере кода в качестве начального значения выбрано 15 миль. Изменение этого значения повлияет на производительность процедуры, но влияние будет зависеть от конкретных случаев использования. Грубый способ установить радиус — предположить, что это позволит найти требуемое LIMIT приблизительно в половине случаев. (Это значение жестко закодировано в процедуре.)

Параметры, передаваемые в FindNearest():

• Ваша широта — -90..90 (не масштабируется — см. жестко закодированное преобразование в процедуре)
• Ваша долгота — -180..180 (не масштабируется)
• Начальное расстояние — (мили или км) — см. обсуждение ниже
• Максимальное расстояние — в милях или км — см. жестко закодированное преобразование в процедуре
• Ограничение — максимальное количество возвращаемых элементов
• Условие — что-то, что нужно добавить после «И» (более подробное обсуждение выше)

Функция найдет ближайшие элементы, не превышая Limit, которые удовлетворяют условию. Но она прекратит поиск на максимальном расстоянии. (Если вы находитесь в Южном полюсе, зачем тратить время на поиск десятой пиццерии очень далеко?)

Из-за «масштабирования», «жесткой кодировки», «условия», имени таблицы и т. д. эта процедура не является полностью универсальной; код должен быть изменен для каждого приложения. Да, я мог бы спроектировать его таким образом, чтобы передавать все эти данные. Но это было бы беспорядочно.

Значение «_start_dist» позволяет контролировать производительность. Слишком маленькое значение приводит к дополнительным итерациям, а слишком большое значение приводит к проверке большего числа строк. Если вы решите настроить хранимую процедуру, выполните следующее. «SELECT @iterations» после вызова хранимой процедуры для ряда типичных значений. Если значение обычно равно 1, уменьшите _start_dist. Если оно обычно равно 2 или более, увеличьте его.

Время выполнения: менее 10 мс для «типичного» использования; любой размер набора данных. Медленнее для патологических случаев (маленькое минимальное расстояние, большое максимальное расстояние, пересечение датской линии, плохая фильтрация, холодный кэш и т. д.)

Крайние случаи:

• Использование расстояния по большой окружности, а не по теореме Пифагора, позволяет получить «правильные» расстояния даже вблизи полюсов.
• Полюса — даже если «ближайший» объект находится почти на 360 градусов (долгота) от вас, его можно найти.
• Линия перемены дат — есть небольшой, «содержимый», фрагмент кода для пересечения линии перемены дат. Пример: вы находитесь в +179 градусах долготы, а ближайший объект находится в -179 градусах.

Процедура возвращает один результирующий набор, SELECT *, distance.

• Возвращаются только строки, которые соответствуют вашему условию и находятся в пределах максимального расстояния.
• Возвращается не более Limit строк.
• Строки будут отсортированы по близости, начиная с наименьшего расстояния.
• «Расстояние» будет указано в милях или километрах (в зависимости от жестко закодированной константы в хранимой процедуре).

Справочный код, предполагая deg*10000 и «мили»

Эта версия основана на масштабировании «Deg*10000 (MEDIUMINT»).

DELIMITER //

drop function if exists GCDist //
CREATE FUNCTION GCDist (
        _lat1 DOUBLE,  -- Scaled Degrees north for one point
        _lon1 DOUBLE,  -- Scaled Degrees west for one point
        _lat2 DOUBLE,  -- other point
        _lon2 DOUBLE
    ) RETURNS DOUBLE
    DETERMINISTIC
    CONTAINS SQL  -- SQL but does not read or write
    SQL SECURITY INVOKER  -- No special privileges granted
-- Input is a pair of latitudes/longitudes multiplied by 10000.
--    For example, the south pole has latitude -900000.
-- Multiply output by .0069172 to get miles between the two points
--    or by .0111325 to get kilometers
BEGIN
    -- Hardcoded constant:
    DECLARE _deg2rad DOUBLE DEFAULT PI()/1800000;  -- For scaled by 1e4 to MEDIUMINT
    DECLARE _rlat1 DOUBLE DEFAULT _deg2rad * _lat1;
    DECLARE _rlat2 DOUBLE DEFAULT _deg2rad * _lat2;
    -- compute as if earth's radius = 1.0
    DECLARE _rlond DOUBLE DEFAULT _deg2rad * (_lon1 - _lon2);
    DECLARE _m     DOUBLE DEFAULT COS(_rlat2);
    DECLARE _x     DOUBLE DEFAULT COS(_rlat1) - _m * COS(_rlond);
    DECLARE _y     DOUBLE DEFAULT               _m * SIN(_rlond);
    DECLARE _z     DOUBLE DEFAULT SIN(_rlat1) - SIN(_rlat2);
    DECLARE _n     DOUBLE DEFAULT SQRT(
                        _x * _x +
                        _y * _y +
                        _z * _z    );
    RETURN  2 * ASIN(_n / 2) / _deg2rad;   -- again--scaled degrees
END;
//
DELIMITER ;

DELIMITER //
-- FindNearest (about my 6th approach)
drop procedure if exists FindNearest6 //
CREATE
PROCEDURE FindNearest (
        IN _my_lat DOUBLE,  -- Latitude of me [-90..90] (not scaled)
        IN _my_lon DOUBLE,  -- Longitude [-180..180]
        IN _START_dist DOUBLE,  -- Starting estimate of how far to search: miles or km
        IN _max_dist DOUBLE,  -- Limit how far to search: miles or km
        IN _limit INT,     -- How many items to try to get
        IN _condition VARCHAR(1111)   -- will be ANDed in a WHERE clause
    )
    DETERMINISTIC
BEGIN
    -- lat and lng are in degrees -90..+90 and -180..+180
    -- All computations done in Latitude degrees.
    -- Thing to tailor
    --   *Locations* -- the table
    --   Scaling of lat, lon; here using *10000 in MEDIUMINT
    --   Table name
    --   miles versus km.

    -- Hardcoded constant:
    DECLARE _deg2rad DOUBLE DEFAULT PI()/1800000;  -- For scaled by 1e4 to MEDIUMINT

    -- Cannot use params in PREPARE, so switch to @variables:
    -- Hardcoded constant:
    SET @my_lat := _my_lat * 10000,
        @my_lon := _my_lon * 10000,
        @deg2dist := 0.0069172,  -- 69.172 for miles; 111.325 for km  *** (mi vs km)
        @start_deg := _start_dist / @deg2dist,  -- Start with this radius first (eg, 15 miles)
        @max_deg := _max_dist / @deg2dist,
        @cutoff := @max_deg / SQRT(2),  -- (slightly pessimistic)
        @dlat := @start_deg,  -- note: must stay positive
        @lon2lat := COS(_deg2rad * @my_lat),
        @iterations := 0;        -- just debugging

    -- Loop through, expanding search
    --   Search a 'square', repeat with bigger square until find enough rows
    --   If the inital probe found _limit rows, then probably the first
    --   iteration here will find the desired data.
    -- Hardcoded table name:
    -- This is the "first SELECT":
    SET @sql = CONCAT(
        "SELECT COUNT(*) INTO @near_ct
            FROM Locations
            WHERE lat    BETWEEN @my_lat - @dlat
                             AND @my_lat + @dlat   -- PARTITION Pruning and bounding box
              AND lon    BETWEEN @my_lon - @dlon
                             AND @my_lon + @dlon   -- first part of PK
              AND ", _condition);
    PREPARE _sql FROM @sql;
    MainLoop: LOOP
        SET @iterations := @iterations + 1;
        -- The main probe: Search a 'square'
        SET @dlon := ABS(@dlat / @lon2lat);  -- good enough for now  -- note: must stay positive
        -- Hardcoded constants:
        SET @dlon := IF(ABS(@my_lat) + @dlat >= 900000, 3600001, @dlon);  -- near a Pole
        EXECUTE _sql;
        IF ( @near_ct >= _limit OR         -- Found enough
             @dlat >= @cutoff ) THEN       -- Give up (too far)
            LEAVE MainLoop;
        END IF;
        -- Expand 'square':
        SET @dlat := LEAST(2 * @dlat, @cutoff);   -- Double the radius to search
    END LOOP MainLoop;
    DEALLOCATE PREPARE _sql;

    -- Out of loop because found _limit items, or going too far.
    -- Expand range by about 1.4 (but not past _max_dist),
    -- then fetch details on nearest 10.

    -- Hardcoded constant:
    SET @dlat := IF( @dlat >= @max_deg OR @dlon >= 1800000,
                @max_deg,
                GCDist(ABS(@my_lat), @my_lon,
                       ABS(@my_lat) - @dlat, @my_lon - @dlon) );
            -- ABS: go toward equator to find farthest corner (also avoids poles)
            -- Dateline: not a problem (see GCDist code)

    -- Reach for longitude line at right angle:
    -- sin(dlon)*cos(lat) = sin(dlat)
    -- Hardcoded constant:
    SET @dlon := IFNULL(ASIN(SIN(_deg2rad * @dlat) /
                             COS(_deg2rad * @my_lat))
                            / _deg2rad -- precise
                        , 3600001);    -- must be too near a pole

    -- This is the "last SELECT":
    -- Hardcoded constants:
    IF (ABS(@my_lon) + @dlon < 1800000 OR    -- Usual case - not crossing dateline
        ABS(@my_lat) + @dlat <  900000) THEN -- crossing pole, so dateline not an issue
        -- Hardcoded table name:
        SET @sql = CONCAT(
            "SELECT *,
                    @deg2dist * GCDist(@my_lat, @my_lon, lat, lon) AS dist
                FROM Locations
                WHERE lat BETWEEN @my_lat - @dlat
                              AND @my_lat + @dlat   -- PARTITION Pruning and bounding box
                  AND lon BETWEEN @my_lon - @dlon
                              AND @my_lon + @dlon   -- first part of PK
                  AND ", _condition, "
                HAVING dist <= ", _max_dist, "
                ORDER BY dist
                LIMIT ", _limit
                        );
    ELSE
        -- Hardcoded constants and table name:
        -- Circle crosses dateline, do two SELECTs, one for each side
        SET @west_lon := IF(@my_lon < 0, @my_lon, @my_lon - 3600000);
        SET @east_lon := @west_lon + 3600000;
        -- One of those will be beyond +/- 180; this gets points beyond the dateline
        SET @sql = CONCAT(
            "( SELECT *,
                    @deg2dist * GCDist(@my_lat, @west_lon, lat, lon) AS dist
                FROM Locations
                WHERE lat BETWEEN @my_lat - @dlat
                              AND @my_lat + @dlat   -- PARTITION Pruning and bounding box
                  AND lon BETWEEN @west_lon - @dlon
                              AND @west_lon + @dlon   -- first part of PK
                  AND ", _condition, "
                HAVING dist <= ", _max_dist, " )
            UNION ALL
            ( SELECT *,
                    @deg2dist * GCDist(@my_lat, @east_lon, lat, lon) AS dist
                FROM Locations
                WHERE lat BETWEEN @my_lat - @dlat
                              AND @my_lat + @dlat   -- PARTITION Pruning and bounding box
                  AND lon BETWEEN @east_lon - @dlon
                              AND @east_lon + @dlon   -- first part of PK
                  AND ", _condition, "
                HAVING dist <= ", _max_dist, " )
            ORDER BY dist
            LIMIT ", _limit
                        );
    END IF;

    PREPARE _sql FROM @sql;
    EXECUTE _sql;
    DEALLOCATE PREPARE _sql;
END;
//
DELIMITER ;
<<code>>

== Sample

Find the 5 cities with non-zero population (out of 3 million) nearest to (+35.15, -90.15). Start with a 10-mile bounding box and give up at 100 miles.

<<code>>
CALL FindNearestLL(35.15, -90.05, 10, 100, 5, 'population > 0');
+---------+--------+---------+---------+--------------+--------------+-------+------------+--------------+---------------------+------------------------+
| id      | lat    | lon     | country | ascii_city   | city         | state | population | @gcd_ct := 0 | dist                | @gcd_ct := @gcd_ct + 1 |
+---------+--------+---------+---------+--------------+--------------+-------+------------+--------------+---------------------+------------------------+
| 3023545 | 351494 | -900489 | us      | memphis      | Memphis      | TN    |     641608 |            0 | 0.07478733189367963 |                      3 |
| 2917711 | 351464 | -901844 | us      | west memphis | West Memphis | AR    |      28065 |            0 |   7.605683607627499 |                      2 |
| 2916457 | 352144 | -901964 | us      | marion       | Marion       | AR    |       9227 |            0 |     9.3994963998986 |                      1 |
| 3020923 | 352044 | -898739 | us      | bartlett     | Bartlett     | TN    |      43264 |            0 |  10.643941157860604 |                      7 |
| 2974644 | 349889 | -900125 | us      | southaven    | Southaven    | MS    |      38578 |            0 |  11.344042217329935 |                      5 |
+---------+--------+---------+---------+--------------+--------------+-------+------------+--------------+---------------------+------------------------+
5 rows in set (0.00 sec)
Query OK, 0 rows affected (0.04 sec)

SELECT COUNT(*) FROM ll_table;
+----------+
| COUNT(*) |
+----------+
|  3173958 |
+----------+
1 row in set (5.04 sec)

FLUSH STATUS;
CALL...
SHOW SESSION STATUS LIKE 'Handler%';

show session status like 'Handler%';
+----------------------------+-------+
| Variable_name              | Value |
+----------------------------+-------+
| Handler_read_first         | 1     |
| Handler_read_key           | 3     |
| Handler_read_next          | 1307  |  -- some index, some tmp, but far less than 3 million.
| Handler_read_rnd           | 5     |
| Handler_read_rnd_next      | 13    |
| Handler_write              | 12    |  -- it needed a tmp
+----------------------------+-------+

Послелоги

Существует алгоритм «Хаверсинуса», который вдвое быстрее функции GCDist здесь. Но у него есть фатальный недостаток — иногда он возвращает NULL для расстояния между точкой и самой собой. (Это происходит из-за вычисления числа, немного большего, чем 1,0, а затем попытки вычислить арккосинус этого числа.)

См. также

• Используемые для тестирования города
• Форумная тема
• Обсуждение на StackOverflow
• Пример
• Z-упорядочение

Рик Джеймс любезно разрешил нам использовать эту статью в базе знаний.

Сайт Рика Джеймса содержит другие полезные советы, пошаговые инструкции, оптимизации и советы по отладке.

Исходный источник: http://mysql.rjweb.org/doc.php/latlng