Задачи на сфере: угловая засечка: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Строка 55: Строка 55:


== Пример программной реализации ==
== Пример программной реализации ==
<syntaxhighlight lang="c">
/*
* Решение угловой засечки
*
* Аргументы исходные:
*    pt1 - {широта, долгота} пункта Q1
*    pt2 - {широта, долгота} пункта Q2
*
* Аргументы определяемые:
*    pt3 - {широта, долгота} точки Q3
*/
int SphereAngular(double pt1[], double pt2[], double azi13, double azi23,
  double pt3[])
{
  double azi12, dist12, azi21, dist13;
  double cos_beta1, sin_beta1, cos_beta2, sin_beta2, cos_dist12, sin_dist12;
  SphereInverse(pt2, pt1, &azi21, &dist12);
  SphereInverse(pt1, pt2, &azi12, &dist12);
  cos_beta1 = cos(azi13 - azi12);
  sin_beta1 = sin(azi13 - azi12);
  cos_beta2 = cos(azi21 - azi23);
  sin_beta2 = sin(azi21 - azi23);
  cos_dist12 = cos(dist12);
  sin_dist12 = sin(dist12);
  if (sin_beta1 == 0. && sin_beta2 == 0.)    // Решение - любая точка
    return -1;       //  на большом круге Q1-Q2.
  else if (sin_beta1 == 0.) {
    pt3[0] = pt2[0];       // Решение - точка Q2.
    pt3[1] = pt2[1];
    return 0;
  } else if (sin_beta2 == 0.) {       // Решение - точка Q1.
    pt3[0] = pt1[0];
    pt3[1] = pt1[1];
    return 0;
  } else if (sin_beta1 * sin_beta2 < 0.) {    // Лучи Q1-Q3 и Q2-Q3 направлены
    if (fabs(sin_beta1) >= fabs(sin_beta2)) { //  в разные полусферы.
      cos_beta2 = -cos_beta2;       // Выберем ближайшее решение:
      sin_beta2 = -sin_beta2;       //  развернём луч Q2-Q3 на 180°;
    } else {       //    иначе
      cos_beta1 = -cos_beta1;       //  развернём луч Q1-Q3 на 180°.
      sin_beta1 = -sin_beta1;
    }
  }
  dist13 = atan2(fabs(sin_beta2) * sin_dist12,
cos_beta2 * fabs(sin_beta1)
+ fabs(sin_beta2) * cos_beta1 * cos_dist12);
  SphereDirect(pt1, azi13, dist13, pt3);
  return 0;
}
</syntaxhighlight>


== Ссылки ==
== Ссылки ==
* [[Вычисление расстояния и начального азимута между двумя точками на сфере]]
* [[Вычисление расстояния и начального азимута между двумя точками на сфере]]
* [http://gis-lab.info/qa/biangulation.html Нахождение точки пересечения двух линий по углам и двум известным точкам (биангуляция)]
* [http://gis-lab.info/qa/biangulation.html Нахождение точки пересечения двух линий по углам и двум известным точкам (биангуляция)]

Версия от 13:17, 12 марта 2014

Эта страница является черновиком статьи.


Линейная засечка — это нахождение положения точки по координатам двух исходных пунктов и значениям азимутов направлений с этих пунктов на определяемую точку.

Общие положения

В качестве модели Земли принимается сфера с радиусом R, равным среднему радиусу земного эллипсоида. Аналогом прямой линии на плоскости является геодезическая линия на поверхности. На сфере геодезическая линия — дуга большого круга.

Введём следующие обозначения:

  • φ — географическая широта,
  • λ — географическая долгота,
  • α — азимут дуги большого круга,
  • σ — сферическое расстояние (длина дуги большого круга, выраженная в долях радиуса шара).

Линейное расстояние по дуге большого круга s связано со сферическим расстоянием σ формулой s = R σ.

Постановка задачи

Исходные данные
координаты пунктов Q₁, Q₂ — φ₁, λ₁, φ₂, λ₂,
начальные направления с пунктов Q₁, Q₂ на точку Q₃ — α₁₃, α₂₃.
Определяемые величины
координаты точки Q₃ — φ₃, λ₃.

Алгоритм

Файл:Sph ang.png
Угловая засечка

Решение любого вида засечек сводится к нахождению полярных координат искомой точки, т.е. начального направления и расстояния на неё с одного из исходных пунктов. На конечном этапе координаты находятся из решения прямой геодезической задачи. Поскольку в угловой засечке направления α₁₃ и α₂₃ уже заданы, остаётся определить расстояние σ₁₃ или σ₂₃.

На рисунке синим цветом выделены заданные элементы сферического треугольника, красным цветом неизвестные, зелёным — вспомогательные элементы. Очевидно, в треугольнике QQQ₃ нет ни одного известного элемента. Однако из решения обратной геодезической задачи для пунктов Q₁, Q₂ могут быть получены расстояние σ₁₂, а также азимуты α₁₂ и α₂₁, после чего углы β₁ и β₂ вычисляются как разности азимутов при соответствующих пунктах. Далее из решения треугольника QQQ₃ найдём сторону σ₁₃.

Последовательность действий:

  1. решить обратную геодезическую задачу для Q₁, Q₂: по φ₁, λ₁, φ₂, λ₂ получить α₁₂, α₂₁, σ₁₂;
  2. вычислить углы β₁, β₂;
  3. в треугольнике QQQ₃ по σ₁₂, β₁, β₂ вычислить σ₁₃;
  4. решить прямую геодезическую задачу для Q₁, Q₃: по φ₁, λ₁, α₁₃, σ₁₃ вычислить φ₃, λ₃.

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

Углы β₁, β₂ и длина σ₁₃ вычисляются по формулам:

Правда, до вычисления длины σ₁₃ необходимо проанализировать полученные значения углов β₁ и β₂. Ниже в коде функции можно увидеть пример такого анализа:

  • если линии QQ₃ и QQ₃ совпадают с QQ₂, решение не определено, т.к. решением может быть любая точка геодезической линии QQ₂;
  • если одна из линий QQ₃ и QQ₃ совпадает с QQ₂, а другая нет, решением является пункт, из которого выходит другая;
  • если линии QQ₃ и QQ₃ уходят в разные полушария от QQ₂, функция находит ближайшее к QQ₂ «ложное пересечение» этих линий.

Здесь необходимо пояснить, что на сфере две несовпадающие геодезические линии всегда пересекаются в двух точках-антиподах. В традиционной постановке задачи направление на нужное пересечение задаётся явно. Если же прямое и обратное направления по условию равнозначны, возникает вопрос выбора одного из антиподов: φ₃, λ₃ или φ₃′ = −φ₃, λ₃′ = λ₃ ± 180°.

Пример программной реализации

/*
* Решение угловой засечки
*
* Аргументы исходные:
*     pt1 - {широта, долгота} пункта Q1
*     pt2 - {широта, долгота} пункта Q2
*
* Аргументы определяемые:
*     pt3 - {широта, долгота} точки Q3
*/
int SphereAngular(double pt1[], double pt2[], double azi13, double azi23,
		  double pt3[])
{
  double azi12, dist12, azi21, dist13;
  double cos_beta1, sin_beta1, cos_beta2, sin_beta2, cos_dist12, sin_dist12;

  SphereInverse(pt2, pt1, &azi21, &dist12);
  SphereInverse(pt1, pt2, &azi12, &dist12);
  cos_beta1 = cos(azi13 - azi12);
  sin_beta1 = sin(azi13 - azi12);
  cos_beta2 = cos(azi21 - azi23);
  sin_beta2 = sin(azi21 - azi23);
  cos_dist12 = cos(dist12);
  sin_dist12 = sin(dist12);

  if (sin_beta1 == 0. && sin_beta2 == 0.)     // Решение - любая точка
    return -1;				      //   на большом круге Q1-Q2.
  else if (sin_beta1 == 0.) {
    pt3[0] = pt2[0];			      // Решение - точка Q2.
    pt3[1] = pt2[1];
    return 0;
  } else if (sin_beta2 == 0.) {		      // Решение - точка Q1.
    pt3[0] = pt1[0];
    pt3[1] = pt1[1];
    return 0;
  } else if (sin_beta1 * sin_beta2 < 0.) {    // Лучи Q1-Q3 и Q2-Q3 направлены
    if (fabs(sin_beta1) >= fabs(sin_beta2)) { //   в разные полусферы.
      cos_beta2 = -cos_beta2;		      // Выберем ближайшее решение:
      sin_beta2 = -sin_beta2;		      //   развернём луч Q2-Q3 на 180°;
    } else {				      //     иначе
      cos_beta1 = -cos_beta1;		      //   развернём луч Q1-Q3 на 180°.
      sin_beta1 = -sin_beta1;
    }
  }
  dist13 = atan2(fabs(sin_beta2) * sin_dist12,
		 cos_beta2 * fabs(sin_beta1)
		 + fabs(sin_beta2) * cos_beta1 * cos_dist12);
  SphereDirect(pt1, azi13, dist13, pt3);

  return 0;
}

Ссылки