brauche hilfe, bei GPS>GK umrechnung

[spunky]

ich hab mich jetzt mal hier angemeldet, da ich woanders nicht viel infos finden konnte, die mir wirklich weiterhelfen...

ich brauche für berufliche und private zwecke eine software (pocketpc) zur flächenberechnung. es gibt billige (bis 100€) software die nur ungenaue ergebnisse liefern, weil sie nur eine messung pro meßpunkt vornehmen. dann gibt es teure software (ab 500€) die sehr gute ergebnisse liefern, weil sie einen durchschnitt von mehreren messungen pro meßpunkt liefern. da mir das aber viel zu teuer ist (da ich es nur gelegentlich nutze aber dennoch eine hohe genauigkeit brauche), wollte ich das selber schreiben.

im prinzip funktioniert schon alles, wenn ich die GPS-koordinaten nach der exel-tabelle von moenk (http://www.geoclub.de/viewtopic.php?t=6233) in GK umrechnen lasse und die fläche dann mit der gaußschen formel berechne.

benutze ich aber folgende einfache formel (http://www.dsdt.info/tipps/?id=627), erhalte ich total falsche GK-koordinaten.

    Dim pi, rho As Single
    pi = 4 * Atn(1)
    rho = 180 / pi
    Dim brDezimal, laDezimal, rm, e2, c, bf, g, co, g2, g1, t, dl, fa, grad, min, sek, sy
    Dim rw, hw
    sy = 3
    e2 = 0.0067192188
    c = 6398786.849

    brDezimal = Text1 ' 51,123456
    laDezimal = Text2 ' 12,123456

    bf = brDezimal / rho
    g = 111120.61962 * brDezimal - 15988.63853 * sIn(2 * bf) + 16.72995 * sIn(4 * bf) - 0.02178 * sIn(6 * bf) + 0.00003 * sIn(8 * bf)
    co = Cos(bf)
    g2 = e2 * (co * co)
    g1 = c / Sqr(1 + g2)
    t = sIn(bf) / Cos(bf) '{=tan(t)}
    dl = laDezimal - sy * 3
    fa = co * dl / rho
    hw = g + fa * fa * t * g1 / 2 + fa * fa * fa * fa * t * g1 * (5 - t * t + 9 * g2) / 24
    rm = fa * g1 + fa * fa * fa * g1 * (1 - t * t + g2) / 6 + fa * fa * fa * fa * fa * g1 * (5 - 18 * t * t * t * t * t * t) / 120
    rw = rm + sy * 1000000 + 500000

    Text3 = Int(rw)
    Text4 = Int(hw)

ich habe schon versucht die formeln aus der exel-tabelle aufzulösen und in den code einzubauen, aber die formel ist einfach zu komplex, ich schaff das leider nicht.

hat vielleicht einer schon die tabelle in code umgewandelt bekommen und kann mir weiterhelfen? bzw. warum funktioniert der code oben nicht, kann man den vielleicht abändern so das die richtigen GKs bei rauskommen?

achso, das ganze wird in VB bzw. eVB geschrieben...

mfg.

 

[rms62]

Hi,

hier meine Umrechnungen Lat/Lon in GK und GK in Lat/Lon

Function gk(lat As Double, lon As Double, h As Single) As String
Dim pi, a, b, e, n, AL, be, ga, de, ep, e2, x, y, z, x1, y1, z1, s, t, L, B1, h1, L1, L0, B3, Bogenlänge
Dim i, BL1, BL2, BL3, BL4, Hochwert, RW1, RW2, RW3, RW4, Rechtswert
'Bessel-Ellipsoid

pi = 3.14159265358979
a = 6377397.155
b = 6356078.962
h = 4.21
e = (a ^ 2 - b ^ 2) / (a ^ 2)
n = (a - b) / (a + b)
AL = (a + b) / 2 * (1 + 1 / 4 * n ^ 2 + 1 / 64 * n ^ 4) 'Alpha
be = -3 / 2 * n + 9 / 16 * n ^ 3 - 3 / 32 * n ^ 5 'Beta
ga = 15 / 16 * n ^ 2 - 15 / 32 * n ^ 4 'Gamma
de = -35 / 48 * n ^ 3 + 105 / 256 * n ^ 5 'Delta
ep = 315 / 512 * n ^ 4 'Epsilon
a = 6378137
b = 6356752.314
e2 = (a ^ 2 - b ^ 2) / a ^ 2
n = a / Sqr(1 - e2 * Sin(lat / 180 * pi) ^ 2)
x = (n + h) * Cos(lat / 180 * pi) * Cos(lon / 180 * pi)
y = (n + h) * Cos(lat / 180 * pi) * Sin(lon / 180 * pi)
z = (n * b ^ 2 / a ^ 2 + h) * Sin(lat / 180 * pi)
x1 = x * 1 + y * -0.0000119021759 + z * -0.000000218166156
y1 = x * 0.0000119021759 + y * 1 + z * 0.000000979323636
z1 = x * 0.000000218166156 + y * -0.0000009793236 + z * 1

x = x1 * 0.9999933 + (-598.095)
y = y1 * 0.9999933 + (-73.707)
z = z1 * 0.9999933 + (-418.197)

a = 6377397.155
b = 6356078.962
e2 = (a ^ 2 - b ^ 2) / a ^ 2
s = Sqr(x ^ 2 + y ^ 2)
t = Atn(z * a / (s * b))
b = Atn((z + e2 * a ^ 2 / b * Sin(t) ^ 3) / (s - e2 * a * Cos(t) ^ 3))
L = Atn(y / x)
n = a / Sqr(1 - e2 * Sin(b) ^ 2)
h = s / Cos(b) - n
B1 = b * 180 / pi
L1 = L * 180 / pi
h1 = h
If Abs(L1 - 6) < 1.5 Then
    L0 = 6
ElseIf Abs(L1 - 9) < 1.5 Then
    L0 = 9
ElseIf Abs(L1 - 12) < 1.5 Then
    L0 = 12
Else
    L0 = 15
End If
b = 6356078.962
i = (L1 - L0) * pi / 180
B3 = B1 / 180 * pi
n = a / Sqr(1 - e2 * Sin(B3) ^ 2)
h1 = Sqr(a ^ 2 / b ^ 2 * e2 * Cos(B3) ^ 2)
t = Tan(B3)
Bogenlänge = AL * (B3 + be * Sin(2 * B3) + ga * Sin(4 * B3) + de * Sin(6 * B3) + ep * Sin(8 * B3))
BL1 = t / 2 * n * Cos(B3) ^ 2 * i ^ 2
BL2 = t / 24 * n * Cos(B3) ^ 4 * (5 - t ^ 2 + 9 * h1 ^ 2 + 4 * h1 ^ 4) * i ^ 4
BL3 = t / 720 * n * Cos(B3) ^ 6 * (61 - 58 * t ^ 2 - 330 * t * h1 ^ 2) * i ^ 6
BL4 = t / 40320 * n * Cos(B3) ^ 8 * (1385 - 3111 * t ^ 2 + 543 * t ^ 4 - t ^ 6) * i ^ 8
Hochwert = Bogenlänge + BL1 + BL2
RW1 = n * Cos(B3) * i
RW2 = n / 6 * Cos(B3) ^ 3 * (1 - t ^ 2 + h1 ^ 2) * i ^ 3
RW3 = n / 120 * Cos(B3) ^ 5 * (5 - 18 * t ^ 2 + t ^ 4 + 14 * h1 ^ 2 - 58 * t ^ 2 * h1 ^ 2) * i ^ 5
RW4 = n / 5040 * Cos(B3) ^ 7 * (61 - 479 * t ^ 2 + 179 * t ^ 4) * i ^ 7
Rechtswert = RW1 + RW2 + 500000 + L0 / 3 * 1000000
gk = Hochwert & "RW" & Rechtswert

End Function

Function LatLon(HW As Double, RW As Double, h As Single) As String
Dim pi, a, b, e, n, AL, be, ga, de, ep, e2, x, y, z, x1, y1, z1, s, t, L, B1, h1, L1, L0, B3, Bogenlänge
Dim i, BL1, BL2, BL3, BL4, Hochwert, RW1, RW2, RW3, RW4, Rechtswert, y0, b0, nf, tf, bf, pif
Dim tf1, tf2, tf3, tf4, bw, l2, l3, l4, lw, lat As Double, lon As Double
'Bessel-Ellipsoid
pi = 3.14159265358979
a = 6377397.155
b = 6356078.962
h = 4.21
e = (a ^ 2 - b ^ 2) / (a * a)
n = (a - b) / (a + b)  'nicht a*b???
AL = (a + b) / 2 * (1 + 1 / 4 * n ^ 2 + 1 / 64 * n ^ 4)
be = 3 / 2 * n - 27 / 32 * n ^ 3 + 269 / 512 * n ^ 5
ga = 21 / 16 * n ^ 2 - 55 / 32 * n ^ 4
de = 151 / 96 * n ^ 3 - 417 / 128 * n ^ 5
ep = 1097 / 512 * n ^ 4
y0 = Int(RW / 1000000)
L0 = y0 * 3
y = RW - y0 * 1000000 - 500000
b0 = HW / AL
bf = (b0 + be * Sin(2 * b0) + ga * Sin(4 * b0) + de * Sin(6 * b0) + ep * Sin(8 * b0))
nf = a / Sqr(1 - e * Sin(bf) ^ 2)
pif = Sqr(a ^ 2 / b ^ 2 * e * Cos(bf) ^ 2)
tf = Tan(bf)
tf1 = tf / 2 / nf ^ 2 * (-1 - pif ^ 2) * y ^ 2
tf2 = tf / 24 / nf ^ 4 * (5 + 3 * tf ^ 2 + 6 * pif ^ 2 - 6 * tf ^ 2 * pif ^ 2 - 4 * pif ^ 4 - 9 * tf ^ 2 * pif ^ 4) * y ^ 4
tf3 = tf / 720 / nf ^ 6 * (-61 - 90 * tf ^ 2 - 45 * tf ^ 4 - 107 * pif ^ 2 + 162 * tf ^ 2 * pif ^ 2 + 45 * tf ^ 4 * pif ^ 2) * y ^ 6
tf4 = tf / 40320 / nf ^ 8 * (1385 + 3663 * tf ^ 2 + 4095 * tf ^ 4 + 1575 * tf ^ 6) * y ^ 8
bw = (bf + tf1 + tf2) * 180 / pi
L1 = 1 / nf / Cos(bf) * y
l2 = 1 / 6 / nf ^ 3 / Cos(bf) * (-1 - 2 * tf ^ 2 - pif ^ 2) * y ^ 3
l3 = 1 / 120 / nf ^ 5 / Cos(bf) * (5 + 28 * tf ^ 2 + 24 * tf ^ 4 + 6 * pif ^ 2 + 8 * tf ^ 2 * pif ^ 2) * y ^ 5
l4 = 1 / 15040 / nf ^ 7 / Cos(bf) * (-61 - 622 * tf ^ 2 - 1320 * tf ^ 4 - 720 * tf ^ 6) * y ^ 7
lw = L0 + (L1 + l2) * 180 / pi
a = 6377397.155
b = 6356078.962
e2 = (a ^ 2 - b ^ 2) / a ^ 2
n = a / Sqr(1 - e2 * Sin(bw / 180 * pi) ^ 2)
x = (n + h) * Cos(bw / 180 * pi) * Cos(lw / 180 * pi)
y = (n + h) * Cos(bw / 180 * pi) * Sin(lw / 180 * pi)
z = (n * b ^ 2 / a ^ 2 + h) * Sin(bw / 180 * pi)
x1 = x * 1 + y * 0.0000119021759 + z * 0.000000218166156
y1 = x * -0.0000119021759 + y * 1 + z * -0.000000979323636
z1 = x * -0.000000218166156 + y * 0.0000009793236 + z * 1

x = x1 * 0.9999933 + (598.095)
y = y1 * 0.9999933 + (73.707)
z = z1 * 0.9999933 + (418.197)
a = 6378137
b = 6356752.314
e2 = (a ^ 2 - b ^ 2) / a ^ 2
s = Sqr(x ^ 2 + y ^ 2)
t = Atn(z * a / (s * b))
b = Atn((z + e2 * a ^ 2 / b * Sin(t) ^ 3) / (s - e2 * a * Cos(t) ^ 3))
L = Atn(y / x)
n = a / Sqr(1 - e2 * Sin(b) ^ 2)
h = s / Cos(b) - n
lat = b * 180 / pi
lon = L * 180 / pi
LatLon = lat & "RW" & lon

End Function

Bis dann

Ralle

 

[Leonderprofi]

Hallo Spunky,

ich nutze aus das Programm. Bei mir kommen die Hochwerte (Latitude) ganz gut hin, doch bei den Rechtswerten (Longitude) haben ich Abweichungen von bis zu drei Kilometern. Hast du eine Lösung für das Problem gefunden. Wäre nett, wenn du dich sonst mal melden könntest!

Schönen Gruß,

Leonderprofi

 

[rms62]

Oben steht doch die Lösung mit den richtigen Ergebnissen.

Bis dann

Ralle

 

[Leonderprofi]

Hallo rms62,

Was meinst du in deinem Prog mit der Zeite y0 = Int(RW/100000). Ich muss das ganzen in Java schreiben.

Gruß,

Malte

 

[rms62]

Hallo,

das ist der Ganzzahlanteil des übergebenen Rechtswertes (RW) geteilt durch 100000.
Aus der Frage entnehme ich, dass du gegebenen Rechts- und Hochwert in Lat/Lon umwandeln willst. Es sind ja zwei Funktionen dargestellt.
1. Lat/Lon in GK
2. GK in Lat/Lon

Bis dann

Ralle

 

[Leonderprofi]

Das gibt es aber ein Problem. In der Zeile:

n = a / Sqr(1 - e2 * Sin(bw / 180 * pi) ^ 2)

wird die Wurzel negativ. Woran könnte das liegen?

 

[rms62]

Da muss vorher ein Fehler vorliegen. Ich habe selber auch lange gebraucht, bis ich die Exceltabellen so umgesetzt hatte, dass alles fehlerfrei lief. Immer wieder hatte ich an der falschen Stelle das Hochzeichen zum quadrieren falsch gesetzt usw.
Wenn du den VB-code 1:1 in Java umsetzt, müsste es eigentlich gehen.

Probier noch mal oder besser noch poste den Java-Code, dann kann ich auch mal gucken.

Bis dann

Ralle

 

[moenk]

Wenn das jemand grad in Java macht: Für die Rotation und Translation der 3D-Kartesischen Koordinaten macht es Sinn wenn die Parameter frei wählbar sind. Ich fiel grad beim durchgucken gesehen, dass keine Matrizenmultiplikation gemacht wird, da werden irgendwie Konstanten eingerechnet, oder hab ich da was übersehen?

 

[spunky]

@rms62: danke dir für deine hilfe, war sicher ein riesen aufwand die tabelle umzuschreiben, hab es ja selbst veruscht, aber dann den überblick verloren, danke nochmal.

vor allen auch an moenk ein danke.

wie ist das nun mit der matrizenmultiplikation, kann man das auch in vb machen? rms62 rechnet ja mit dem satz 2 (D2001) und soweit stimmen auch die ergebnisse überein.