Программирование перевода координат из одной системы в другую.

пространственная (прямоугольная) система координат не подразумевает какую либо проекцию, и какого либо геоида она просто не зависит от этих параметров.

 

тогда у нас просто аффиные преобразования, а получить коэффициенты аффиного преобразования - элементарная задача

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

 

Оффтоп

Вот эта как раз задача и элементарная, стала для меня не элементарной. Именно в этом и вопрос.

--- Сообщения объединены, 22 июл 2015, Оригинальное время сообщения: 22 июл 2015 ---

Я вполне точно сформулировал проблему. Тем самым определил задачу которую я решаю. Вы уточнили до абсолютного минимума мою задачу . (Найти коэффициенты аффинного преобразования). А затем по непонятным причинам переопределили задачу как поиска датума по точкам. Это сходная, но иная задача. Которую я не решаю.

Оффтоп

Понятие «Датум» используется в геодезии и картографии для наилучшей аппроксимации к геоиду в данном месте. Датум задается смещением референц-эллипсоида по осям: X, Y, Z, а также поворотом декартовой системы координат в плоскости осей на угол rX, rY, rZ. Также необходимо знать параметры референц-эллипсоида а и f, где а — размер большой полуоси, f — сжатие эллипсоида.

 

http://gis-lab.info/qa/affine-math.html

 

Ок. Сужаю проблему. Эту всю теорию я знаю. Параметры находятся по минимально необходимым 3 точкам. Но у меня массив данных. Расчленять массив на три локальных искать среднее и работать как с тремя точками это не совсем верно. (вернее не подходит для меня)

 

либо Полиномиальные преобразования
либо
1. находим коэф по первым 3м точкам
2. получаем отклонение для остальных
3. берём 3 точки с максимальным отклонением и находим коэф по ним
как то так...

 

Параметры связи двух пространственных (прямоугольных) систем координат. Причем в общем случаи. Не привязываясь к ИГД и проекциям. Найти коэффициенты аффинного преобразования. Для примера ситуация в офтопе

Оффтоп

Допустим мы тахом отсняли облако точек, но в конце работы определили, что он не был нивелирован по горизонту, далее мы отсняли такое же облако точек тахом в другом положении, теперь нам надо найти разницу между положением таха в первом и втором случаи то есть найти параметры dX,dY,dZ wX,wY,wZ (ну и масштабный коэффициент) такую ситуацию трудно представить на практике и ее практическое применение. Но она целиком объясняет поставленную мной задачу.

 

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

 

Метод и алгоритмы определения параметров преобразования между различными системами координат применительно к задачам обработки спутниковых измерений / Е.П. Алексашин, А.М. Ширенин // Геодезия и картография. - 2002. - №4. - С. 4-26.

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

 

Добрый день!
Реализовал алгоритм, приведенный в ГОСТе, для пересчета координат плоских прямоугольных из СК-42 в геодезические WGS-84.
Результат проверил программкой archaeoSYS: вычисленные широта и долгота совпадают с точностью до 9 знака после запятой. Но вот с высотой явно что-то не то. Изначально она задавалась равной нулю, после пересчета же она приняла отрицательное значение, причем довольно большое (это видно на скриншоте).
В алгоритме я ошибок не нашел, проверил несколько раз. Кто-нибудь сталкивался с подобным? Может быть где-то в самом ГОСТе ошибка?

 

Без куска кода понять в чём дело невозможно.

 

Код:

import static java.lang.Math.PI;
import static java.lang.Math.abs;
import static java.lang.Math.asin;
import static java.lang.Math.cos;
import static java.lang.Math.pow;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;
import static java.lang.Math.tan;
import static java.lang.Math.toDegrees;
import static java.lang.Math.toRadians;
 
public class TransCord3 {
public static void main(String[] args) {
// Плоские прямоугольньные координаты определяемой точки системы СК-42 в проекции Гаусса-крюгера [м]
double x = 6099177;
double y = 7309763;
 
// 1 - преобразование плоских прямоугольных координат системы СК-42 в геодезические
// (в проекции Гаусса-Крюгера на эллипсоиде Красовского)
 
// 1я вспомогательная величина
double beta = x/6367558.4968;
// геодезическая широта точки, абсцисса которой равна абсциссе заданной точки, а ордината равна нулю
double b0 = beta + sin(2*beta)*(0.00252588685 - 0.00001491860*pow(sin(beta),2) + 0.00000011904*pow(sin(beta),4));
// номер 6-градусной зоны
int zoneNumber = (int)(y*pow(10,-6));
// 2я вспомогательная величина
double z0 = (y - (10*zoneNumber+5)*pow(10,5))/(6378245*cos(b0));
// поправка к геодезической широте
double deltaB = -pow(z0,2)*sin(2*b0)*(0.251684631 - 0.003369263*pow(sin(b0),2) + 0.000011276*pow(sin(b0),4)
-pow(z0,2)*(0.10500614 - 0.04559916*pow(sin(b0),2) + 0.00228901*pow(sin(b0),4) - 0.00002987*pow(sin(b0),6)
-pow(z0,2)*(0.042858 - 0.025318*pow(sin(b0),2) + 0.014346*pow(sin(b0),4) - 0.001264*pow(sin(b0),6)
-pow(z0,2)*(0.01672 - 0.00630*pow(sin(b0),2) + 0.01188*pow(sin(b0),4) - 0.00328*pow(sin(b0),6)))));
// расстояние от заданной точки до осевого меридиана зоны
double l = z0*(1 - 0.0033467108*pow(sin(b0),2) - 0.0000056002*pow(sin(b0),4) - 0.0000000187*pow(sin(b0),6)
-pow(z0,2)*(0.16778975 + 0.16273586*pow(sin(b0),2) - 0.00052490*pow(sin(b0),4) - 0.00000846*pow(sin(b0),6)
-pow(z0,2)*(0.0420025 + 0.1487407*pow(sin(b0),2) + 0.0059420*pow(sin(b0),4) - 0.0000150*pow(sin(b0),6)
-pow(z0,2)*(0.01225 + 0.09477*pow(sin(b0),2) + 0.03282*pow(sin(b0),4) - 0.00034*pow(sin(b0),6)
-pow(z0,2)*(0.0038 + 0.0524*pow(sin(b0),2) + 0.0482*pow(sin(b0),4) + 0.0032*pow(sin(b0),6))))));
// вычисление геодезической широты
double lat = b0 + deltaB;
// вычисление геодезической долготы
double lon = 6*(zoneNumber - 0.5)/57.29577951 + l;
// высота
double height = 0;
System.out.printf("Геодезические координаты в системе СК-42:%n%.9f%n%.9f%n%.9f%n%n",
toDegrees(lat), toDegrees(lon), height);
 
// 2 - преобразование геодезических координат системы СК-42 в пространственные прямоугольные координаты
 
// параметры эллипсоида Красовского
int a = 6378245; // большая полуось
double alfa = 1/298.3; // сжатие
// квадрат эксцентриситета эллипсоида
double e2 = 2*alfa - pow(alfa, 2);
// радиус кривизны первого вертикала
double n = a/sqrt(1 - e2*pow(sin(lat), 2));
// пространственные прямоугольные координаты системы СК-42
double x1 = (n+height)*cos(lat)*cos(lon);
double y1 = (n+height)*cos(lat)*sin(lon);
double z1 = (n*(1-e2)+height)*sin(lat);
System.out.printf("Пространственные прямоугольные координаты в системе СК-42:%n%f%n%f%n%f%n%n", x1, y1, z1);
 
// 3 - преобразование пространственных прямоугольных коодинат из системы СК-42 в систему ПЗ-90.02
 
// элементы трансформирования систем координат при переходе между системами СК-42 и ПЗ-90.02
// линейные 
double deltaX = 23.93;
double deltaY = -141.03;
double deltaZ = -79.98;
// угловые
double omegaX = 0;
double omegaY = toRadians(-0.35/3600);
double omegaZ = toRadians(-0.79/3600);
// масштабный элемент трансформирования 
double m = -0.22*pow(10, -6);
 
// матрица угловых элементов трансформирования систем координат
double[] mx1 = {	  1,  omegaZ, -omegaY,
-omegaZ,	   1,  omegaX,
 omegaY, -omegaX,	   1};
for(int i = 0; i < mx1.length; i++)
mx1 = mx1*(1+m);
// пространственные прямоугольные координаты системы ПЗ-90.02
double x2 = mx1[0]*x1 + mx1[1]*y1 + mx1[2]*z1 + deltaX;
double y2 = mx1[3]*x1 + mx1[4]*y1 + mx1[5]*z1 + deltaY;
double z2 = mx1[6]*x1 + mx1[7]*y1 + mx1[8]*z1 + deltaZ;
 
// 4 - преобразование пространственных прямоугольных коодинат из системы ПЗ-90.02 в систему WGS-84
 
// линейные элементы трансформирования систем координат при переходе между системами ПЗ-90.02 и WGS-84
deltaX = -0.36;
deltaY = 0.08;
deltaZ = 0.18;
// пространственные прямоугольные координаты системы WGS-84
double x3 = x2 + deltaX;
double y3 = y2 + deltaY;
double z3 = z2 + deltaZ;
System.out.printf("Пространственные прямоугольные координаты в системе WGS-84:%n%f%n%f%n%f%n%n", x3, y3, z3);
 
// 5 - преобразование пространственных прямоугольных координат системы WGS-84 в геодезические координаты
 
// параметры общеземного эллипсоида в системе WGS-84
a = 6378137; // большая полуось
alfa = 1/298.257223563; // сжатие
// квадрат эксцентриситета эллипсоида
e2 = 2*alfa - pow(alfa, 2);
// вспомогательная величина
double d = sqrt(pow(x3, 2) + pow(y3, 2));
// анализ значения d
if(d == 0) {
lat = PI*z3/(2*abs(z3));
lon = 0;
height = z3*sin(lat) - a*sqrt(1 - e2*pow(sin(lat), 2));
}
else {
double lonA = abs(asin(y3/d));
if(y3 < 0 & x3 > 0) lon = 2*PI - lonA;
if(y3 < 0 & x3 < 0) lon = PI + lonA;
if(y3 > 0 & x3 < 0) lon = PI - lonA;
if(y3 > 0 & x3 > 0) lon = lonA;
if(y3 == 0 & x3 > 0) lon = 0;
if(y3 == 0 & x3 < 0) lon = PI;
// аналих хначения z3
if(z3 == 0) {
lat = 0;
height = d - a;
}
else {
// вычисление вспомогательных величин
double r = sqrt(pow(x3, 2) + pow(y3, 2) + pow(z3, 2));
double c = asin(z3/r);
double p = e2*a/(2*r);
double s1 = 0;
double b = c + s1;
double s2 = asin(p*sin(2*b)/sqrt(1 - e2*pow(sin(b), 2)));
d = abs(s2 - s1);
// допуск прекращения итеративного процесса
double precision = 0.000000001;
// итеративный процесс для вычисления геодезической широты и высоты
while(d >= precision) {
s1 = s2;
b = c + s1;
s2 = asin((p*sin(2*b))/sqrt(1 - e2*pow(sin(b), 2)));
d = abs(s2 - s1);
}
lat = b;
height = d*cos(lat) + z3*sin(lat) - a*sqrt(1 - e2*pow(sin(lat), 2));
}
}
System.out.printf("Геодезические координаты в системе WGS-84:%n%.9f%n%.9f%n%.9f%n%n",
toDegrees(lat), toDegrees(lon), height);
} 
}

 

Ошибка в операторе

нужное значение
double d = sqrt(pow(x3, 2) + pow(y3, 2));
чуть выше вы затираете в цикле в операторе
d = abs(s2 - s1);

Так, к слову
вместо double m = -0.22*pow(10, -6);
проще double m = -0.22E-6;
Вместо pow(x3, 2) и т.п.
лучше x3*x3

Вместо e2 = 2*alfa - pow(alfa, 2);
лучше e2 = alfa*(2 - alfa);
Куча if – это тяжкое наследие FORTRAN II (середина 60-х прошлого века)
когда не было встроенной функции atan2

Я бы рекомендовал алгоритм Морозова из методички Валеры и Тамары Луповок
Паскаль процедура

Код:

procedure XYZ2BLH(X, Y, Z: extended; SemiMajorAxis, Flattening: extended;
  var B, L, H: extended);
var
  D, a, f, ee, c, N: extended;
  p, p1: extended;
  k, k1: extended;
  t, t0, t1: extended;
const
  eps: extended=1E-13;
begin
  a := SemiMajorAxis;
  f := Flattening;
  ee := f*(2-f);
  c := a/(1-f);
  D := sqrt(sqr(X)+sqr(Y));
  if D>Abs(Z) then
  begin
	t0 := Z/D;
	p := c*ee/D;
	k := 1/(1-ee);
	t := t0;
	repeat
	  t1 := t;
	  t := t0+p*t1/sqrt(k+sqr(t1));
	until (Abs(t1-t)<eps);
	N := c*sqrt((1+sqr(t))/(k+sqr(t)));
	H := D*sqrt(1+sqr(t))-N;
	B := arctan(t);
  end
  else
  begin
	t0 := D/Z;
	p1 := -a*ee/Z;
	k1 := 1-ee;
	t := t0;
	repeat
	  t1 := t;
	  t := t0+p1*t1/sqrt(k1+sqr(t1));
	until (Abs(t1-t)<eps);
	N := a*sqrt((1+sqr(t))/(k1+sqr(t)));
	H := Abs(Z)*sqrt(1+sqr(t))-k1*N;
	B := Pi/2-arctan(t);
  end;
  L := ArcTan2(Y, X);
 
end; // procedure XYZ2BLH

 

stout, еще раз спасибо )
С высотой разобрался. Реализовал оба алгоритма в итоге, расхождения там совсем небольшие получаются.
Если из гостовского алгоритма эту кучу if выкинуть, он даже проще морозовского будет. Но это уже вопрос к тем, кто ГОСТ составлял.

 

Вообще-то алгоритмов перехода от XYZ к BLH дюжины две, есть среди них и аналитически строгие, без итераций. Основная идея последних – сведение нелинейного уравнения относительно широты к алгебраическому уравнению четвёртой степени. Вот только особенность современных массовых процессоров такова, что итерационные алгоритмы оказываются быстрее. (Кроме того, операция деления всё ещё долгая. На большинстве современных процессоров умножение чисел double занимает один такт, как и сложение. А вот деление – 19 тактов. Ровно столько, сколько извлечение квадратного корня)
Забавно также и то, что алгоритм Морозова за бугром не известен. В алгоритме Морозова особое внимание уделено точности вычислений. То есть, когда широта меньше 45° (не строгое значение) то находится тангенс, когда больше, то котангенс.
Заслуженной популярностью за рубежом пользуется итерационный алгоритм Боуринга. Это такая хитровывернутая запись алгоритма Ньютона для решения нелинейного уравнения. (см. статью Fukushima во вложении)

http://www.mygeodesy.id.au/documents/Transforming Cartesian Coordinates.pdf
https://geodesy.curtin.edu.au/local/docs/1Featherstone295.pdf

 

С учетом того, что преобразовывать координаты нужно будет в реальном времени, это то что нужно.

 

Уважаемый stout, а Вы не в курсе, как они это делают? Дело в том, что я попробовал реализовать программно такое же поведение: полные матрицы разворотов вокруг каждой из осей при 7 параметрическом преобразовании, честные синусы и косинусы. Но увы, получил облом. В отладке вижу, что то, что не упростил я, за меня округлил процессор. Такое впечатление, что мощности встроенного типа для хранения чисел с плавающей точкой (Delphi Extended 10 байт, 20 значащих разрядов мантиссы) не хватает, и процессор округляет результаты до того же угла, и той же единицы. Может у Вас есть какая либо информация на эту тему, способная натолкнуть на мысль, как это можно победить?
Спасибо.

 

80 бит (long double, мантисса 52 бит) достаточно для самых точных геодезических расчетов в пределах планеты Земля.
Для sin() или tg() малых углов вместо значения угла (в радианах, естественно) вы можете взять разложение тригонометрических функций в ряд и учесть второй член разложения.

(added)
По моей оценке, чтобы при использовании упрощенного синуса вместо точного вылезла ошибка в 0.000001м , разворот осей должен быть > 20", но т.к. это единичная операция, а в матрицах таких операций много, возможно для суперточных расчетов второй член ряда Тейлора имеет смысл.

 

Хватает, хватает. Возьмите параметры, где углы разворота не доли секунды, а секунды, тогда увидите разницу.

Чуток поменьше. 19.2
В отличии от double (64 бита, стандарт IEE-754), где в нормализованном числе первая единица всегда подразумевается, поэтому точность 53 бита мантиссы, а не 52, в long double (или паскалевский extended, extended precision в терминах IEE 754) мантисса 64 бита.
Но работать с такими числами в операционных системах с вытесняющей многозадачностью (а это и Windows и Linux) надо очень осторожно. Попробуйте периодически проверять точность с помощью дельфийской GetPrecisionMode – будете сильно удивлены. Особенно после работы различных кодеров-декодеров видио и аудио.
Не зря MS для своих трансляторов устанавливает, что точность long double≡double.
То, что в ляйковском софте вычисления идут по точным формулам, заметил лет 20 назад мой друг Володя Третьяков. Тогда написали вместо преобразования Бурша-Вольфа строгие формулы, и результаты стали полностью совпадать с ляйковскими.
З.Ы. log(2⁶⁴)=19.26