Реализации алгоритмов/Перевод географических координат в прямоугольные в прямоугольные координаты проекции Гаусса-Крюгера

Здесь приведена реализация перевода географических координат в прямоугольные координаты проекции Гаусса-Крюгера.

w:Python

# Перевод географических координат (широты и долготы) точки в прямоугольные
# координаты проекции Гаусса-Крюгера (на примере координат Москвы).

# Географические координаты точки (в градусах)
dLon = 37.618  # Долгота (положительная для восточного полушария)
dLat = 55.752  # Широта (положительная для северного полушария)

# Номер зоны Гаусса-Крюгера (если точка рассматривается в системе
# координат соседней зоны, то номер зоны следует присвоить вручную)
zone = int(dLon/6.0+1)

# Импорт математических функций
from math import sin, cos, tan, pi

# Параметры эллипсоида Красовского
a = 6378245.0          # Большая (экваториальная) полуось
b = 6356863.019        # Малая (полярная) полуось
e2 = (a**2-b**2)/a**2  # Эксцентриситет
n = (a-b)/(a+b)        # Приплюснутость

# Параметры зоны Гаусса-Крюгера
F = 1.0                   # Масштабный коэффициент
Lat0 = 0.0                # Начальная параллель (в радианах)
Lon0 = (zone*6-3)*pi/180  # Центральный меридиан (в радианах)
N0 = 0.0                  # Условное северное смещение для начальной параллели
E0 = zone*1e6+500000.0    # Условное восточное смещение для центрального меридиана

# Перевод широты и долготы в радианы
Lat = dLat*pi/180.0
Lon = dLon*pi/180.0

# Вычисление переменных для преобразования
v = a*F*(1-e2*(sin(Lat)**2))**-0.5
p = a*F*(1-e2)*(1-e2*(sin(Lat)**2))**-1.5
n2 = v/p-1
M1 = (1+n+5.0/4.0*n**2+5.0/4.0*n**3)*(Lat-Lat0)
M2 = (3*n+3*n**2+21.0/8.0*n**3)*sin(Lat-Lat0)*cos(Lat+Lat0)
M3 = (15.0/8.0*n**2+15.0/8.0*n**3)*sin(2*(Lat-Lat0))*cos(2*(Lat+Lat0))
M4 = 35.0/24.0*n**3*sin(3*(Lat-Lat0))*cos(3*(Lat+Lat0))
M = b*F*(M1-M2+M3-M4)
I = M+N0
II = v/2*sin(Lat)*cos(Lat)
III = v/24*sin(Lat)*(cos(Lat))**3*(5-(tan(Lat)**2)+9*n2)
IIIA = v/720*sin(Lat)*(cos(Lat)**5)*(61-58*(tan(Lat)**2)+(tan(Lat)**4))
IV = v*cos(Lat)
V = v/6*(cos(Lat)**3)*(v/p-(tan(Lat)**2))
VI = v/120*(cos(Lat)**5)*(5-18*(tan(Lat)**2)+(tan(Lat)**4)+14*n2-58*(tan(Lat)**2)*n2)

# Вычисление северного и восточного смещения (в метрах)
N = I+II*(Lon-Lon0)**2+III*(Lon-Lon0)**4+IIIA*(Lon-Lon0)**6
E = E0+IV*(Lon-Lon0)+V*(Lon-Lon0)**3+VI*(Lon-Lon0)**5

print ('Широта:            ', dLat)
print ('Долгота:           ', dLon)
print ('Северное смещение: ', N)
print ('Восточное смещение:', E)

w:Java

// Перевод географических координат (широты и долготы) точки в прямоугольные
// координаты проекции Гаусса-Крюгера (на примере координат Москвы).

// Географические координаты точки (в градусах)
double dLon = 37.618;
double dLat = 55.752;

// Номер зоны Гаусса-Крюгера (если точка рассматривается в системе
// координат соседней зоны, то номер зоны следует присвоить вручную)
int zone =  (int) (dLon/6.0+1);

// Параметры эллипсоида Красовского
double a = 6378245.0;          // Большая (экваториальная) полуось
double b = 6356863.019;        // Малая (полярная) полуось
double e2 = (Math.pow(a,2)-Math.pow(b,2))/Math.pow(a,2);  // Эксцентриситет
double n = (a-b)/(a+b);        // Приплюснутость


// Параметры зоны Гаусса-Крюгера
double F = 1.0;                   // Масштабный коэффициент
double Lat0 = 0.0;                // Начальная параллель (в радианах)
double Lon0 = (zone*6-3)* Math.PI/180;  // Центральный меридиан (в радианах)
double N0 = 0.0;                  // Условное северное смещение для начальной параллели
double E0 = zone*1e6+500000.0;    // Условное восточное смещение для центрального меридиана

// Перевод широты и долготы в радианы
double Lat = dLat*Math.PI/180.0;
double Lon = dLon*Math.PI/180.0;

// Вычисление переменных для преобразования
double sinLat = Math.sin(Lat);
double cosLat = Math.cos(Lat);
double tanLat = Math.tan(Lat);

double v = a * F * Math.pow(1-e2* Math.pow(sinLat,2),-0.5);
double p = a*F*(1-e2) * Math.pow(1-e2*Math.pow(sinLat,2),-1.5);
double n2 = v/p-1;
double M1 = (1+n+5.0/4.0* Math.pow(n,2) +5.0/4.0* Math.pow(n,3)) * (Lat-Lat0);
double M2 = (3*n+3* Math.pow(n,2) +21.0/8.0* Math.pow(n,3)) * Math.sin(Lat - Lat0) * Math.cos(Lat + Lat0);
double M3 = (15.0/8.0* Math.pow(n,2) +15.0/8.0* Math.pow(n,3))*Math.sin(2 * (Lat - Lat0))*Math.cos(2 * (Lat + Lat0));
double M4 = 35.0/24.0* Math.pow(n,3) *Math.sin(3 * (Lat - Lat0)) * Math.cos(3 * (Lat + Lat0));
double M = b*F*(M1-M2+M3-M4);
double I = M+N0;
double II = v/2 * sinLat * cosLat;
double III = v/24 * sinLat * Math.pow(cosLat,3) * (5-Math.pow(tanLat,2)+9*n2);
double IIIA = v/720 * sinLat * Math.pow(cosLat,5) * (61-58*Math.pow(tanLat,2)+Math.pow(tanLat,4));
double IV = v * cosLat;
double V = v/6 * Math.pow(cosLat,3) * (v/p-Math.pow(tanLat,2));
double VI = v/120 * Math.pow(cosLat,5) * (5-18*Math.pow(tanLat,2)+Math.pow(tanLat,4)+14*n2-58*Math.pow(tanLat,2)*n2);

// Вычисление северного и восточного смещения (в метрах)
double N = I+II* Math.pow(Lon-Lon0,2)+III* Math.pow(Lon-Lon0,4)+IIIA* Math.pow(Lon-Lon0,6);
double E = E0+IV*(Lon-Lon0)+V* Math.pow(Lon-Lon0,3)+VI* Math.pow(Lon-Lon0,5);

System.out.println("Широта:            "+ dLat);   //Широта:            55.752
System.out.println("Долгота:           "+ dLon);   //Долгота:           37.618
System.out.println("Северное смещение: "+ N);      //Северное смещение: 6181924.245933299
System.out.println("Восточное смещение:"+ E);      //Восточное смещение:7413223.481449484