Gesucht: Umrechnung WGS84 nach Lambert (Österreich)

[lameth]

Hallo

TransDat werd ich morgen mal versuchen, wenn ich wieder in der Firma bin.

Zur Änderung der Abplattung bzw. 1/f:
Die Änderung verändert "e" erst an der 10ten Kommastelle - die Auswirkung ist somit leider in den (zumindest für mich) relevanten 3 Kommastellen von Rechts/Hochwert nicht merkbar.
Trotzdem Danke für den Hinweis

Schwieriges Thema muss ich sagen, aber irgendwie interessant

 

[Wallraff]

Hallo,

ich hab' mich tatsächlich ans Codieren gemacht. Weiß aber ehrlich gesagt nicht, was ich da mache ... diese log diese sin, cos ...
Komme inzwischen zwar im Ostwert auf 30 m hin, im Nordwert fehlen mir noch so 5000 m. Die werden auch noch irgendwoher kommen

Dass die Abplattung scheinbar nichts ausmacht, kann ich mit Vorbehalt bestätigen. Wenn Du aber eine professionelle Routine schreibst, die sich auf die Landesvermessung bezieht, würde ich die Werte des
GRS80 nehmen.
Da ist die Abplattung 1 : 298,257 2221 01

In den Zentimetern (um die es noch nicht geht) bringt das mitnehmen einiger Nachkommastellen bei Lambda0, 13. 3333333333.


Grüße Wallraff

 

[lameth]

So,

zumindest für GPS -> Lambert (neu) stimmen die Daten der Funktion mit denen von TransDat nahezu überein.
TransDat: 451.102,259 | 504.583,949
SQL: 451.102,328 | 504.584,038

Immerhin ein Fortschritt - auch wenn die TransDat Werte nicht mit Geoland.at überein stimmen ...
Und der Rest (aka Lambert(alt)) bringt leider auch keine befriedigenden Näherungen

Werd mich am Wochenende mal über die 3-Zonen-GK Berechnung schmeißen und schaun, ob's dort weniger frustrierend ist

 

[Wallraff]

Hallo,

Die üblichen Programmierfehler ... die man nicht findet.
Ichhabsgewusstichhabsgewusst

Bei mir ist n aus sin(phi) nicht gleich n aus den Logarithmen.
Dem Nordwert fehlen immer noch 5000 m. Ostwert - driftend ...

Grüße Wallraff

 

[lameth]

Bei mir ist n aus sin(phi) nicht gleich n aus den Logarithmen.

oje, das hab ich noch gar nicht geprüft

 

[JonnyMueller]

Hallo Lameth,

ich hatte einige Zeit meine e-Mails nicht gelesen, daher meine verspätete Antwort. Hier ist meine Vorgehensweise, ich weiß nicht, ob sie richtig ist. Ich schrieb ja schon, dass insbersondere meine Datum-Transformation von WGS84 ins Datum Austria verbesserungsbedürftig ist.

Ich mache die Helmert-Transformation mit folgenden 3 (nicht 7) Parametern, die ich (glaube ich) irgendwo in Wikipedia gefunden habe:
x-Verschiebung = 595
y-Verschiebung = 87
z-Verschiebung = 473
Weiterhin lege ich dem Datum Austria das Bessel-Ellipsoid zugrunde.

Die von mir verwendeten Parameter zur Lambert-Projektion habe ich in meinen vorangehenden Postings erwähnt.

Ich hoffe, das hilft Dir weiter.

Gruß,
Jonny

 

[Wallraff]

Hallo,

also bei mir war's ein Klammerfehler (bei mehr als zweien verheddere ich mich...).
Jetzt komme ich zwar aufs gleiche wie Transdat beim WGS84.
Ich war schon nahe dran, aus der Veröffentlichung von Fröhlich neu zu programmieren, da darin eine Größe klar als "isometrische Breite" zu erkennen war, die in der anderen Lösung aus Pennsylvania total vermurkst ausschaut. Passt aber mit Transdat.

Passt nicht mit geoland ...

Eine Transformation mit drei statt sieben Parametern kann nicht die Lösung sein. Aber soweit sind wir ja noch gar nicht.

Grüße Wallraff

 

[KajakFun]

Hi,
mein kleiner Vorschlag zur Fehlersuche. Postet bzw. vergleicht doch mal eure Zwischenergebnisse
z.B. WGS84 kartesisch xyz vor und nach der Helmerttransformation um zu checken ob der Rechenweg bis dorthin wenigstens ok war um etwailige Unstimmigkeiten auf die nachfolgende Lambert Projektion einzugrenzen zu können.

Ich habe mit einem BasicProgramm mal die ersten Schritte rechnen lassen
phi=47.07472 lam=12.69417 h=0 (Höhe)
WGS84
x=4245242.080
y=956252.7291
z=4647426.0087
Helmert Transformation x'y'z' im Bessel Ellipsoid
x' =4244663.122
y'=956153.559
z'=4646944.307
Grade im Bessel Ellipsoid
phi=47.075036
lam=12.694571

(Zur Lambert Projektion habe ich auch noch ein paar Berechnungsformeln gefunden...im Anhang Seite98)
http://www.icao.int/pbn/docs/eurocontrolwgsman24.pdf

Gruss KajakFun

 

[Wallraff]

Hallo,

das zieht ja Kreise

Ich meine, beim Datumsübergang sind wir ja noch nicht.
Ich hab's mir jetzt einfach gemacht und nicht nochmal die Formeln von
Fröhlich
programmiert, habe die 1.Auflage, sondern das dort gegebene Beispiel gerechnet: passt auf den Zentimeter.
Muss allerdings mit dem sin(Phi0) rechnen, sonst so 80 m Differenz, hmm....

Unsere Lambert-Koordinatenberechnung stimmt - nur nicht in Österreich :hilfe:

Gruß nach Wien...
Wallraff

 

[lameth]

Holla, da hab ich wohl einen Stein ins rollen gebracht

Also:

Bei mir ist n aus sin(phi) nicht gleich n aus den Logarithmen.

Das passt bei mir zum Glück. Mehr oder weniger.
Aus den log's: 0.7373626303
sin(phi0): 0,737277336811944
Ich führ das jetzt mal auf Kommastellen zurück - eventuell werd ich die Anzahl noch erhöhen.

Hallo Lameth,
ich hatte einige Zeit meine e-Mails nicht gelesen, daher meine verspätete Antwort.

Trotzdem Danke, dass du dich hier gemeldet hast

Hi,
mein kleiner Vorschlag zur Fehlersuche. Postet bzw. vergleicht doch mal eure Zwischenergebnisse
z.B. WGS84 kartesisch xyz vor und nach der Helmerttransformation um zu checken ob der Rechenweg bis dorthin wenigstens ok war um etwailige Unstimmigkeiten auf die nachfolgende Lambert Projektion einzugrenzen zu können.
Ich habe mit einem BasicProgramm mal die ersten Schritte rechnen lassen
phi=47.07472 lam=12.69417 h=0 (Höhe)
WGS84
x=4245242.080 4246165.3212218741
y=956252.7291 956460.6916856348
z=4647426.0087 4648443.5268492708
Helmert Transformation x'y'z' im Bessel Ellipsoid
x' =4244663.122 4245081.4942760000
y'=956153.559 956753.0399810000
z'=4646944.307 4648342.4378560000
Grade im Bessel Ellipsoid
phi=47.075036 47.080049300787728000
lam=12.694571 12.701061066082460000

Gute Idee - hier mal meine Werte (in Rot).
Wie man sieht eigentlich nur leicht unterschiedlich, was lustigerweise schon mit den kartesischen Koordinaten im Basisellipsoid beginnt
Hast du hier andere Formeln als ich?
Meine:

E² = A² - B² / A²
N = A / sqrt(1 - E² * sin²(LAT))
X = (N + H) * cos(LAT) * cos(LONG)
Y = (N + H) * cos(LAT) * sin(LONG)
Z = ((1 - E²) * N + H) * sin(LAT)
mit
A = 6378137.000
B = 6356752.314
H = 1389.592
F = 1 / 298.257223563

Die beiden Links werd ich mir mal in der Mittagspause zu Gemüte führen

 

[lameth]

Gut, ich kann mal soweit ein Update geben:

-) Die Berechnung für Lambert (neu) - sprich rein die L.-Projektion liefert bis auf einige Dezimeter das selbe wie TransDat.
Kann man soweit als erledigt betrachten und bedeutet für mich, dass der Projektionsalgo passt.
-) Für Lambert(alt) hab ich eine Korrektur in den Projektionsparametern gefunden. Ändert man PHI0 von 47.5 auf 48.0 fällt auch die Differenz von rd. 55000km an der Y-Achse weg.

Mein - derzeit ( ) - letztes Problem liegt somit in der Umformung der geografischen zu kartesischen Koordinaten und wieder zurück.
Der o.a. Punkt von KajakFun liefert mir mit meiner Umformung rd. 500m Differenz - die Projektion mit seinen transformierten Koordinaten rd. 25m.
Da mein Vorgang exakt der selbe wie bei Fröhlich/Körner ist steh ich im Moment etwas in der Luft und es wär nett, wenn du mir sagen könntest,
welche Formeln du dabei verwendet hast.
Ich denke, der Fehler liegt gleich in der ersten Umrechnung von geografisch(wgs84) -> kartesisch(wgs84) - da haben wir ~900m Diff.,
was sich dann wohl über den ganzen Rest weiterzieht.

 

[KajakFun]

Hi lameth,
ich habe mit Höhe 1389.592 die selben x,y,z Werte wie du bekommen aber nach der Helmert Transformation unterscheiden sie sich.
x=4246165.321
y=956460.691
z=4648443.526
Nach der Transformation mit den Parametern
cx=-577.326
cy=-90.129
cz=-463.919
s=-2.423E-6
rx=5.137/3600/180*pi rem in radians
ry=1.474/3600/180*pi
rz=5.297/3600/180*pi
x'=cx+(1+s)*(x-rz*y+ry*z)
y'=cy+(1+s)*(rz*x+y-rx*z)
z'=cz+(1+s)*(-ry*x+rx*y+z)
bekomme ich jedoch andere Werte

x'=4245586.36
y'=956361.52
z'=4647961.82
Für die Lambert Projektion erhalte ich Rechts=351508.004 Hoch=352974.3
Die Formel für die Proj. muss ich aber erst noch durchchecken.

 

[lameth]

Oh man, oh man, oh man ....
Ich hab vergessen die Drehwinkel der Helmertparameter in Radiant umzurechnen


Jetzt passts bei 10 Testpunkten meist auf <1m

Ich danke euch ganz herzlich

 

[KajakFun]

:silvester_mrgreen:

 

[pfeffer]

Aber ich poste mal gerne den Code der Projektionsberechnung, möglicherweise fällt ja jemandem etwas auf:

[...]
/* dependent Variables */
set @t = tan(pi()/4 - @LAT / 2) / power((1.0 - (@e * sin(@LAT))) / (1.0 + (@e * sin(@LAT))), @e/2)
set @m = cos(@LAT) / sqrt(1.0 - square(@e) * square(sin(@LAT)))
set @R = @a * @F0 * power(@t, @n)

Schon lustig: @m wird gar nicht verwendet. Das ist in der Veröffentlichung, von der Du die Formeln genommen hast, genauso... komisch.
Es mir nur aufgefallen, weil ich es grad in Cachewolf programmiere...
Hattest Du noch irgendwelche Fehler in dieser Projektionsroutine entdeckt?

Gruß,
Pfeffer.

 

[pfeffer]

o Mann, Mann, mann.
Mir scheint, die Vorzeichen der Transformationsparamter, die von der EU hier ( http://www.crs-geo.eu/nn_124226/crseu/EN/CRS__Description/crs-national__node.html?__nnn=true ) zu finden sind, sind irgendwie nicht konsequent.
Wenn ich die shift-Parameter von Gauß-Küger DHDN nehme, um nach WGS84 bzw ETRS89 zu konvertieren, passt alles. Wenn ich die von Österreich nehme, dann muss ich das Vorzeichen aller Shiftparameter umdrehen.
Sollte man meinen, die haben ne ordentliche, einheitliche Datenbank da aufgebaut...

EDIT: Die Ergebnisse von Transdat stimmen nicht mit denen hier von KajakFun angegebenen überein :-( KajakFun scheint die Skalierung anderherum gemacht zu haben als Transdat das macht. Wenn man die Koos von KajakFun 12.69417 / 47.07472 einsetzt und in Breite/Länge in Bessel Hermanskogel umrechnet kommt 12.6947940705 / 47.0752078589 raus, in Lambert(alt) 351525.077 297425.104

Gruß,
Pfeffer.

 

[pfeffer]

das ist ja alles seltsam.
Um die Lambert-Werte zu erhalten, die geoland.at bei Klick auf den x/y-Button umrechnet, muss man als zentrale Breite 48° verwenden.
In EPSG:31287 ist jedoch 47,5° angegeben und der WMS-Server von http://www.centropemap.org/ verwendet auch 47,5° zur Projektion.

Noch was: Ich habe nochmal Formeln gefunden, für (mutßamlich) alle Projektionen der Welt und zwar hin- und zurück: http://www.epsg.org/guides/docs/G7-2.pdf

EDIT: und noch eine schöne Quelle für Projektionsformeln: http://www.remotesensing.org/geotiff/proj_list/

Gruß,
Pfeffer.

 

[Wallraff]

Hallo,

angenehm, sich hier am Schreibtisch die Quellen liefern zu lassen.

Sehe, dass geoland.at seine Präsentation geändert hat. Hatte damals gerätselt, warum Lambert(neu) nicht passt. Hab im Moment nicht so das Interesse daran, Druckfehler zu finden und eigene Programmierfehler zu machen.

Grüße Wallraff

 

[KajakFun]

Ich habe momentan leider auch wenig Zeit mich wieder in die Materie einzuarbeiten.
Schau mal hier nach wegen der Referenzdaten für die Umrechnung WGS84->Austria und den Lambert Projektion Referenzmeridianen.
http://www.bev.gv.at/pls/portal/doc...ENT_ALLGEMEIN/0200_PRODUKTE/PDF/KOORD-SYS.PDF
Wenn die Helmert Transformation mit den korrekten Parametern trotzdem andere Ergebnisse nach der Transformation liefert kann evtl. ein Fehler noch an den teilweise negativen Vorzeichen in der Abbildungsmatrix liegen. Checke mal die Rechenvorschrift in der Transformationsmatrix Berechnung nach.
Unter http://vogis.cnv.at/ kann man unter Vorarlberg Atlas
für x/Y Koordinaten Umrechnung durchführen (Leiste X/Y)
Die Eingabe der WGS84 Grade12,694169 und 47,074719 liefert genau dieselben Lambert(alt) Projektions Werte 351525 und 297425.
Ich nehme dann an dass das auch so korrekt ist.

Projektion Rechtswert Hochwert
GK M28 179325,65 217783,85
GK M31 -48490,63 215274,96
Lambert 351525,04 352993,21
Lambert (alt) 351525,04 297425,04
BMN M31 401509,37 5215274,96
WGS84 12,694169 47,074719
UTM 32N 780445,68 5220092,06

Viel Glück

 

[pfeffer]

Vielen Dank für die Antworten!

Postet bzw. vergleicht doch mal eure Zwischenergebnisse (...) Ich habe mit einem BasicProgramm mal die ersten Schritte rechnen lassen
phi=47.07472 lam=12.69417 h=0 (Höhe)
WGS84
x=4245242.080
y=956252.7291
z=4647426.0087
Helmert Transformation x'y'z' im Bessel Ellipsoid
x' =4244663.122
y'=956153.559
z'=4646944.307
Grade im Bessel Ellipsoid
phi=47.075036
lam=12.694571

ich bekomme einschließlich 2 Nachkommastellen (=cm) exakt das gleiche Ergbnis wie Du.
in Lat/Lon umgerechnet stimmen die Ergebnisse bis auf alle von Dir angegeben Stellen überein.
Allerdings liefert TransDat ein abweichendes Ergebnis:
12.6947940705 und 47.0752078589 (WGS84 auf geografische Koordinaten mit "MGI (AT/CU), Hermannskogel, Bessel"
Das ist eine Abweichung von ca. 2/10.000 Grad, was etwa 20m entspricht.

Diese Abweichung findet sich dann auch in den Lambert-projizierten Koordinaten:
351525.07 297425.10 (Transdat)
351508.04 297406.16 (ich)
-----------------------
-17m -19m

meine sind also um 19m südlicher und 17m westlicher als die von Transdat.

Wenn man unsere Lon/Lat (E 12.694571358127053 N 47.07503610847066) in TransDat auf Lambert projizieren lässt, erhält man bis auf cm exakt das gleiche Ergebnis wie ich, wenn man mit 48° als Zentralem Breitengrad rechnet. In der von Dir in Deinem letzen Beitrag angegebenen Quelle steht jedoch 47,5°.

Für mich stellen sich jetzt vorallem 3 Fragen:
1. Woher kommt die Abweichug bei der Helmert-Transformation von ca. 20m?
2. wird in Österreich nun mit 47,5° oder mit 48° gerechnet? - stimmt die Spezifikation mit 47,5° oder die Implementierungen von den Landesvermessungsämtern (zb. geoland.at oder http://vogis.cnv.at/ oder TransDat)?
3. in welchen Koordinatensystem sind gedruckte Karten in Österreich?

wegen 1. habe ich mal die einzelnen Komponenten ex, ey, ez auf 0 gesetzt, um die Stärke und Richtung festzustellen:
ex = -0,0000249048787985968
mit ex = 0 ergibt sich eine Verschiebung nach
Osten um 120m und nach Norden um -26m
ey = -7,14615365955456E-06
mit ey = 0 ergibt sich eine Verschiebung nach
Osten um 7m und nach Norden um 44m
ez = -0,0000256805806883721
mit ez = 0 ergibt sich eine Verschiebung nach
Osten um -112m und nach Norden um 1m

Daraus entnehme ich, dass eine einfach Änderung von Vorzeichen vermutlich nicht die Korrektur um die nötigen ca. 20m bringen wird.
Sind meine ex, ey, ez falsch?

Vielen Dank!
Pfeffer.