Как стать автором
Поиск
Написать публикацию
Обновить

Комментарии 13

Спасибо за комментарий. Spatial index не пробовал использовать, но мне очень интересно сделать замеры, этот индекс кажется перспективным для данной задачи.

PS Первоначально простое решение удовлетворяло, именно поэтому проект развивался по шагам описанным выше.

Рейтинг скрыт

Настоятельно рекомендую попробовать... мне как-то они помогли даже в задаче, далёкой от работы с картой, но требовавшей R-tree index для оптимизации - и здОрово, надо сказать, помогло.

Только учитывайте, что далеко не каждая spatial-функция умеет использовать такой индекс. И что самое противное, нет нигде вменяемого полного списка функций, которые используют (или не используют) такой индекс, всё приходится определять опытным путём. По счастью, в наборе геофункций есть сплошь и рядом полу-дубликаты, как правило, одну и ту же проверку можно выполнить несколькими способами.

Рейтинг скрыт

Все станет намного проще, если хоть чуть-чуть погрузиться в предметную область.

Тогда откроется сокровенное знание что геодезические координаты можно перевести в прямоугольную проекцию на плоскость. Например, проекцию Меркатора.

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

Рейтинг скрыт

Разве при таком преобразовании сохраняется расстояние между любыми двумя произвольными точками?

Рейтинг скрыт

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

Есть и другие проекции - Гаусса-Крюгера, UTM... Но там сложнее за счет зонирования.

Меркатор для данной задачи будет оптимальным решением.

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

Достаточно написать UDF и повесить триггер на добавление/изменение записи где будут автоматом считаться поля с проекцией.

Рейтинг скрыт

Если интересны подробности:

/*
 * 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 и хранить там прямоугольные координаты наравне с геодезическими.

Если учесть что

R^{2} = (X_{1} - X_{2})^{2} + (Y_{1} - Y_{2})^{2}

Что экивалентно

R^{2} = X_{1}^{2} - 2*X_{1}*X_{2} + X_{2}^{2} + Y_{1}^{2} - 2*Y_{1}*Y_{2} + Y_{2}^{2}

Или

R^{2} = (X_{1}^{2} + Y_{1}^{2}) - 2*(X_{1}*X_{2} + Y_{1}*Y_{2}) + X_{2}^{2} + Y_{2}^{2}

То можно сразу еще хранить сумму квадратов прямоугольных координат точки (что чуть-чуть ускорит расчет - первая сумма в скобках уже предрасчитана).

Ну и не связываться с квадратными корнями, сразу сравнивать квадрат расстояния.

Рейтинг скрыт

Ладно, закажу пиццу на станцию Амундсен-Скотт) Посмотрю, какой "квадрат" получится посчитать у не-плоскоземельщика)

Мне как-то давно пришлось с sqlite извращаться. Из коробки запроса по земным координатам не было. Поэтому мы в базе хранили не широту и долготу, а пересчитывали 2 полярные координаты в 3 прямоугольные. Запрашивали у базы, что внутри куба со стороной 2R, потом уточняли евклидово расстояние.

Смысл в том, чтобы в базу заранее внести какую-то избыточность, которая сможет ускорить поиск. 3 координаты вместо 2х, или вот про spatial index тут прочитал, или континент/страну/область указать.

Рейтинг скрыт

А при чем тут это? Если вы работаете с ближайшей окрестностью - Меркатор не даст существенной для вашей задачи погрешности. И вполне себе используется во многих приложениях

Для более дальних - надо Гаусса-Крюгера или 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 в помощь. Там все это реализовано на высоком идейно-художественном уровне.

Рейтинг скрыт

Да, без велосипедов редко обходимся)

При пересчете, надеюсь, параметры геоида учитывали? Земля, вообще-то, не шарик...

Не, на это забили. В той задаче нас не волновала погрешность измерения расстояний в 0,3%. Более того, в БД пихали XYZ на единичной сфере, просто синусы-косинусы сферических координат. А средний радиус земли уже учитывали в пред- и пост-обработке запросов.

Только это все медленно работает. Ну или синусы-косинусы-тангенсы-арктангенсы на таблицы переводить.

Так смысл изысканий автора статьи как раз в том и был, чтобы ускорить запрос к БД. За счёт упрощения запроса. UDF не будет работать быстрее, чем ST_Distance_Sphere у автора.

Есть и другие проекции - Гаусса-Крюгера, UTM... Но там сложнее за счет зонирования.

Меркатор для данной задачи будет оптимальным решением.

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

А вот зонирование вполне может ускорить предварительную фильтрацию. Из общепринятых как-нибудь доберусь до Plus Code, если задача будет. Но принцип тот-же - в БД к точкам добавляется избыточная информация, ускоряющая фильтрацию.

Рейтинг скрыт

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

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

Речь идет о том, чтобы выдать список точек в небольшой (там было упоминание о том, что наиболее частым является радиус 0.5км) окрестности от заданной точки. На таких расстояниях координаты в проекции Меркатора дают крайне малую погрешность. Пренебрежимо малую для данной задачи.

Преимущество в данном случае в том, что количество тяжелых вычислений (тригонометрические функции) минимизируется - один раз считаем для точек в БД (при добавлении или изменении записи) и один раз при переводе заданной точки в прямоугольные координаты. Дальше идет массовая операция расчета расстояний между заданной точкой и точками в БД - там уже используются более легкие арифметические операции умножения, сложения и вычитания.

Вообще, тут просится UDF для перевода координат в проекцию, триггер на добавление/изменение записи и UDTF для выдачи списка точек по заданным координатам и радиусу.

А вот зонирование вполне может ускорить предварительную фильтрацию.

Нет. Зонирование придумано для уменьшения искажений, в т.ч. и расстояний при работе с проекциями. Очевидно, что любая карта есть прямоугольная проекция геоида на плоскость. Отсюда и переход к прямоугольным координатам.

Да, проекция Меркатора даст достаточно ощутимые искажения на больших расстояниях (причем, искажения идут исключительно по долготе, обусловлены уменьшением радиуса секущей окружности по мере удаления от экватора). Поэтому делят на зоны шириной 6 градусов и все прямоугольные координаты считаются внутри зоны (т.е. расстояние не от нулевого мередиана, а от начала зоны т.е. появляется третий параметр "центральный мередиан зоны").

При этом для двух точке, лежащих в разных зонах расстояние может быть меньше чем для двух точек лежащих в одной зоне (например - одна точка рядом с концом предыдущей зоны, вторая - в самом начале новой зоны - они окажутся рядом). Но возрастает сложность расчета для точек в разных зонах - это будет сумма расстояний от точки до конца зоны и от начала зоны до второй точки).

Для работы с проекциями и датумами есть очень серьезная библиотека - PROJ Но требует погружения во все это дело. Я когда-то для своих целей делал ограниченный набор функций (мне нужна была только конвертация WGS84-СК42 и конвертации WGS84 в Меркатор или UTM - поперечный Меркатор и СК42 - проекция Гаусса-Крюгера для работы с картами ГШ). Потому нужда отпала, код я отдал другому человеку, который выложил все это на гитхаб. Там же есть функции для расчета номеров тайлов (по координатам точки и масштабному уровню) для запроса их с серверов OSM. Это если не хочется самому заморачиваться с рендерингом - берете готовые картинки (тайл - 256х256 пикселей), склеиваете и выводите на экран.

Рейтинг скрыт

Более того, в БД пихали XYZ на единичной сфере, просто синусы-косинусы сферических координат. А средний радиус земли уже учитывали в пред- и пост-обработке запросов.

Гм... у геоданных есть параметр SRID вообще-то..

Рейтинг скрыт

Не знаю, как работает spatial index в mysql, но, если хотите правильный метод для выборки пространственно локализованных данных, то можно написать какой-нибудь вид R деревьев, например, https://en.m.wikipedia.org/wiki/Priority_R-tree

Рейтинг скрыт
Зарегистрируйтесь на Хабре, чтобы оставить комментарий

Публикации

Истории

Ваш аккаунт

Разделы

Информация

Услуги

Техническая поддержка
© 2006–2024, Habr