Konvertierung von *.pth Dateien

[routeconverter]

Die Umrechnung erledige ich derzeit wie in der angehängten Excel-Datei beschrieben. In Java sieht das aus:

Code:

private static final double aBessel = 6377397.155;
    private static final double bBessel = 6356078.962;
    private static final double e2Bessel = (Math.pow(aBessel, 2) - Math.pow(bBessel, 2)) / Math.pow(aBessel, 2);
    private static final double nBessel = (aBessel - bBessel) / (aBessel + bBessel);

    private static final double aWgs84 = 6378137;
    private static final double bWgs84 = 6356752.314;
    private static final double e2Wgs84 = (Math.pow(aWgs84, 2) - Math.pow(bWgs84, 2)) / Math.pow(aWgs84, 2);

    private static final double alphaGk2Wgs84 = (aBessel + bBessel) / 2 * (1 + Math.pow(nBessel, 2) / 4 + Math.pow(nBessel, 4) / 64);
    private static final double betaGk2Wgs84 = nBessel * 3 / 2 - Math.pow(nBessel, 3) * 27 / 32 + Math.pow(nBessel, 5) * 269 / 512;
    private static final double gammaGk2Wgs84 = Math.pow(nBessel, 2) * 21 / 16 - Math.pow(nBessel, 4) * 55 / 32;
    private static final double deltaGk2Wgs84 = Math.pow(nBessel, 3) * 151 / 96 - Math.pow(nBessel, 5) * 417 / 128;
    private static final double epsilonGk2Wgs84 = Math.pow(nBessel, 4) * 1097 / 512;

    public static double[] gaussKruegerRightHeightToWgs84LongitudeLatitude(double right, double height) {
        if (!(right > 1000000))
            throw new IllegalArgumentException("Invalid Gauss-Krueger right value given: " + right);
        if (!(height > 1000000))
            throw new IllegalArgumentException("Invalid Gauss-Krueger height value given: " + height);

        double h = 4.21;

        // Umrechnung GK nach B, L
        int y0 = new Double(right / 1000000).intValue();
        double L0 = y0 * 3;
        int yInt = new Double(right - y0 * 1000000 - 500000).intValue();
        double B0 = height / alphaGk2Wgs84;
        double Bf = (B0 + betaGk2Wgs84 * Math.sin(2 * B0) + gammaGk2Wgs84 * Math.sin(4 * B0) + deltaGk2Wgs84 * Math.sin(6 * B0) + epsilonGk2Wgs84 * Math.sin(8 * B0));
        double Nf = aBessel / Math.sqrt(1 - e2Bessel * Math.pow(Math.sin(Bf), 2));
        double pif = Math.sqrt(Math.pow(aBessel, 2) / Math.pow(bBessel, 2) * e2Bessel * Math.pow(Math.cos(Bf), 2));
        double tf = Math.tan(Bf);
        double tf1 = tf / 2 / Math.pow(Nf, 2) * (-1 - Math.pow(pif, 2)) * Math.pow(yInt, 2);
        double tf2 = tf / 24 / Math.pow(Nf, 4) * (5 + 3 * Math.pow(tf, 2) + 6 * Math.pow(pif, 2) - 6 * Math.pow(tf, 2) * Math.pow(pif, 2) - 4 * Math.pow(pif, 4) - 9 * Math.pow(tf, 2) * Math.pow(pif, 4)) * Math.pow(yInt, 4);
        // double tf3 = tf / 720 / Math.pow(Nf, 6) * (-61 - 90 * Math.pow(tf, 2) - 45 * Math.pow(tf, 4) - 107 * Math.pow(pif, 2) + 162 * Math.pow(tf, 2) * Math.pow(pif, 2) + 45 * Math.pow(tf, 4) * Math.pow(pif, 2)) * Math.pow(yInt, 6);
        // double tf4 = tf / 40320 / Math.pow(Nf, 8) * (1385 + 3663 * Math.pow(tf, 2) + 4095 * Math.pow(tf, 4) + 1575 * Math.pow(tf, 6)) * Math.pow(yInt, 8);
        double B = (Bf + tf1 + tf2) * 180 / Math.PI;
        double l1 = 1 / Nf / Math.cos(Bf) * yInt;
        double l2 = (1 / Math.pow(Nf, 3) / 6 / Math.cos(Bf)) * (-1 - 2 * Math.pow(tf, 2) - Math.pow(pif, 2)) * Math.pow(yInt, 3);
        // double l3 = 1 / Math.pow(Nf, 5) / 120 / Math.cos(Bf) * (5 + 28 * Math.pow(tf, 2) + 24 * Math.pow(tf, 4) + 6 * Math.pow(pif, 2) + 8 * Math.pow(tf, 2) * Math.pow(pif, 2)) * Math.pow(yInt, 5);
        // double l4 = 1 / Math.pow(Nf, 7) / 15040 / Math.cos(Bf) * (-61 - 622 * Math.pow(tf, 2) - 1320 * Math.pow(tf, 4) - 720 * Math.pow(tf, 6)) * Math.pow(yInt, 7);
        double L = L0 + (l1 + l2) * 180 / Math.PI;

        // Ell. Koordinaten auf dem Bessel-Ellipsoid
        double N = aBessel / Math.sqrt(1 - e2Bessel * Math.pow(Math.sin(B / 180 * Math.PI), 2));
        double x1 = (N + h) * Math.cos(B / 180 * Math.PI) * Math.cos(L / 180 * Math.PI);
        double y1 = (N + h) * Math.cos(B / 180 * Math.PI) * Math.sin(L / 180 * Math.PI);
        double z1 = (N * Math.pow(bBessel, 2) / Math.pow(aBessel, 2) + h) * Math.sin(B / 180 * Math.PI);

        // Rotierte Vektoren
        double x2 = x1 * 1 + y1 * 0.0000119021759 + z1 * 0.000000218166156;
        double y2 = x1 * -0.0000119021759 + y1 * 1 + z1 * -0.000000979323636;
        double z2 = x1 * -0.000000218166156 + y1 * 0.0000009793236 + z1 * 1;

        // Translationen anbringen
        double x = x2 * 0.9999933 + (598.095);
        double y = y2 * 0.9999933 + (73.707);
        double z = z2 * 0.9999933 + (418.197);

        // Vektoren (in ETRF89)
        double s = Math.sqrt(Math.pow(x, 2) + Math.pow(y, 2));
        double T = Math.atan(z * aWgs84 / (s * bWgs84));
        double B2 = Math.atan((z + e2Wgs84 * Math.pow(aWgs84, 2) / bWgs84 * Math.pow(Math.sin(T), 3)) / (s - e2Wgs84 * aWgs84 * Math.pow(Math.cos(T), 3)));
        double L2 = Math.atan(y / x);
        // double N2 = aWgs84 / Math.sqrt(1 - e2Wgs84 * Math.pow(Math.sin(B2), 2));
        // h = s / Math.cos(B2) - N2;
        double latitude = B2 * 180 / Math.PI;
        double longitude = L2 * 180 / Math.PI;
        return new double[]{longitude, latitude};
    }

Fehler möchte ich da nicht ausschließen. Getestet habe ich das überwiegend mit einem Magic Maps-Nutzer aus Berlin.