Швейцарские координаты проекции на эллипсоидальные координаты - PullRequest
Новая
       43

Швейцарские координаты проекции на эллипсоидальные координаты

0 голосов
/ 13 сентября 2018

Я пытаюсь преобразовать координаты CH1903 + в координаты WGS84 на основе следующего примера, предоставленного швейцарским "Bundesamt für Statistik" с использованием C #.

Пример:

Example

Насколько я могу рассчитать все значения, как показано в примере.Но в конце, когда я пытаюсь вычислить переменную «S» на основе значений в примере, я получаю неправильные результаты. 

double E = 2679520.05; 
double alpha = 1.00072913843038; 
double phi = 0.820535226; 
double b = 0.820535226;
double K = 0.0030667323772751

Я попробовал обе реализации: 

double S = Math.Log(Math.Tan(Math.PI / 4.0 + phi / 2.0)); --> result: 0.931969600979248

или 

double S = 1/alpha * (Math.Log(Math.PI/4 + b/2) - K) + E * Math.Log(Math.Tan((Math.PI/4) + (Math.Asin(E * Math.Sin(phi))/2.0))); --> result: NaN

Может кто-нибудь сказать, что не так в моей реализации, что я получаю неправильные результаты?Если я правильно понимаю пример, оба вычисления должны вернуть 0.933114264192610 как результат для заданных значений.

1 Ответ

0 голосов
/ 11 декабря 2018

Это моя реализация для преобразования швейцарских координат CH1903 + -> десятичных координат (WGS84), которая до сих пор хорошо работает

Пример CH1903 + Координаты N = 1249251.544, E = 2645320.487; 

// swiss CH1903+ coordinates
double N = 1249251.544,
       E = 2645320.487;

const double K = 0.0030667323772751; // Konstante der Breitenformel

const double BESSEL_GH = 6377397.155; // grosse Halbachse des Bessel-Ellipsoids 
const double BESSEL_SH = 6356078.963; // kleine Halbachse des Bessel-Ellipsoids 
double GRS_GH = 6.378137 * Math.Pow(10, 6); // grosse Halbachse des GRS84-Ellipsoids 
double GRS_SH = 6.356752314 * Math.Pow(10, 6); ; // kleine Halbachse des GRS84-Ellipsoids 

double E2_BESSEL = (Math.Pow(BESSEL_GH, 2) - Math.Pow(BESSEL_SH, 2)) / Math.Pow(BESSEL_GH, 2), // 1.numerische Exzentrizität (im Quadrat) des Bessel-Ellipsoids 
       E_BESSEL = Math.Sqrt(E2_BESSEL), // 1.numerische Exzentrizität des Bessel-Ellipsoids
       E2_GRS = (Math.Pow(GRS_GH, 2) - Math.Pow(GRS_SH, 2)) / Math.Pow(GRS_GH, 2); // 1.numerische Exzentrizität (im Quadrat) des GRS84-Ellipsoids 

const double TOLERANCE = 0.0000000000000001;

///// swiss coordinates -> ellipsoid coordinates /////
double
    Y = E - 2600000.0, // swiss CH1903 coordinates
    X = N - 1200000.0, // swiss CH1903 coordinates

    PHI_NULL = (46 + 57 / 60.0 + 8.66 / 3600.0) * Math.PI / 180, // geogr. Breite des Nullpunkts in Bern
    LAMBDA_NULL = (7 + 26 / 60.0 + 22.5 / 3600.0) * Math.PI / 180, // geogr. Länge des Nullpunkts in Bern

    R = (BESSEL_GH * Math.Sqrt(1 - E2_BESSEL)) / (1 - E2_BESSEL * Math.Pow(Math.Sin(PHI_NULL), 2)), // Radius der Projektionskugel
    ALPHA = Math.Sqrt(1 + E2_BESSEL / (1 - E2_BESSEL) * Math.Pow(Math.Cos(PHI_NULL), 4)), // Verhältnis Kugellänge zu Ellipsoidlänge

    B_GLOBAL_NULL = Math.Asin(Math.Sin(PHI_NULL) / ALPHA), // Breite des Nullpunkts auf der Kugel

    L_PSEUDO = Y / R, // Kugelkoordinaten bezüglich Pseudoäquatorsystem in Bern
    B_PSEUDO = 2 * (Math.Atan(Math.Pow(Math.E, (X / R))) - (Math.PI / 4)), // Kugelkoordinaten bezüglich Pseudoäquatorsystem in Bern

    L = Math.Atan(Math.Sin(L_PSEUDO) / (Math.Cos(B_GLOBAL_NULL) * Math.Cos(L_PSEUDO) - Math.Sin(B_GLOBAL_NULL) * Math.Tan(B_PSEUDO))), // Kugelkoordinaten bezüglich Nullpunkt Bern 
    B = Math.Asin(Math.Cos(B_GLOBAL_NULL) * Math.Sin(B_PSEUDO) + Math.Sin(B_GLOBAL_NULL) * Math.Cos(B_PSEUDO) * Math.Cos(L_PSEUDO)), // Kugelkoordinaten bezüglich Nullpunkt Bern 

    LAMBDA = LAMBDA_NULL + L / ALPHA,
    PHI = B,
    S = (Math.Log(Math.Tan(Math.PI / 4.0 + B / 2.0)) - K) / ALPHA + E_BESSEL * Math.Log(Math.Tan(Math.PI / 4.0 + Math.Asin(E_BESSEL * Math.Sin(PHI)) / 2.0));

bool cont = true;

// iterate to tolerance
while (cont)
{
    double PHI_TMP = 2 * Math.Atan(Math.Pow(Math.E, S)) - Math.PI / 2;
    double S_TMP = (Math.Log(Math.Tan(Math.PI / 4 + B / 2)) - K) / ALPHA + E_BESSEL * Math.Log(Math.Tan(Math.PI / 4 + Math.Asin(E_BESSEL * Math.Sin(PHI_TMP)) / 2));

    if (Math.Abs(PHI - PHI_TMP) < TOLERANCE)
    {
        cont = false;
    }
    else
    {
        PHI = PHI_TMP;
        S = S_TMP;
    }
}

///// ellipsoid coordinates (CH1903) -> karth. coordinates (CH1903) /////
double RN_CH1903 = BESSEL_GH / (Math.Sqrt(1 - E2_BESSEL * Math.Pow(Math.Sin(PHI), 2))), // Normalkrümmungsradius
       X_KARTH = RN_CH1903 * Math.Cos(PHI) * Math.Cos(LAMBDA),
       Y_KARTH = RN_CH1903 * Math.Cos(PHI) * Math.Sin(LAMBDA),
       Z_KARTH = (RN_CH1903 * (1 - E2_BESSEL)) * Math.Sin(PHI);

///// karth. coordinates (CH1903) -> karth. coordinates (ETRS89/CHTRS95) /////
double X_CH1930P = X_KARTH + 674.374,
       Y_CH1930P = Y_KARTH + 15.056,
       Z_CH1930P = Z_KARTH + 405.346;

///// karth. coordinates (ETRS89/CHTRS95) -> ellipsoid coordinates (ETRS89/CHTRS95) /////
double PHI_ETRS = Math.Atan(Z_CH1930P / Math.Sqrt(Math.Pow(X_CH1930P, 2) + Math.Pow(Y_CH1930P, 2))),
       RN_ETRS = BESSEL_GH / (Math.Sqrt(1 - E2_GRS * Math.Pow(Math.Sin(PHI_ETRS), 2)));

cont = true;

// iterate to tolerance
while (cont)
{
    double PHI_ETRS_TMP = Math.Atan((Z_CH1930P / (Math.Sqrt(Math.Pow(X_CH1930P, 2) + Math.Pow(Y_CH1930P, 2)))) / (1 - (RN_ETRS * E2_GRS) / (RN_ETRS + 0))),
           RN_ETRS_TMP = BESSEL_GH / (Math.Sqrt(1 - E2_GRS * Math.Pow(Math.Sin(PHI_ETRS_TMP), 2)));

    if (Math.Abs(PHI_ETRS - PHI_ETRS_TMP) < TOLERANCE)
    {
        cont = false;
    }
    else
    {
        PHI_ETRS = PHI_ETRS_TMP;
        RN_ETRS = RN_ETRS_TMP;
    }
}

double LAMBDA_ETRS = Math.Atan(Y_CH1930P / X_CH1930P);

///// ellipsoid coordinates (ETRS89/CHTRS95) -> decimal coordinates (WGS84) /////
double lat = PHI_ETRS * 180 / Math.PI;
double lon = LAMBDA_ETRS * 180 / Math.PI;
...