• Willkommen im Geoclub - dem größten deutschsprachigen Geocaching-Forum. Registriere dich kostenlos, um alle Inhalte zu sehen und neue Beiträge zu erstellen.

Gesucht: Umrechnung WGS84 nach Lambert (Österreich)

JonnyMueller

Geonewbie
Liebe Geoclub-Freunde,

in einer Datenbank habe ich die WGS84-Koordinaten (Brteitengrad, Längengrad) mehrerer Orte in Österreich gespeichert. Aus diesen WGS84-Koordinaten möchte ich einen Link auf "geoland.at" erzeugen, mit welchem ich den jeweiligen Ort in einer Online-Karte anzeigen kann. Dazu benötige ich die österreichischen Lambert-Koordinaten dieses Ortes.

Beispiel für den Großglockner:
(1) WGS84-Koordinaten: Längengrad: 12,69417°, Breitengrad: 47,07472°
(2) österreichische Lambert-Koordinaten: x=351525, y=352993
(3) Link zum Anzeigen: http://www.geoland.at/geoland2/init.aspx?karte=geo&hotspot=351525|352993||25000

Nun habe ich keine Ahnung, wie ich von (1) nach (2) komme.

Ich habe versucht, mich schlau zu machen, konnte aber außer folgenden Angaben über die österreichischen Lambert-Koordinaten nichts in Erfahrung bringen:
- Bessel-Ellipsoid, a = 6377397,155 m, b = 6356078,963 m
- Schnittkegelprojektion
- 1. paralleler Breitengrad 46°
- 2. paralleler Breitengrad 49°
- Breitengrad des Projektionsursprungs: 48°
- Zentraler Längengrad: 13,333333°
- False Easting: 400000 m
- False Northing: 400000 m

Leider finde ich nirgendwo eine konkrete Umrechnungsformel von WGS84-Koordinaten nach österreichische Lambert-Koordinaten - also Schritt (1) nach (2) wie oben.

Hat jemand von Euch so eine Umrechnung schon einmal programmiert oder hat jemand den Algorithmus als mathematische Formel? Ich habe keine Ausbildung in Geographie oder Vermessungstechnik, daher benötige eine konkrete mathematische Umrechnungsformel oder einen Programmcode (oder diesbezügliche Internet-Links).

Über hilfreiche Antworten würde ich mich freuen.

Viele Grüße,
Jonny
 

Wallraff

Geocacher
Hallo Jonny,

noch keine Antwort ?
Ich hatte mich voriges Jahr mal dafür interessiert und ein paar versteckte Hinweise auf der Homepage des
BEV
gefunden ... die Seite ist aber nicht die alte sondern wurde neu strukturiert. :(

Hilft wenigstens die
Liste
Demnach z.B. besser 47°30' statt 48°.

Noch besser wären ein paar Hinweise, die ich aber auch nur bruchstückhaft geben kann.
Es gibt zwei Arten Lambert-Koordinaten:
die "alten" mit dem Bessel-Ellipsoid
die "neuen" mit dem WGS84-Ellipsoid. Übersichtskarte 1:500.000.

Die Formeln, wenn's die richtigen sind (ich hüte mich vor Garantien). Ich hab sie heute -endlich- gefunden:
http://surveying.wb.psu.edu/psu-surv/Projects/PASingleZone.pdf
Ziert sich etwas beim Laden.
Steht zwar Pennsylvania drauf, ist in Anhang I aber vermtl. Lambert drin.
Fürs Ellipsoid nicht für die Kugel, wobei Du halt die Bessel-Größen einsetzen musst.
Wie gesagt ohne Gewähr.

Grüße Wallraff
 
OP
J

JonnyMueller

Geonewbie
Hallo Wallraff,

vielen Dank für Deine Antwort. Ja, der Breitengrad-Ursprung liegt wohl (wie von Dir angemerkt) bei 47,5° - und nicht (wie oben von mir geschrieben) bei 48°. Danke für Deine Klarstellung.

Ich werde das Dokument mit den Formeln mal durchstudieren und mich wieder melden, was ich damit anfangen konnte.

Gruß,
Jonny
 

jmsanta

Geoguru
hier gibt es eine Umrechnungsmöglöichkeit
http://www.xs4all.nl/~estevenh/1/
aber auch dein GPS-Empfänger sollte das auf die Reihe bekommen - bei den engl. Garmin-Empfängern Kartendatum/Positionsformat:"Austrian Grid"
 

KajakFun

Geowizard
Ich hatte auch mal das Problem Vorarlberger GK Koordinaten auf WGS84 umzurechnen.
Die Garmin GPS Empfänger Routine war übrigens nicht genau genug.
Ich hatte das damals so gelöst.
1.) Umrechnung in das Swiss Grid CH1903
Die Umrechnungsparameter findet man hier (Helmerttransformation)
http://www.gis.univie.ac.at/karto/lehre/hgex/retro/hgex98/infos/bregenz.htm
2.) Umrechnung Swiss Grid nach WGS84
Die Infos zu den Berechnungen dazu findet man auf dieser Seite
http://www.swisstopo.admin.ch/internet/swisstopo/de/home/apps/calc.html
http://www.swisstopo.admin.ch/internet/swisstopo/de/home/topics/survey/sys/refsys/switzerland.parsysrelated1.24280.downloadList.65964.DownloadFile.tmp/swissprojectionde.pdf
http://www.swisstopo.admin.ch/internet/swisstopo/de/home/topics/survey/sys/refsys/switzerland.parsysrelated1.24280.downloadList.7077.DownloadFile.tmp/ch1903wgs84de.pdf
Vermutlich geht es auch direkt und einfacher...ich bin da ganz pragmatisch...Hauptsache es klappt.
Man kann dort auch Koords Listen Online umrechnen lassen.
Vielleicht hilft das etwas weiter.

Hätte ich fast vergessen..
Es gibt einen Österreichischen Luftbildserver der in der Menuleiste mit dem Button X/Y auch Koordinatenumrechnungen durchführt.
Nich schlecht zur Kontrolle der eigenen Rechenergebnisse
http://vogis.cnv.at/dva04/init.aspx
 
OP
J

JonnyMueller

Geonewbie
Lieber Wallraff, jmsanta und KajakFun,

gleich mal vielen Dank für Euere Hilfsbereitschaft. Mein Problem ist nun weitegehend gelöst. Wallraffs Tipp mit dem Pennsylvania-Dokument brachte mich in entscheidender Weise weiter. Im genannten Dokument fand ich nämlich die Umrechnung von geographischen Koordinaten auf eine konforme Kegelprojektion (d.h. sogenannte Lambert-Projektion). Ich habe den beschriebenen Algorithmus unverändert übernommen und setze lediglich die oben beschriebenen Werte ein. Die Resultate sind annähernd (meistens plus/minus 5 Meter) diejenigen, die ich in "Geoland.at" finde.

Die Umrechnung von WGS84-Koordinaten (Längengrad, Breitengrad) in österreichische Lambert-Koordinaten (Rechtswert, Hochwert) vollziehe ich in zwei Schritten:
1. "Datum Shift", d.h. Umrechnung der Koordinaten von WGS84 ins "Datum Austria"
2. Konforme Kegelprojektion, d.h. Umrechnung der "Datum Austria"-Koordinaten in Lambert-Koordinaten

Beim "Datum Shift" habe ich eine bereits vorliegende Funktion, mit der ich WGS84 ins "Potsdam-Datum" umrechne, entsprechend angepasst. Hierbei - so vermute ich - habe ich noch leichten Optimierungsbedarf. Ich werde dies aber in einer neuen Diskussion klären, falls ich diesbezüglich weitermachen möchte.

Das Feature, welches nun auf meiner "Geofinder.ch"-Website verfügbat ist, möchte ich anhand dreier Beispiele zeigen:

Wie Ihr seht, liegen die von "Geofinder.ch" (= von mir entwickelte Website) und die von "Geoland.at" (= Referenz) berechneten Lambert-Koordinaten relativ nahe zusammen.

Die jeweils entsprechende "Geoland.at"-Landkarte erhaltet Ihr, wenn Ihr auf den genannten "Geofinder.ch"-Seiten (dort auf der rechten Seite über der "Google Map") auf den Link "in Geoland.at öffnen" klickt. Probiert's einfach aus...

Wenn Ihr Verbesserungsvorschläge habet, freue ich mich auf weitere Antworten.

Viele Grüße,
Jonny
 

Wallraff

Geocacher
Hallo Jonny,

Du hast das echt programmiert ?
Ach Du SCH :zensur:
Keine Druckfehler gefunden ?

Und @ jmsanta:
Der Link ist gut, aber wohl nur auf Positionen zwischen 3 und 7,45° ö.L. ausgelegt.
Habe aber noch nicht umfassend getestet.

Grüße Wallraff
 

lameth

Geocacher
Hallo,

ich hab ein ähnliches Problem, und häng mich hier am besten gleich dran:

Als Erweiterung in unserer SQL Datenbank möchte ich eine Funktion schreiben, die GPS wgs84 Koordinaten in Lambert umrechnet.
Einerseits zum Import ins GIS und für Abstandsberechnungen, andererseits weil auch alle anderen Koordinaten in der DB in Lambert Werten vorliegen.

Folgende Schritte werden in der Funktion durchgeführt:

1.) Umrechnung der WGS84 Koordinaten auf Kartesische Koordinaten
2.) Helmert Transformation (7 Parameter) der WGS84 XYZ in Kartesische Koordinaten des Bessel Ellipsoids (sprich 'Datum Austria')
3.) Rückrechnen der Bessel XYZ in geografische Koordinaten im Bessel Ellipsoid
4.) Lambert Kegelprojektion der neuen Koordinaten

An sich funktionierts - leider sind die Ergebnisse etwas mangelhaft.
Und seltsamerweise sind die Ergebnisse ohne die Transformation ins Datum Austria - sprich rein die Lambert Projektion der WGS84 Koordinaten -
deutlich besser als mit Transformation.

Beispiel:
Das Koordinatenpaar 14.0241486° / 48.4388267°
soll lt. geoland.at 451162 / 504640 als Äquivalent haben.
Ergebnis mit Helmert: 451970.922 / 505005.306 (mehrere 100m Differenz)
Ergebnis ohne Helmert: 451113.639 / 504608.006 (<50m Differenz)

Meine Frage jetzt - kann sich das jemand erklären bzw. mir helfen?
Ich kann den SQL Code gerne online stellen, falls ihn sich jemand ansehen will.

Ich bin in dem Thema an sich ein Laie - kann aber folgendes zu den Punkten sagen:
1.) Die Formeln hab ich aus einem Mathebuch und im www überprüft - ist an sich ein simpler Schritt.
2.) ist eine einfache Matrizenauflösung. Die Parameter und die Matrix hab ich vom BEV bzw. wikipedia (übereinstimmend)
3.) ist im Prinzip ähnlich wie 1.)
zu 4.) hab ich 2 Wege programmiert, einmal einen Rechenweg aus einer Formelsammlung einer dt. Uni (stimmt an sich mit dem Weg
aus Wallraff's Link in #2 überein) und einmal mit dem extrahierten Code aus proj4. Beide liefern annähern gleiche Ergebnisse.

Nach meinem Verständnis müsste es eigentlich funktionieren, und wie ich lese funktioniert der Weg bei JonnyMueller scheinbar mit Fehlern <5m.

Wär wirklich nett, wenn mir hier jemand weiter helfen könnte. Es ist zwar nicht dringend - aber solche Sachen stören mich ;)
 

Wallraff

Geocacher
Hallo,

mehr aus dem Bauch heraus, der gleich eine Portion Spargel bekommt :

Es gibt in Österreich das alte "Bessel- Lambert" und das neue "WGS84-Lambert".
Hast Du schon das neue in der DB ?

Grüße Wallraff
 

lameth

Geocacher
Hallo,

nein - an sich rechne ich's (noch) in den Bessel Ellipsoid um.
Allerdings stimmen die Werte des o.a. Koordinatenpaares auch nicht mit den Angaben für "Lambert (alt)" auf geoland.at überein.

Ich hab heute mal testweise die Ellipsoidwerte in der Funktion auf die des WGS84 geändert - ohne wirklichen Erfolg.
Allerdings stimmen dann möglicherweise die Transformationsparameter für Helmert nicht mehr - kann das sein?
Lt. Angaben vom BEV sind die ja für den Bessel.

Zum aus der Haut fahren ist das :kopfwand:
 

Wallraff

Geocacher
Hallo,

ich kann nicht schlafen ...

Wenn Deine ersten beiden Koordinatenpaare zusammengehören, benutzt Du doch Lambert-neu.
Eine Umrechnung aufs Bessel-Ellipsoid ist doch gar nicht nötig. Du hast einen Punkt auf dem WGS84-Ellipsoid und willst dessen Lambert-Koordinaten auf dem WGS84-Ell. berechnen. Also kannst Du aus den geographischen Koordinaten ohne Umwege direkt die Lambert-Koordinaten berechnen ?!

Jetzt verrate mir aber mal, wo auf der neuen BEV-Seite die Transformationsparameter stehen. Die finde ich nicht mehr :???:

Grüße Wallraff
 

lameth

Geocacher
Hallo Wallraff,

nun, meinem - zugegebenermaßen bescheidenen - Verständnis der Materie nach müsste es 2 Szenarien geben:

Gegeben Koordinatenpaar 14.0241486° / 48.4388267°
Szen. 1.) Helmert Transformation & Lambert Projektion: liefert Lambert (alt)
Szen. 2.) Lambert Projektion alleine: liefert Lambert (neu)

Das würde zwar erklären, warum die Funktion ein besseres Ergebnis liefert, wenn ich die Helmert Transformation ausklammere,
allerdings sind ~50m Differenz auf beiden Achsen immer noch zu viel.
Was mMn. bedeuten würde, dass der Fehler im Projektionsalgorithmus liegt.
Der allerdings stimmt mit diesen:
-) http://surveying.wb.psu.edu/psu-surv/Projects/PASingleZone.pdf
-) http://www.rschlichte.de/geogrid/coord_conv.html
-) proj4 - http://proj.maptools.org/
überein.

Im schlimmsten Falle könnte ich sicher damit leben, da bei Abstandsberechnungen von mehreren km die paar m wohl weniger ins Gewicht fallen werden.
Aber ganz glücklich wär ich nicht damit.
 

Wallraff

Geocacher
Hallo,

alles korrekt. Ich wollte nur nicht, dass Du einen blöden Denkfehler machst - wie er auch mir passieren kann.

Hm, soll ich Dich weiter in den Wahnsinn treiben ?
Die "reguläre" Lambert-Projektion wird mit den 46. und 49. Breitenkreis gerechnet, s.a.
BEV

Ich fand schon früher mal einen
Link zur österreichischen ÜKÖ500
Dort werden auf S. 8 der 43. und 48. Breitenkreis genannt ? Arbeitest Du denen zu ? Oder haben die sich vertippt ?

Sonst weiß ich keinen Rat mehr.

@KajakFun:
Jaja, ich wollte aber wissen, wie man sich durchklicken muss.
Es geht u.a. mit Support-Download-Formatbeschreibungen und sonstige Informationen - was sofort einleuchtet ...


Grüße Wallraff
 

lameth

Geocacher
Hallo,
danke, dass du dich so um mein Problem kümmerst ;)

Also ich rechne an sich mit den Breiten 46/49.
Aber ich poste mal gerne den Code der Projektionsberechnung, möglicherweise fällt ja jemandem etwas auf:

Code:
/* Ellipsoid Parameters */
declare @a as decimal(20,10)							/* Semi-major axis of ellipsoid, in meters */
declare @f as decimal(20,10)								/* Flattening of ellipsoid */

/*
/* Ellipsoid Parameters, WGS 84  */
set @a = 6378137.0
set @f = 1 / 298.257223563
*/

/* Ellipsoid Parameters, Bessel  */
set @a = 6377397.155
set @f = 1 / 299.1528128


 /* Lambert Conformal Conic projection Parameters */
declare @Origin_Lat as decimal(20,10)				/* Latitude of origin, in radians */
declare @Origin_Long as decimal(20,10)				/* Longitude of origin, in radians */
declare @Std_Parallel_1 as decimal(20,10)			/* Lower std. parallel, in radians */
declare @Std_Parallel_2 as decimal(20,10)			/* Upper std. parallel, in radians */
declare @False_Northing as decimal(20,10)			/* False northing, in meters */
declare @False_Easting as decimal(20,10)			/* False easting, in meters */
declare @Scale_Factor as decimal(20,10)			/* Scale Factor */

set @Origin_Lat = radians(47.5)
set @Origin_Long = radians(13.333333)
set @Std_Parallel_1 = radians(46.0)
set @Std_Parallel_2 = radians(49.0)
set @False_Northing = 400000.0
set @False_Easting = 400000.0
set @Scale_Factor = 1.0

declare @e as decimal(20,10)							/* Eccentricity of ellipsoid */
declare @m1 as decimal(20,10)
declare @m2 as decimal(20,10)
declare @t0 as decimal(20,10)
declare @t1 as decimal(20,10)
declare @t2 as decimal(20,10)
declare @n as decimal(20,10)							/* Ratio of angle between meridians*/
declare @F0 as decimal(20,10)
declare @Rb as decimal(20,10)
declare @t as decimal(20,10)
declare @m as decimal(20,10)
declare @R as decimal(20,10)
declare @gamma as decimal(20,10)

/* Zone Konstants */
set @e = sqrt(2.0 * @f - square(@f))
set @m1 = cos(@Std_Parallel_1) / sqrt(1.0 - square(@e) * square(sin(@Std_Parallel_1)))
set @m2 = cos(@Std_Parallel_2) / sqrt(1.0 - square(@e) * square(sin(@Std_Parallel_2)))
set @t0 = tan(pi()/4 - @Origin_Lat / 2) / power((1.0 - (@e * sin(@Origin_Lat))) / (1.0 + (@e * sin(@Origin_Lat))), @e/2)
set @t1 = tan(pi()/4 - @Std_Parallel_1 / 2) / power((1.0 - (@e * sin(@Std_Parallel_1))) / (1.0 + (@e * sin(@Std_Parallel_1))), @e/2)
set @t2 = tan(pi()/4 - @Std_Parallel_2 / 2) / power((1.0 - (@e * sin(@Std_Parallel_2))) / (1.0 + (@e * sin(@Std_Parallel_2))), @e/2)
set @n = (log(@m1) - log(@m2)) / (log(@t1) - log(@t2))
set @F0 = @m1 / (@n * power(@t1, @n))
set @Rb = @a * @F0 * power(@t0, @n)

/* 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)

/* Solution */
set @gamma = @n * (@LONG - @Origin_Long)
set @Easting = @R * sin(@gamma) + @False_Easting
set @Northing = @Rb - @R * cos(@gamma) + @False_Northing

--select E = @Easting, N = @Northing
@LAT / @LONG sind die Inputparameter - werden natürlich Anfangs in Rad umgerechnet.
 

KajakFun

Geowizard
Hi lameth,
ich bin ja nicht der absolute Geo- und Programmierfreak aber mir ist aufgefallen,
dass im Listing z.Teil nur halbe Winkel auftauchen (wenn ich das mit den Formeln in dem geposteten Link vergleiche).
z.B.
set @t1 = tan(pi()/4 - @Std_Parallel_1 / 2 .....
set @t2 = tan(pi()/4 - @Std_Parallel_2 / 2....
set @t = tan(pi()/4 - @LAT / 2......
ist das so korrekt?
Ausserdem vermisse ich eine Iterationsschleife
Ich nehme an dass man mit den Hilfswerten und Startwerten für eine Näherungslösung eine Iterationsschleife 5-6 mal durchrattern lassen muss um exakte Werte zu erhalten?
 

lameth

Geocacher
Hallo KajakFun,

ja, die halben Winkel tauchen aus folgendem Grund auf:
In 4 Quellen für den Lambert Algorithmus hab ichs 2x mit "phi/2" und 2x mit "phi" gefunden.
Bei ganzen Winkeln "phi" ergibt sich dann leider:
-) die Werte von t/t0/t1/t2 werden negativ
-) und daraus gefolgt ergeben n, F0 und Rb komplexe Zahlen

Rechnet man mit Absolutwerten von t/t0/t1/t2 wird leider der Hochwert (Northing) negativ,
da log(t1) > log(t2) wird und (da beide negativ sind) n und alles daraus folgende tief ins negative rutscht.

Habs dann verworfen, und bin bei phi/2 geblieben ;)

Die Iterationen tauchen meines Wissens nur bei Berechnung von Winkeln auf, sprich invers bei der Berechnung von lat/long aus Lambert Rechts/Hoch.

Trotzdem Danke :)
 

Die Baumanns

Geowizard
Zur Kontrolle kannst Du es mal mit TransDat versuchen.
Der ist ziemlich gut und kann auch die Systeme in östereich.
Da kann man dann die Berechnung mal gegenchecken.

Gruß Guido
 

Wallraff

Geocacher
Hallo,

nächste Möglichkeit:
da e im Exponent auftaucht und was weiß ich anrichtet, das dazu benutzte f aufpeppen ?

1/f = 299,15281535

Ich weiß, wo mein alter QBasic-Interpreter schlummert ... aber auch erst über alle Fallstricke (Phi,Phi/2 ...) stolpern :/


Grüße Wallraff
 
Oben