На главную... Продукты | Технологии | Классификаторы | Проекты | Скачать | Цены| Форум | Статьи | Обучение | Контакты

Функции mapDistance, mapRealDistance, mapInversePositionComputation, mapConventionalDistance

Поиск  Пользователи  Правила  Войти
Форум » Linux » Средства разработки ГИС-приложений для Linux
Страницы: Пред. 1 2
RSS
Функции mapDistance, mapRealDistance, mapInversePositionComputation, mapConventionalDistance
 
Еще раз спасибо!
Всё заработало, но - работает не на всем пространстве Земли почему-то. Не могли бы Вы посмотреть, что даёт Ваша функция для следующих данных:
Точка 1: широта 57.873058, долгота 27.918778
Точка 2: широта 53,677778, долгота 75,706944

Моя трансляция описанных выше по теме функций здесь дает какое-то ошибочное расстояние в 5206 километра. Референсное расстояние по сфере 2963 километра.

Так же прикладываю табличку, наглядно демонстрирующую, как с определенного расстояния моя трансляция Вашей функции начинает давать что-то непонятное. В этой табличке даны: Широта и Долгота точки2 (точка 1 не меняется), референсное расстояние по сфере, рассчитанное функцией ReverseGeoTask расстояние, разница с референсным, рассчитанный азимут. Я просто последовательно увеличиваю долготу на 2 градуса - т.е. по идее азимут должен плавно уменьшаться. А расстояние поначалу плавно уменьшаться, а потом плавно увеличиваться. Видно, что это работает до определенного момента

Изменено: Иван - 27.06.2017 13:02:48
 
Здравствуйте !

Согласен, в коде была ошибка которую я  исправил очень давно но тут на форуме не обновлял

В сообщении выше я внёс правку. Заклбчаетс яона в следующем
строчку в функции ReversGeoTask
Код
aLen:=D.X/(P0.b*cos(dA12));
нужно исправить вот так  
Код
aLen:=Abs(D.X/(P0.b*cos(dA12))); // <<< + Abs
Не тот глуп кто не знает, а тот, кто не знает где искать.
 
Да, так стало лучше. Но где-то всё-равно чего-то не хватает. Возможно, у меня - поэтому и прошу проверить, какие данные дает Ваша реализация для "ошибочных" точек. К примеру, беру такие точки:
Точка 1: широта 57.873058, долгота 27.918778
Точка 2: широта 28,918922, долгота -13,783472

Расчет дает 4064 км и азимут 142 на вторую точку. А референс (да и карты) дают 4536 км и 243 градуса. Или вот иллюстрация - у обеих точек ставлю широту 10 градусов и долготу 0 градусов. Начинаю двигать вторую точку по долготе, получается что при нахождении обеих точек на широте 10, с разницей по долготе 180 градусов - имеем расстояние 19735 км и азимут 90, чего по идее не может быть - есть ведь более короткий путь (азимут).



P.S. Извините, что беспокою - делаю программу для своего хобби, радиоприема.
Изменено: Иван - 27.06.2017 19:36:40
 
Цитата
P.S. Извините, что беспокою - делаю программу для своего хобби, радиоприема.
Чесно говоря работу функции на расстояниях больших чем 1000 км я не проверял, поэтому могут быть какие то неточности с дальностью

Хотя я могу пересмотреть работу алгоритма для разных полушарий и больших расстояний и выложить сюда новый код..
Наверное на выходных так и сделаю
Не тот глуп кто не знает, а тот, кто не знает где искать.
 
Иван, предлагаю второй вариант решения обратной геодезической задачи.
Код полностью переписан
Код
function ReversGeoTask(const P1,P2 : T3dPoint): T3dPoint;
var L,lam, dlam : double;
    uSq, A, B, C, sig, dsig : double;
    sinsig, cossig, sinSqsig : double;
    cos2sigM  : double;
    iterationLimit : integer;
    cosU1, cosU2 : double;
    sinU1, sinU2 : double;
begin
  L :=  P2.Y - P1.Y;
  A := (1-eWGS84.rf) * tan(P1.B);
  B := (1-eWGS84.rf) * tan(P2.B);
  cosU1 := 1 / sqrt((1 + A*A));
  cosU2 := 1 / sqrt((1 + B*B));
  sinU1 := A * cosU1;
  sinU2 := B * cosU2;
  lam := L;  dlam := 0;
  iterationLimit  := 250;
  repeat
    sinSqsig := (cosU2*sin(lam)) * (cosU2*sin(lam)) + (cosU1*sinU2-sinU1*cosU2*cos(lam)) * (cosU1*sinU2-sinU1*cosU2*cos(lam));
    if (sinSqsig = 0) then exit; 
    sinsig   := sqrt(sinSqsig);
    cossig   := sinU1*sinU2 + cosU1*cosU2*cos(lam);
    sig      := arctan2(sinsig, cossig);
    B     := cosU1 * cosU2 * sin(lam) / sinsig;
    A        := 1 - B*B;
    if (A <> 0) then
      cos2sigM := cossig - 2*sinU1*sinU2/A
    else
      cos2sigM := 0;
    C    := eWGS84.rf/16*A*(4+eWGS84.rf*(4-3*A));
    dlam := lam;
    lam  := L + (1-C) * eWGS84.rf * B * (sig + C*sinsig*(cos2sigM+C*cossig*(-1+2*cos2sigM*cos2sigM)));
    dec(iterationLimit);
  until (abs(lam-dlam)<=1E-12) or (iterationLimit<1);

  if (iterationLimit=0) then  exit;

  uSq  := (1 - B*B) * (eWGS84.a*eWGS84.a - eWGS84.b * eWGS84.b) / (eWGS84.b*eWGS84.b);
  A    := 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  B    := uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
  dsig := B*sinsig*(cos2sigM+B/4*(cossig*(-1+2*cos2sigM*cos2sigM)- B/6*cos2sigM*(-3+4*sinsig*sinsig)*(-3+4*cos2sigM*cos2sigM)));

  result.X := eWGS84.b*A*(sig-dsig);
  result.Y := arctan2(cosU2*sin(lam),  cosU1*sinU2-sinU1*cosU2*cos(lam));
  result.Z := arctan2(cosU1*sin(lam), -sinU1*cosU2+cosU1*sinU2*cos(lam));
  if result.Y<0 then result.Y := result.Y+2*pi;
  if result.Z<0 then result.Z := result.Z+2*pi;
end;

Изменено: KFF - 28.06.2017 10:28:31
Не тот глуп кто не знает, а тот, кто не знает где искать.
 
прямая задача
Код
function TrueGeoTask(Bo, Lo, Azim, Dist: double): T3dPoint;
var sinAzim, cosAzim  : double;
    tanU1,cosU1,sinU1 : double;
    sinA,cosSqA,A,B    : double;
    sg1,sig,dsg,psig,uSq  : double;
    cos2sM ,sinSg, cosSg: double;

begin
  sinAzim := sin(Azim); 
  cosAzim := cos(Azim);
  tanU1   := (1-eWGS84.rf) * tan(Bo);
  cosU1   := 1 / sqrt(1 + tanU1*tanU1);
  sinU1   := tanU1 * cosU1;
  sg1     := arctan2(tanU1, cosAzim);
  sinA    := cosU1 * sinAzim;
  cosSqA  := 1 - sinA*sinA;
  uSq     := cosSqA*(eWGS84.a*eWGS84.a-eWGS84.b*eWGS84.b)/(eWGS84.b*eWGS84.b);
  A       := 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  B       := uSq/1024*(256+uSq*(-128+uSq*(74-47*uSq)));
  sig     := Dist/(eWGS84.b * A); psig := 0;
  repeat
    cos2sM := cos(2*sg1 + sig);
    sinSg  := sin(sig);
    cosSg  := cos(sig);
    dsg    := B*sinSg*(cos2sM+0.25*B*(cosSg*(-1+2*cos2sM*cos2sM)-
              B/6*cos2sM*(-3+4*sinSg*sinSg)*(-3+4*cos2sM*cos2sM)));
    psig   := sig;
    sig    := Dist / (eWGS84.b*A) + dsg;
  until (abs(sig-psig) <= 1E-12);
  A   := eWGS84.rf/16*cosSqA*(4+eWGS84.rf*(4-3*cosSqA));
  B   := sinU1*sinSg - cosU1*cosSg*cosAzim;
  result.B := arctan2(sinU1*cosSg + cosU1*sinSg*cosAzim, (1-eWGS84.rf)*sqrt(sinA*sinA + B*B));
  result.L  := Lo+ arctan2(sinSg*sinAzim, cosU1*cosSg - sinU1*sinSg*cosAzim) -
               (1-A) * eWGS84.rf * sinA * (sig + A*sinSg*(cos2sM+A*cosSg*(-1+2*cos2sM*cos2sM)));
  result.Z := arctan2(sinA, -B);
end;
Не тот глуп кто не знает, а тот, кто не знает где искать.
 
Приветствую!

Тоже возникла потребность в расчетах без открытого HMAP.
Может есть смысл в АПИ добавить функции расчетов расстояний и прочего для HUSER, например. Велосипедостроительством заниматься очень не хочется..

С уважением,
Матвеев П.В.
 
Расчет расстояния - это сложная функция. Она строит ортодромию для больших расстояний с учетом эллипсоида карты,
затем каждый отрезок трансформирует в свою проекцию и считает длины.
HMAP можно получить для временной карты и служебного классификатора (service.rsc).
 
Приветствую!

Спасибо за ответ
Но ведь все необходимые параметры для расчета вроде в HUSER есть.. Задача посчитать расстояние не для какой-то произвольной системы, а для общепринятых. Проще говоря, по epsg коду получил параметры - рассчитал расстояние. Ведь преобразование координат для huser реализовано. Понятно, что можно создать временную карты по параметрам полученным по epsg коду, но не хотелось лишних действий выполнять.

С уважением,
Матвеев П.В.
Страницы: Пред. 1 2
Читают тему (гостей: 1)



© КБ Панорама, 1991-2024

Регистрируясь или авторизуясь на форуме, Вы соглашаетесь с Политикой конфиденциальности