В настоящее время часто встречается задача - на карте расположено множество точек и необходимо выделить их подмножество (отфильтровать) по критерию расстояния от некоторой “особой” точки. В качестве примеров укажем сайты маркетплейсов, где точки пункты выдачи, в сетевых компаниях это офисы на карте, логистические компании отслеживают передвижение транспортных средств, операторы такси находят на карте ближайшие автомобили и можно перечислить еще много областей применения. При кажущейся простоте постановки задачи её решение требует значительных ресурсов, так как запросы к базе данных, как правило, очень тяжелые. Здесь на примере задачи ‘Пункты выдачи’ с фильтрацией по расстоянию (радиусу) от накопительного узла товаров показывается как можно сократить количество обработанных строк в запросе. Реализация задачи на БД mySQL.
Координаты исходной точки могут быть получены несколькими способами, первый если устройство пользователя оборудовано GPS и пользователь разрешил передачу своих координат, то мы получаем координаты пользователя без затруднений, второй способ менее точный мы получаем координаты пользователя из его IP, для этого используется библиотека maxmind.
В БД mySQL 5.7 появилась встроенная функция ST_Distance_Sphere для расчета (евклидово) расстояния между двумя точками. Данная функция принимает широту и долготу двух точек, а возвращает расстояние между ними в метрах, далее мы производим сравнение с критической величиной радиуса и определяем попадает ли данная точка в круг с этим радиусом.
WHERE (ST_Distance_Sphere(point(longitude, latitude), point(longitude2, latitude2)) ) <= radius
На этом решение задачи можно было бы закончить, но если выполнить данный запрос на миллионе строк, то запрос будет заметно тормозить.
На лицо необходимость предварительной фильтрации. Было решено производить фильтрацию при помощи квадрата, внутрь которого будет вписана окружность критического радиуса. При этом пункты выдачи, которые попадают в квадрат, но не попадают в окружность в последующем будут отфильтрованы как раз при помощи функции ST_Distance_Sphere , то есть львиная доля фильтрации придется именно на квадрат. Итак, как сделать фильтрацию по квадрату более быстрой, чем фильтрация при помощи функции ST_Distance_Sphere.
Вычисляем координаты вершин квадрата.
Сторона квадрата будет равняться диаметру окружности, то есть радиус*2.
Формула вычисления координат вершин квадрата:
Distance - искомая нами дистанция в километрах (радиус).
Экваториальная окружность Земли составляет около 40 075 км (не будем использовать меридианную 40 008).
Окружность 360 градусов.
Для полного понимания данных расчетов необходимо понимать:
Широта и долгота измеряются в градусах:
градусом географической широты является 1/180 часть меридиана (latitude - широта значения от -90 до 90)
градусом географической долготы является 1/360 часть экватора (longitude - долгота значения от -180 до 180)
Небольшое отхождение от темы. Данная формула не работает для тех кто считает что земля плоская))), а таких с каждым днем все больше(((
Итак, мы вычислили координаты вершин квадрата, в который вписан наш круг и можем применить промежуточную фильтрацию + основную. Фильтрация по WHERE between работает быстрее, чем встроенная фильтрация по ST_Distance_Sphere
WHERE longitude between longitudeMin and longitudeMax
and latitude between latitudeMin and latitudeMax
and (ST_Distance_Sphere(point(longitude, latitude),
point(longitude2, latitude2)) ) <= radius // radius в метрахВ качестве оптимизации добавляем индекс на столбцы:longitude, latitude
CREATE INDEX coordinate ON table_name (latitude, longitude);
Данная оптимизация позволила без тормозов обрабатывать 2 млн. записей.
Но было решено на этом не останавливаться и произвести еще одну оптимизацию.
В связи с тем что пользователи чаще всего используют фильтр с радиусом 0-5 км, было решено добавить еще один столбец в таблицу для фильтрации, этот столбец будет целочисленными и его значение будет от -90 000 до 90000.
ALTER TABLE table_name ADD COLUMN latitudeInt INT DEFAULT NULL
Заполнять этот столбцы мы будем по формуле: latitudeInt = Math.round(latitude*1000). Итак, что же мы делаем, по факту имея в БД значение latitude хх.хххххххх мы его приводим к целочисленному виду и округляем до 5 знаков ххххх. При текущей фильтрации с радиусом 0-5 км пяти знаков достаточно для произведения фильтрации, тем самым мы упрощаем процесс фильтрации, тем более фильтрацию с помощью функции ST_Distance_Sphere мы оставляем.
Представим, что координаты исходной точки latitude = 10.111 и longitude = 41.000 по формуле выше вычисляем:
longitude min - 40.98951994951971
longitude max - 41.01048005048029
latitude min - 10.10051684341859
latitude max - 10.12131565814199
По latitude min и latitude max вычисляем искомый диапазон latitudeInt.
getRange(latMin, latMax) {
const latitudeMinInt = Math.floor(latMin * 1000)
const latitudeMaxInt = Math.ceil(latMax * 1000)
const mas = [latitudeMinInt];
for (let index = 1; index < latitudeMaxInt - latitudeMinInt + 1; index++) {
mas.push(latitudeMinInt+index)
}
return mas;
}В нашем случае диапазон будет:
[10102,10103,10104,10105,10106,10107,10108,10109,10110,10111,10112,10113,10114,10115,10116,10117,10118,10119,10120,10121,10122]
В случае если несколько пунктов выдачи попали например, в диапазон 10107 - 10108, то фактически мы их склеили одним значением поля latitudeInt
Теперь фильтры будут представлены следующим образом:
latitudeIntMas = [10102,10103,10104,10105,10106,10107,10108,10109,10110,10111,10112, 10113,10114,10115,10116,10117,10118,10119,10120,10121,10122]
longitudeMin - 40.98951994951971
longitudeMax -41.01048005048029
WHERE latitudeInt in (latitudeIntMas )
AND longitude between longitudeMin and longitudeMax
AND (ST_Distance_Sphere(point(longitude, latitude), point(longitude2, latitude2)) ) <= radius // radius в метрах
Так же добавляем индекс
CREATE INDEX latitudeInt ON table_name (latitudeInt);
Заметим, что данные расчеты не учитывают кривизну земли, но для основного запроса потребителя 0-5 км погрешность не является существенной.
Я постарался максимально подробно рассказать о моем кейсе. Надеюсь этот опыт поможет другим разработчикам, так как данный функционал используется на многих проектах. Если вы решали подобные задачи, пишите в комментариях, будет интересно посмотреть решения других разработчиков.
Комментарии (13)
SpiderEkb
24.10.2023 20:16-1Все станет намного проще, если хоть чуть-чуть погрузиться в предметную область.
Тогда откроется сокровенное знание что геодезические координаты можно перевести в прямоугольную проекцию на плоскость. Например, проекцию Меркатора.
И тогда задача сводится к тривиальной - расстояние между двумя точками с известными прямоугольными координатами на плоскости. И даже волшебная функция не потребуется.
Andy_U
24.10.2023 20:16+1Разве при таком преобразовании сохраняется расстояние между любыми двумя произвольными точками?
SpiderEkb
24.10.2023 20:16На тех расстояниях, которые имеют смысл в данной задаче, погрешность будет пренебрежимо малой.
Есть и другие проекции - Гаусса-Крюгера, UTM... Но там сложнее за счет зонирования.
Меркатор для данной задачи будет оптимальным решением.
Из геодезических координат в меркатор пересчет простой (обратно сложнее - там итеративный процесс), но обратное преобразование не требуется тут.
Достаточно написать UDF и повесить триггер на добавление/изменение записи где будут автоматом считаться поля с проекцией.
SpiderEkb
24.10.2023 20:16+1Если интересны подробности:
/* * Mercator transformation * accounts for the fact that the earth is not a sphere, but a spheroid */ #define D_R (M_PI / 180.0) #define R_D (180.0 / M_PI) #define R_MAJOR 6378137.0 #define R_MINOR 6356752.3142 #define RATIO (R_MINOR/R_MAJOR) #define ECCENT (sqrt(1.0 - (RATIO * RATIO))) #define COM (0.5 * ECCENT) double deg_rad(double ang) { return ang * D_R; } double rad_deg(double ang) { return ang * R_D; } double fmin(double f1, double f2) { return (f1 > f2) ? f2 : f1; } double fmax(double f1, double f2) { return (f1 > f2) ? f1 : f2; } double Lon2MercatorX(double lon) { return R_MAJOR * deg_rad (lon); } double Lat2MercatorY(double lat) { lat = fmin (89.5, fmax (lat, -89.5)); double phi = deg_rad(lat); double sinphi = sin(phi); double con = ECCENT * sinphi; con = pow((1.0 - con) / (1.0 + con), COM); double ts = tan(0.5 * (M_PI * 0.5 - phi)) / con; return 0.0 - R_MAJOR * log(ts); } double MercatorX2Lon(double x) { return rad_deg(x) / R_MAJOR; } double MercatorY2Lat(double y) { double ts = exp(-y / R_MAJOR); double phi = M_PI_2 - 2.0 * atan(ts); double dphi = 1.0; for (int i = 0; (fabs(dphi) > 0.000000001 && i < 15); i++) { double con = ECCENT * sin (phi); dphi = M_PI_2 - 2.0 * atan (ts * pow((1.0 - con) / (1.0 + con), COM)) - phi; phi += dphi; } return rad_deg (phi); }Координаты в проекции Меркатора тут в метрах. Геодезические - в градусах. Расстояния в пределах нескольких километров будут с достаточной для решения задачи точностью.
Достаточно добавить два поля: MerX и MerY и хранить там прямоугольные координаты наравне с геодезическими.
Если учесть что
Что экивалентно
Или
То можно сразу еще хранить сумму квадратов прямоугольных координат точки (что чуть-чуть ускорит расчет - первая сумма в скобках уже предрасчитана).
Ну и не связываться с квадратными корнями, сразу сравнивать квадрат расстояния.
alcotel
24.10.2023 20:16Ладно, закажу пиццу на станцию Амундсен-Скотт) Посмотрю, какой "квадрат" получится посчитать у не-плоскоземельщика)
Мне как-то давно пришлось с sqlite извращаться. Из коробки запроса по земным координатам не было. Поэтому мы в базе хранили не широту и долготу, а пересчитывали 2 полярные координаты в 3 прямоугольные. Запрашивали у базы, что внутри куба со стороной 2R, потом уточняли евклидово расстояние.
Смысл в том, чтобы в базу заранее внести какую-то избыточность, которая сможет ускорить поиск. 3 координаты вместо 2х, или вот про spatial index тут прочитал, или континент/страну/область указать.
SpiderEkb
24.10.2023 20:16А при чем тут это? Если вы работаете с ближайшей окрестностью - Меркатор не даст существенной для вашей задачи погрешности. И вполне себе используется во многих приложениях
Для более дальних - надо Гаусса-Крюгера или UTM (с зонами и то и другое). Например, сетка и метрическая разметка на картах ГШ идет в проекции Г-К (правда, там датум не WGS, а СК-42).
Кстати, для полюсов отдельная тема.
Мне как-то давно пришлось с sqlite извращаться. Из коробки запроса по земным координатам не было.
Написать UDF, я так понимаю, задача непосильная?
static const double Pi = 3.14159265358979; static const double dEarthRadius = 6378.135; static const double dFlat = 0.993291; double DistanceLL(double LatA, double LonA, double LatB, double LonB) { double TanA, TanB, GeoA, GeoB, ColA, ColB; double Sin_ColA, Cos_ColA, Sin_ColB, Cos_ColB; double DLon, CosDelta, DeltaRad; double Cos_DLon, CoLat, Cos_CoLat, DeltaM; // Conversation to radiants LatA = LatA * Pi / 180.0; LonA = LonA * Pi / 180.0; LatB = LatB * Pi / 180.0; LonB = LonB * Pi / 180.0; /* Determning latitutes in co latitudes Point A and Sine & Cosine values */ TanA = tan(LatA) * dFlat; GeoA = atan(TanA); ColA = M_PI_2 - GeoA; Sin_ColA = sin(ColA); Cos_ColA = cos(ColA); /* Determning latitutes in co latitudes Point A and Sine & Cosine values */ TanB = tan(LatB) * dFlat; GeoB = atan(TanB); ColB = M_PI_2 - GeoB; Sin_ColB = sin(ColB); Cos_ColB = cos(ColB); // Determening Distance between A and B DLon = LonB - LonA; Cos_DLon = cos(DLon); CosDelta = Sin_ColA * Sin_ColB * Cos_DLon + Cos_ColA * Cos_ColB; if(CosDelta > 1.0) CosDelta = 1.0; else if(CosDelta < -1.0) CosDelta = -1.0; DeltaRad=acos(CosDelta); // Determening distance in meter CoLat = M_PI_2 - (LatA + LatB) / 2.0; Cos_CoLat = cos(CoLat); DeltaM = DeltaRad * ((1.0/3.0 - Cos_CoLat * Cos_CoLat) * 0.00335278 + 1.0) * dEarthRadius * 1000; return DeltaM; }И никакого пересчета не надо. Подозреваю, что в вашей функции этот алгоритм и реализован.
Только это все медленно работает. Ну или синусы-косинусы-тангенсы-арктангенсы на таблицы переводить.
Поэтому мы в базе хранили не широту и долготу, а пересчитывали 2 полярные координаты в 3 прямоугольные.
Ну велосипедостроение никто не отменял. При пересчете, надеюсь, параметры геоида учитывали? Земля, вообще-то, не шарик...
или континент/страну/область указать
Лучше почитать про проекции UTM и Г-К. Что такое зоны, центральный меридиан зоны и т.п.
Хоти что бы совсем точно-точно - libproj в помощь. Там все это реализовано на высоком идейно-художественном уровне.
alcotel
24.10.2023 20:16Да, без велосипедов редко обходимся)
При пересчете, надеюсь, параметры геоида учитывали? Земля, вообще-то, не шарик...
Не, на это забили. В той задаче нас не волновала погрешность измерения расстояний в 0,3%. Более того, в БД пихали XYZ на единичной сфере, просто синусы-косинусы сферических координат. А средний радиус земли уже учитывали в пред- и пост-обработке запросов.
Только это все медленно работает. Ну или синусы-косинусы-тангенсы-арктангенсы на таблицы переводить.
Так смысл изысканий автора статьи как раз в том и был, чтобы ускорить запрос к БД. За счёт упрощения запроса. UDF не будет работать быстрее, чем ST_Distance_Sphere у автора.
Есть и другие проекции - Гаусса-Крюгера, UTM... Но там сложнее за счет зонирования.
Меркатор для данной задачи будет оптимальным решением.
Не увидел, чем помогают измерению расстояний проекции на плоскость, кроме лёгкого искажения тех-же сферических координат. Точнее дают - дополнительный геморрой на краях карты.
А вот зонирование вполне может ускорить предварительную фильтрацию. Из общепринятых как-нибудь доберусь до Plus Code, если задача будет. Но принцип тот-же - в БД к точкам добавляется избыточная информация, ускоряющая фильтрацию.
SpiderEkb
24.10.2023 20:16Не увидел, чем помогают измерению расстояний проекции на плоскость, кроме лёгкого искажения тех-же сферических координат. Точнее дают - дополнительный геморрой на краях карты.
В данной задаче, насколько я понял, речь не идет о точном расчете расстояний между двумя сильно удаленными точками.
Речь идет о том, чтобы выдать список точек в небольшой (там было упоминание о том, что наиболее частым является радиус 0.5км) окрестности от заданной точки. На таких расстояниях координаты в проекции Меркатора дают крайне малую погрешность. Пренебрежимо малую для данной задачи.
Преимущество в данном случае в том, что количество тяжелых вычислений (тригонометрические функции) минимизируется - один раз считаем для точек в БД (при добавлении или изменении записи) и один раз при переводе заданной точки в прямоугольные координаты. Дальше идет массовая операция расчета расстояний между заданной точкой и точками в БД - там уже используются более легкие арифметические операции умножения, сложения и вычитания.
Вообще, тут просится UDF для перевода координат в проекцию, триггер на добавление/изменение записи и UDTF для выдачи списка точек по заданным координатам и радиусу.
А вот зонирование вполне может ускорить предварительную фильтрацию.
Нет. Зонирование придумано для уменьшения искажений, в т.ч. и расстояний при работе с проекциями. Очевидно, что любая карта есть прямоугольная проекция геоида на плоскость. Отсюда и переход к прямоугольным координатам.
Да, проекция Меркатора даст достаточно ощутимые искажения на больших расстояниях (причем, искажения идут исключительно по долготе, обусловлены уменьшением радиуса секущей окружности по мере удаления от экватора). Поэтому делят на зоны шириной 6 градусов и все прямоугольные координаты считаются внутри зоны (т.е. расстояние не от нулевого мередиана, а от начала зоны т.е. появляется третий параметр "центральный мередиан зоны").
При этом для двух точке, лежащих в разных зонах расстояние может быть меньше чем для двух точек лежащих в одной зоне (например - одна точка рядом с концом предыдущей зоны, вторая - в самом начале новой зоны - они окажутся рядом). Но возрастает сложность расчета для точек в разных зонах - это будет сумма расстояний от точки до конца зоны и от начала зоны до второй точки).
Для работы с проекциями и датумами есть очень серьезная библиотека - PROJ Но требует погружения во все это дело. Я когда-то для своих целей делал ограниченный набор функций (мне нужна была только конвертация WGS84-СК42 и конвертации WGS84 в Меркатор или UTM - поперечный Меркатор и СК42 - проекция Гаусса-Крюгера для работы с картами ГШ). Потому нужда отпала, код я отдал другому человеку, который выложил все это на гитхаб. Там же есть функции для расчета номеров тайлов (по координатам точки и масштабному уровню) для запроса их с серверов OSM. Это если не хочется самому заморачиваться с рендерингом - берете готовые картинки (тайл - 256х256 пикселей), склеиваете и выводите на экран.
Akina
24.10.2023 20:16Более того, в БД пихали XYZ на единичной сфере, просто синусы-косинусы сферических координат. А средний радиус земли уже учитывали в пред- и пост-обработке запросов.
Гм... у геоданных есть параметр SRID вообще-то..
wintsa
24.10.2023 20:16+1Не знаю, как работает spatial index в mysql, но, если хотите правильный метод для выборки пространственно локализованных данных, то можно написать какой-нибудь вид R деревьев, например, https://en.m.wikipedia.org/wiki/Priority_R-tree
Akina
Чем рожать костыль, лучше было бы построить и использовать spatial index.
https://dev.mysql.com/doc/refman/5.7/en/using-spatial-indexes.html
Gotlieb Автор
Спасибо за комментарий. Spatial index не пробовал использовать, но мне очень интересно сделать замеры, этот индекс кажется перспективным для данной задачи.
PS Первоначально простое решение удовлетворяло, именно поэтому проект развивался по шагам описанным выше.
Akina
Настоятельно рекомендую попробовать... мне как-то они помогли даже в задаче, далёкой от работы с картой, но требовавшей R-tree index для оптимизации - и здОрово, надо сказать, помогло.
Только учитывайте, что далеко не каждая spatial-функция умеет использовать такой индекс. И что самое противное, нет нигде вменяемого полного списка функций, которые используют (или не используют) такой индекс, всё приходится определять опытным путём. По счастью, в наборе геофункций есть сплошь и рядом полу-дубликаты, как правило, одну и ту же проверку можно выполнить несколькими способами.