Kartografie

Einfache Rechnungen

Auf dieser Seite finden sie einige Algorithmen aus Kartografie und Orientierung. Die Auswahl bevorzugt einfache Verfahren, die auch für Amateur-Programme oder sogar im Gelände in Frage kommen. Die Algorithmen dieser Seite werden in VBA-Syntax angegeben, weil Rechnungen häufig mit OpenOffice oder Excel durchgeführt werden. Sie sollten sich ohne Problem in andere Programmiersprachen übertragen lassen.
Alle Angaben sind unverbindlich, die Verwendung erfolgt auf eigenes Risiko !
Kartografie Orientierung auf der Erdoberfläche
Länge & Breite
Geografische Länge und Breite für große Maßstäbe
UTM
 Universal Transverse Mercator - einfach & schnell für kleine Maßstäbe
Grad und Bogen Grad, Minuten, Sekunden, Bogenmaß, Pi, ...
Maßstab Karte und Realität, Meilen
Ebene & Kugel Rechnung mit ebenen (UTM) oder Kugel-Koordinaten (Länge & Breite)
Abstand (Ebene) Distanz zwischen 2 Punkten einer Ebene (UTM-Koordinaten)
Abstand (Kugel) Distanz zwischen 2 Punkten einer Kugelfläche (geogr. Länge & Breite)
Richtung (Ebene) Visierlinie (Horizontalrichtung) von einem Standpunkt zu einem Zielpunkt
Richtung (Kugel) Horizontalrichtung (Azimuth) aus geogr. Länge & Breite
Himmelsrichtung Praxis-freundliche Formulierung der Richtung N-O-S-W
Bewegung Reise von einem Startpunkt um eine Strecke in eine Richtung
Höhenwinkel In der Entfernung versinken Schiffe, Türme und Gipfel
Koordinaten Verwaltung in Datenbanken, Transport-Formen
Erd-Modelle Kugel, Ellipsoid, Standard-Modelle
Verwandte Themen UTM-System, österreichische Karten, Visual Basic (VBA)
Links Ausgewählte Links zum Thema 'Kartografie'

Grad, Minuten, Sekunden - Bogenmaß - Pi (π)
Grad, Minuten, Sekunden
Geografische Winkel wie Länge (Abstand vom Null-Meridian) und Breite (Abstand vom Äquator) werden oft in Grad angegeben.
Ein voller Kreis hat 360 Grad (°), die meist als ganze Zahl angegeben werden. Kleinere Teile werden in Minuten (') und Sekunden (") angegeben, meist ebenfalls als ganze Zahlen.
1 ° = 60'1 ' = 60 "
Bei noch höherer Genauigkeit kann man die Sekunden als Gleitkommazahl angeben. In diesem Fall wird ein Winkel jedoch besser in Dezimalgrad (s.rechts) bezeichnet.
Dezimalgrad
Winkel können auch mit einer einzigen Dezimalzahl angegeben werden. Ein voller Kreis hat weiterhin 360°, allerdings werden die Grade als Gleitkommazahl bezeichnet, und zwar ohne Minuten oder Sekunden.
dezgrad = grad + (minuten + sekunden/60) / 60
Ein allfälliges Vorzeichen wird meist nur für die Grad-Zahl angegeben.

Beispiel:
grad=12; minuten=34; sekunden=56;
dezgrad = 12.58222
Die Umrechnung von Dezimalgrad in Grad, Minuten und Sekunden erfolgt am besten in mehreren Schritten (Algorithmus statt Formel).

Vor der Berechnung muss die Genauigkeit festgelegt werden:
• Wenn nur ganzzahlige Grad berechnert werden, dann genügt es, Dezimalgrad zu runden.
• Wenn Grad ... Sekunden berechnet werden, dann werden die Grad und Minuten abgeschnitten (Funktion int ), die Sekunden gerundet.
Im Beispiel rechts werden ganzzahlige Grad, Minuten und Sekunden berechnet. Der Rundungsfehler bei Berechnung ganzzahliger Sekunden beträgt maximal 0.5" = 0.0001389° (In Mitteleuropa 10-15m).

Sonderfälle wie 60' oder 60" müssen abgefangen werden. Werte wie dgrad=12.1 ergeben sonst Resultate wie 12°5'60". Die Funktion korrigiert diesen Fall auf 12°6'

Negative Winkel werden vor der Rechnung in den Absolutwert umgewandelt. Das negative Vorzeichen wird zuletzt (nur für Grad) wieder eingefügt. VBA kann auch -0 zurückgeben. Wenn das nicht möglich ist, dann muss das Vorzeichen zusätzlich im Array zurückgegeben werden.

Die berechneten {Grad, Minuten, Sekunden} werden zweckmäßig als Array zurückgegeben. Wenn die Funktion einer Programmiersprache (z.B. VBA) kein Array zurückgeben kann, muss sie mehrfach ausgeführt werden und kann jeweils nur eines der Resultate zurückgeben: Argument mode wählt einen der 3 berechneten Werte zur Rückgabe aus:
mode=1 (Standardwert) gibt Grad zurück, mode=2 die Minuten und mode=3 die Sekunden. 		Zerlegung von DezimalGrad mit VBA:
Function degdec_to_dms(deg_dec, mode)
Dim d, m, s, v As Integer
Dim deg, rest As Double
deg = deg_dec
v = Sgn(deg)
deg = Abs(deg)
d = Int(deg)
rest = (deg - CDbl(d)) * 60
m = Int(rest)
rest = (rest - CDbl(m)) * 60
s = Round(rest, 0)
If s = 60 Then
s = 0
m = m + 1
End If
If m = 60 Then
m = 0
d = d + 1
End If If v < 0 Then
d = -d
End If
' Rückgabe
If mode = 2 Then
degdec_to_dms = m
ElseIf mode = 3 Then
degdec_to_dms = s
Else
degdec_to_dms = d
End If
End Function

♦ Details zu Funktionen ganzer Zahlen (Abschneiden, Runden, Rest, ...)
▶ Bei der Verarbeitung negativer Winkel ist besondere Vorsicht angebracht - In den hier gezeigten Beispielen wird das Vorzeichen nur auf den Grad-Wert angewendet.
Programme und Programmiersprachen reagieren unterschiedlich. Insbsondere der Bereich -1°...0° ist anfällig gegen Fehler ! 		Beispiel: Zusammensetzung der Bestandteile zu DezimalGrad in Excel
=(ABS(xgrad)+(xmin+xsec/60)/60)*VORZEICHEN(xgrad)
(Ersetzen sie die Namen xgrad, xmin, xsec) durch die Namen oder Adressen der entsprechenden Zellen)
Die Funktion rechts zeigt die Zusammensetzung von Grad, Minuten und Sekunden zu Dezimalgrad.

VBA kennt im Gegensatz zu Excel nicht den Wert -0 (für Grad)
Daher wurde die Möglichkeit eingebaut, für den Bereich -1...0° ein optionales Vorzeichen (-1 oder +1) anzugeben. 		Zusammensetzung von Dezimalgrad in VBA:
Function dms_to_deg(deg, min, sec, Optional vz As Integer = 0) As Double
Dim v As Integer
Dim d As Double
If vz <> 0 Then
v = Sgn(vz)
Else
v = Sgn(deg)
End If
d = (Abs(deg) + (min + sec / 60) / 60) * v
dms_to_deg = d
End Function
Bogenmaß
In allen gängigen Programmiersprachen, in Kalkulationsprogrammen (OpenOffice, Excel, .. ) usw. werden Winkelfunktionen nicht im Gradmaß sondern im Bogenmaß berechnet.
Im Bogenmaß wird ein voller Kreis mit dem Zahlenwert 2π angegeben, für π=3.14159...
In der Regel werden alle Winkel für Rechnungen in das Bogenmaß umgewandelt, für die Ein- und Ausgabe (user interface) in das Gradmaß. Allenfalls wird das Gradmaß (Dezimalgrad) zusätzlich in Grad, Minuten und Sekunden umgewandelt (s.o.) und umgekehrt.
In Kalkulationsprogrammen (OpenOffice, Excel, ..) sollten sie jedoch Minuten und Sekunden nicht berechnen, sondern nur das Anzeige-Format ändern ! 		Umwandlung Dezimalgrad→Bogenmaß und Bogenmaß→Dezimalgrad in VBA:
Const pi = 3.14159265358979

Function grad_to_bog(grad As Double) As Double
grad_to_bog = grad / 180 * pi
End Function

Function bog_to_grad(bog As Double) As Double
bog_to_grad = bog / pi * 180
End Function
Die Kreiszahl Pi (Archimedes-Konstante, Ludolph-Zahl, π)
Pi (π):
Wenn eine Programmiersprache (VBA) die Konstante π nicht kennt, dann wird sie am Anfang des VBA-Moduls als Konstante definiert:
Const pi = 3.14159265358979
Wenn ein Programm mit größerer Genauigkeit rechnen kann:
π = 3.1415926535897932384626438328

● Achtung - Die meisten Beispiele dieser Seite verwenden die globale Konstante oder Variable pi ! 		Alternativ kann man π mit Hilfe der Arcustangens-Funktion berechnen:
pi = 4 * arctan(1)
Damit man π (z.B. in VBA) nicht jedesmal neu berechnen muss, wird diese Variable am Anfang des VBA-Moduls global definiert und danach nur bei der ersten Verwendung berechnet, z.B.
Dim pi As Double

Function grad_to_bog(grad As Double) As Double
If pi = 0 Then pi = 4 * Atn(1)
grad_to_bog = grad / 180 * pi
End Function
Neugrad
Die französische Revolution gab den Anstoß zur Durchsetzung vieler neuer Ideen, insbesondere auch des Dezimalsystems.
Ein derartiger Denkansatz sah vor, den Winkel eines vollen Kreises in 400 Teile (Neugrad, gon) zu teilen. Bis auf wenige Ausnahmen hat sich dieses System jedoch nicht durchgesetzt.
Landkarten mit Skalierung in Neugrad verwenden meist Paris und nicht Greenwich als Bezugs-Meridian. 		Unterteilung:
1 Neugrad (gon) = 100 Neuminuten (c) = 1000 Milligon (mgon)
Umrechnung:
gon = dezimalgrad / 360 * 400
1 gon = pi / 200
dezimalgrad = gon / 400 * 360

Maßstab:   Karte und Realität, Meilen
Der Maßstab gibt das Verhältnis an, in welchem die realen Maße auf die Karte übersetzt sind. Auf jeder Landkarte ist der Maßstab angeführt. Gute Karten haben bei der Legende auch ein kleines Lineal aufgedruckt.

Praxis-taugliche Landkarten haben z.B. die Maßstäbe 1:200000 (Planung, Sichtweite, Fahrrad, ...) oder 1:50000 (Detail, Wandern, ...).
Maßstab 1:200000
1cm auf der Karte ↔ 2km in der Realität.
Der Raster der UTM-Karten (5cm) entspricht 10km Realität.

Maßstab 1:50000
1cm auf der Karte ↔ 500m in der Realität.
Der Raster der UTM-Karten (2cm) entspricht 1km Realität.
Lineal und Planteiler
Maße werden von Landkarten mit Linealen oder mit Planteilern abgenommen.
Planteiler sind Plättchen aus transparentem Kunststoff mit eingeritzten Maßlinien. Sie sind zusammen mit Landkarten erhältlich oder können im Internet bestellt werden. Planteiler kann man auch leicht selbst herstellen: Muster am PC zeichnen, und auf Transparent-Folie drucken oder kopieren.
Wenn man im Gelände kein Lineal hat, dann lässt sich mit Bleistift und Papier (Rückseite einer Karte) leicht eines herstellen:
Papier scharf falten (ergibt eine gerade Linie), an den Karten-Maßstab anlegen und mit Bleistift markieren.
Karten-Maße
Die meisten älteren Karten sind entlang von Meridianen und Breitenkreisen geschnitten. Die Breitenkreise werden vom Äquator zu den Polen immer kleiner: Daher sind solche Karten nicht rechteckig sondern leicht trapezförmig. Die kürzere Seite (in Europa Nord = oben) weist zum Pol.

Berechnung der S-N-Ausdehnung ('Höhe') in der Realität und auf der Landkarte:
y = r_pol * (lat_n - lat_s) * pi / 180
ymm = y / maßstab * 1000

Beispiel: Karte ÖK200 BMN67 Graz:
Geogr.Breite Süd (unten) und Nord (oben): lat_s=46°30', lat_n=47°30'
Erdradius (polar): r_pol=6356752m
Maßstab = 1:200000
y = 6356752m * (47.5 - 46.5) * pi / 180 = 110946m = 110.95km
ymm = 110946m / 200000 * 1000 = 554.7mm
Ergebnis: Die Karte ist 555mm hoch und zeigt ein Gebiet von 111km S-N-Ausdehnung.
Berechnung der W-O-Ausdehnung ('Breite') einer Karte am südlichen (unteren) Rand:
x = r_eq * cos(lat*pi/180) * (lon_e - lon_w) * pi / 180
xmm = x / maßstab * 1000

Beispiel: Karte ÖK200 BMN67 Graz:
Breite unten: lat_s=46°30'
Geogr.Länge West (links) und Ost (rechts): lon_w=14°50', lon_e=15°50'
Erdradius (äquatorial): r_eq=6378137m
Maßstab = 1:200000
x = 6378137m * cos(46.5 * pi / 180) * (15.83 - 14.83) * pi / 180 = 76627m = 76.6km
xmm = 76627 / 200000 * 1000 = 383.1mm
Ergebnis: Die Karte ist 383mm breit und zeigt ein Gebiet von 77km W-O-Ausdehnung.
Am nördlichen (oberen) Rand ändert sich nur
Breite 'oben' lat_n=47°30'
daraus ergibt sich xmm = 376.0mm
Auf neueren UTM-Karten ist ein Raster angegeben und direkt in Meter oder km bezeichnet.
Die Angabe 598000mE bedeutet UTM-Ostwert=598km, die Angabe 5340000mN UTM-Nordwert=5340km.
Die Umrechnung ist einfach:
mm_karte = m_real / maßstab * 1000
m_real = mm_karte * maßstab / 1000
UTM-Karten sind meist nicht durch ihre Ränder sondern durch ihren Mittelpunkt definiert. Auf die exakt rechteckige Kartenfläche wird das Bild projiziert. Breitenkreise liegen zueinander parallel, jedoch nicht immer parallel zum oberen / unteren Kartenrand. Meridiane laufen zum Pol (in Europa nach Norden) trapezförmig aufeinander zu.
Eine Abschätzung der Karten-Abmessungen lässt sich jedoch mit den oben angegebenen Algorithmen immer ausführen.
Meilen, Fuß & Zoll
Der Name 'Meile' stammt vom lateinischen Wort für 1000, entsprechend der Entfernung von 1000 Doppelschritten, ca. 1480m.
Die historische Entwicklung brachte insbesondere im kleinräumigen Europa viele verschiedene lokale Meilen-Definitionen zwischen 1.5km und 11km hervor. Meist wurden dazu Unter-Einheiten im Hexagesimalsystem verwendet, z.B. Klafter, Fuß, Zoll ..
1875 wurde dieses Chaos erfolgreich durch den einheitlichen internationalen Standard 'Meter' mit dezimaler Unterteilung abgelöst.
Lediglich im US-Einflußgebiet halten sich hartnäckig Meilen, Fuß und Zoll. Das betrifft die meisten US-Landkarten ebenso wie wirtschaftliche Einflußzonen. In der US-Raumfahrt führte die Verwechslung von km und Meilen bereits zu teuren Verlusten.
Der Meilen-Maßstab auf US-Landkarten ist nicht dezimal unterteilt, sondern binär, d.h. 1/2, 1/4, 1/8, usw.
Nautische Meilen (historisch zahlreiche Varianten) sind meist durch 1 Minute des Erdumfangs definiert, d.h. ca.
1 histor.Seemeile =ca. 6378000m * 2 * π / 360 / 60 = 1855m
Diese Definitionen sind jedoch an den uneinheitlichen Erdradius gebunden und daher überholt.
Heute ist die Internationale Seemeile so definiert:
1 Seemeile = 1 sm = 1852m
1 km = 0,539956803455723 sm
An der Seemeile orientiert sich auch die Angabe der Geschwindigkeit in Knoten.
1 Knoten = 1kn = 1 sm/h = 1.852 km/h
Die (Wasser)-Tiefe wird in Faden gemessen:
1 Faden = 1/1000 sm = 1.852 m
Erst 1959 haben sich jene Staaten, die immer noch Meilen verwenden, wenigstens auf eine einheitliche Definition geeinigt:
Britische Landmeile (Statute Mile):
1 mile = 8 furlongs = 1609.344 m
1 furlong = 10 chains = 201.168 m
1 chain = 4 rods = 20.1168 m
1 rod = 5.5 yards = 5.0292 m
1 yard = 3 foot = 0.9144 m
1 foot = 12 inch = 0.3048 m
1 inch (Zoll) = 25.4 mm

1 km = 0.621371192237334 Meilen
1 cm = 0.393700787401575 Zoll
Angeblich soll die US-Landvermessung eine eigene, um ca. 3mm längere Meile und entsprechende Unter-Einheiten verwenden.
Als absurdes Kuriosum bieten ausgerechnet moderne UTM-Landkarten wieder (zusätzlich) Maßstäbe und Randmarken in Meilen an.

Kartographische Rechnungen
Hinweise zur (vorsichtigen !) Anwendung von Formeln und Algorithmen:

♣ In der Kartografie ist der Nullpunkt von Richtungs-Winkeln immer Nord, der Winkel wird im Uhrzeiger-Sinn gemessen (N→O→S→W).
In der Geometrie wird der Winkel-Nullpunkt meist rechts angenommen und entgegen dem Uhrzeigersinn gemessen.
In der Bild-Auswertung wird die Y-Koordinate meist von oben nach unten angegeben.

♣ Die geogr. Länge wird sowohl im Bereich 0..360° angegeben als auch im Bereich -180°...+180° östlich Greenwich. In den USA wird oft ohne Hinweis westliche Länge verwendet, d.h. umgekehrtes Vorzeichen der geogr. Länge. Viele Algorithmen verarbeiten die Länge nur in einem begrenzten Bereich korrekt.

♣ Viele kartographische Algorithmen sind leider für einen bestimmten Kulturkreis vereinfacht, z.B. für Europa (nur Nordhalbkugel), die USA (nur westliche Länge) oder Australien (nur Südhalbkugel). Verlässliche Formeln findet man u.a. in der Hochsee-Navigation, die nicht auf bestimmte Regionen begrenzt ist.
♣ In Kartografie und Orientierung braucht man oft die Tangens-Funktion. Dabei muss man die Spezialfälle der Winkel 0°, 90°, 180° und 270° abfangen und getrennt berechnen.
Die Tangens-Funktion kann entgegengesetzte Richtungen nicht unterscheiden. Das ist natürlich in der Praxis nicht tragbar, deshalb muss jeder Algorithmus dieses Problem lösen (Berechnung der ↓ Richtung).

♣ Als ↑ Längen-Einheit werden in Europa ausnahmslos Meter verwendet. In den USA und in ihrem Einflussgebiet findet man nur Meilen, in der Luft- und Seefahrt Seemeilen. US-Landkarten sind in Zoll (inch) geschnitten, der Maßstab meist nicht dezimal sondern binär (1/2, 1/4, 1/8...) unterteilt.

♣ Bei der Berechnung bzw. Auswertung von vertikalen Sichtwinkeln (↓ Sichthöhe) muss man schon bei Entfernungen von wenigen km die Erdkrümmung berücksichtigen.
UTM-System
Das UTM-System ist eine kartographische Projektion, bei welcher die Fläche eines Meridian-Streifens ("Orangenspalte") in eine ebene Fläche umgerechnet wird.
UTM-Koordinaten können wie cartesische Koordinaten (X,Y) verwendet werden.

Der große Vorteil des modernen UTM-Koordinatensystems zeigt sich, wenn man damit rechnen will:
Die Algorithmen sind einfach und rasch auszuführen.

GPS, digitale Landkarten und die neuen österr. Karten verwenden UTM bzw. bieten UTM als Option.
♣ Tipp: Alle österreichischen Landkarten gibt es digitalisiert auf CD (Verlag BEV) und im Internet: Daraus kann man die Koordinaten aller Punkte in Östrerreich entnehmen.
Auch aus den meisten anderen digitalen Karten (z.B. GoogleMap) kann man Koordinaten entnehmen.
Die UTM-Daten eines Punktes der Erdoberfläche bestehen aus UTM-Zone, UTM-Ost-Wert (Easting) und UTM-Nord_Wert (Northing).

Die UTM-Zone ist eine ganze Zahl im Bereich 1...60. Österreich liegt in den UTM-Zonen 32 (West) und 33 (Mitte, Ost).
Ost-Wert (Abstand vom Zentral-Meridian) und Nord-Wert (Abstand vom Äquator) werden in Metern angegeben.

♦ Details zum UTM-Koordinatensystem und zur Transformation von Koordinaten (Geogr. ↔ UTM).
Die hier gezeigten UTM-Algorithmen können problemlos (ohne Gewähr) für eigene Programme verwendet werden. Einzige Bedingung ist, dass alle verwendeten Koordinaten aus einer einzigen UTM-Zone stammen.
Fehlende Formeln kann man selbst relativ einfach aus der Analytischen Geometrie entnehmen, sie müssen jedoch für die Regeln der Kartographie angepasst werden. 		Schwerpunkt der UTM-Anwendungen ist der 'kleine' Maßstab, der zutrifft, wenn alle verwendeten Punkte in einer UTM-Zone oder zumindest in 2 benachbarten Zonen liegen. Das betrifft in Mitteleuropa ein Gebiet von einigen 100km Ausdehnung.
Kugel-Oberfläche
Wenn die verwendeten Punkte nicht aus der gleichen UTM-Zone stammen (große West-Ost-Entfernung), dann muss man die wesentlich aufwändigeren Formeln der Kugel-Oberfläche verwenden.
Vorteil: Diese Algorithmen sind allgemeiner, d.h. sie können auf beliebige Punkte der Erd-Oberrfläche angewendet werden.
Nachteil: Die Rechnungen sind im Gelände ohne Hilfsmittel nur schwer auszuführen.
Rotations-Ellipsoid
Wenn das Kugel-Modell der Erde nicht genau genug ist, muss man ein Ellipsoid verwenden. Die Algorithmen sind um eine weitere Stufe komplizierter. Sie enthalten das Kugel-Modell als Sonderfall.
Professionelle Programme erlauben die Wahl eines von mehreren standardisierten Erd-Modellen, sowie die Transformation von Koordinaten von einem Modell auf ein anderes.

Abstand in der Ebene (UTM)
Funktion abstand berechnet den Abstand zwischen 2 Punkten auf einer (flachen) Ebene. Sie eignet sich für kurze Abstände (Sichtweite, Wandern, Fahrrad, ...).
• Die beiden Punkte P1(x1,y1) und P2(x2,y2) sind durch ihre xy-Koordinaten gegeben, ihre Reihenfolge kann vertauscht werden.
• Modernen Landkarten oder GPS-Geräte liefern solche Koordinaten direkt in Meter als x=UTM-Ost (Easting) und y=UTM-Nord (Northing).
• Das Ergebnis der Funktion verwendet die gleiche Längen-Einheit wie die XY-Koordinaten, also am besten Meter (m). 		Abstand zweier Punkte auf einer Fläche mit VBA:
Function abstand(x1 As Double, y1 As Double, x2 As Double, y2 As Double) As Double
abstand = Sqr((x2 - x1)^2 + (y2 - y1)^2)
End Function
Diese einfache Rechnung kann man nur mit UTM-Koordinaten ausführen. Beide Orte müssen in der gleichen UTM-Zone liegen. Ein Taschenrechner genügt, mit etwas Geschick auch Bleistift und Papier.

Beispiel:
Auf der österr.Karte 5320-Wien kann man die UTM-Koordinaten direkt in Meter ablesen. Alle Orte der Karte liegen in der gleichen UTM-Zone 33
Korneuburg (east=x1=598945m, north=y1=5355620m)
Großenzersdorf (east=x2=615240m, north=y2=5339985)
Distanz Korneuburg - Großenzersdorf:
Die Daten werden von der Karte in Meter abgelesen und eingesetzt:
abstand = sqr((615240-598945)^2 + (5339985-5355620)^2)
= sqr(16295^2 -15635^2)
= sqr(265527025 + 244453225
= sqr(509980250)
= 22582.7m = 22.58km
▶ Wenn die beiden Orte in unterschiedlichen UTM-Zonen liegen, dann müssen die Koordinaten eines Ortes in die Zone des anderen Ortes transformiert werden. Das ist ohne Fachkenntnisse und PC kaum möglich.
• In diesem Fall verwendet man geogr.Länge und Breite (nächstes ↓ Kapitel) oder man behilft sich: Beide Karten aneinander legen (gegen Nord-Richtung etwas gedreht) und mit Lineal oder Papierstreifen messen. 		• In Österreich (allgemein in Mitteleuropa) liegen alle Orte zwischen 12° Ost und 18° Ost in UTM-Zone 33 - Das betrifft alle Orte östlich von Kundl oder Wörgl im Inntal / Tirol.
• Alle Orte zwischen 6° Ost und 12° Ost liegen in UTM-Zone 32 - Das betrifft alle österreichischen Orte in Voralrberg und in Tirol westlich von Kundl.
Wenn man keine UTM-Koordinaten hat, dann muss man die Distanz mit Funktion abstand_kugel (↓ nächstes Kapitel) aus der geogr.Länge und Breite berechnen:
Korneuburg (lon1=16.33531°, lat1=48.34602°)
Großenzersdorf (lon2=16.55106°, lat2=48.20264°)
Gerade durch 2 Punkte
Ein Vorteil von UTM ist die Möglichkeit, mit Methoden der analytischen Geometrie zu rechnen.
Berechnung einer Geraden (Steigung, Ordinaten-Abschnitt) durch 2 Punkte:
k = (y2-y1) / (x2-x1)
d = y1 - k * x1
d = y2 - k * x2

Beispiel
Gerade durch die Orte Korneuburg und Großenzersdorf:
k = (5339985-5355620) / (615240-598945) = -0.9595
d = 5355620 + 0.9595 * 598945 = 5930306

Abstand auf einer Kugel-Oberfläche
Wenn man keine UTM-Koordinaten hat, oder wenn die beiden Orte in West-Ost-Richtung zu weit entfernt sind (nicht in der gleichen UTM-Zone liegen), dann muss man die Distanz aus geogr. Länge und Breite berechnen.
Die hier gezeigten Beispiele verwenden für die Erd-Oberfläche das Modell einer Kugel. Die Algorithmen für ein Rotations-Ellipsoid sind noch aufwändiger. 		Zunächst werden einige Hilfs-Funktionen vorgestellt, danach die Anwendung zur Berechnung der Entfernung auf einer Kugel-Oberfläche.
Arcus-Sinus
VBA bietet leider keine arcsin Funktion.
Das Beispiel zeigt, wie man die fehlende Umkehrung der Sinus-Funktion durch eine eigene Funktion ersetzt. Damit keine Kollision mit der Excel-Funktion ARCSIN auftrtitt, wurde ein anderer Name gewählt.
Function my_arcsin(x As Double) As Double
Select Case x
Case Is = 1
my_arcsin = pi / 2
Case Is = -1
my_arcsin = -pi / 2
Case Else
my_arcsin = Atn(x / Sqr(-x * x + 1))
End Select
End Function
Arcus-Cosinus
VBA bietet leider keine arccos Funktion.
Das Beispiel zeigt, wie man die fehlende Umkehrung der Cosinus-Funktion durch eine eigene Funktion ersetzt. Damit keine Kollision mit der Excel-Funktion ARCCOS auftrtitt, wurde ein anderer Name gewählt.
Function my_arccos(x As Double) As Double
if x=-1 Then
my_arccos = 4 * Atn(1)
Else
my_arccos = Atn(-x / Sqr(-x * x + 1)) + 2 * Atn(1)
End If
End Function
Zentral-Winkel
Für die meisten hier vorgestellten Rechnungen auf der Kugel-Oberfläche braucht man den Zentral-Winkel: Das ist jener Winkel, unter dem man 2 Punkte der Oberfläche vom Erdmittelpunkt aus sieht.

Die vorgestellte Funktion berechnet diesen Winkel. Die Koordinaten der Punkte (geogr.Länge, Breite) werden im Bogenmaß erwartet, das Ergebnis wird im Bogenmaß zurückgegeben.

Die Argumente müssen (nur in VBA) als Werte (ByVal) übergeben werden, damit die Funktion auch VBA-intern eingesetzt werden kann. 		Berechnung des Zentral-Winkels mit VBA
Function kugel_zw(ByVal lon1 As Double, ByVal lat1 As Double, _
ByVal lon2 As Double, ByVal lat2 As Double) As Double
Dim zw As Double
If (lon1 = lon2 And lat1 = lat2) Then
kugel_zw = 0
Else
zw = Sin(lat1) * Sin(lat2)
zw = zw + Cos(lat1) * Cos(lat2) * Cos(lon2 - lon1)
zw = my_arccos(zw)
kugel_zw = zw
End If
End Function
Funktion abstand_kugel berechnet den kürzesten Abstand zwischen 2 Punkten auf der Oberfläche einer Kugel. Sie eignet sich für große Abstände in globalem Maßstab.
Die beiden Punkte P1(lon1,lat1) und P2(lon2,lat2) sind durch östl. Länge 0..2π oder -π..+π und nördl. Breite -π..+π im Bogenmaß gegeben.

• Der Verbindungs-Kreis (Orthodrome) liegt in einer Ebene mit den beiden Punkten und dem Erdmittelpunkt.
• Der mittlere Erd-Radius in Meter wird am Anfang des VBA-Moduls als Konstante definiert.
• Der Zentralwinkel zw im Bogenmaß wird mit der oben vorgestellten Hilfs-Funktion kugel_zw berechnet. 		Abstand zweier Punkte auf einer Kugel mit VBA:
Const r_erde_kugel = 6378000

Function abstand_kugel(lon1 As Double, lat1 As Double, lon2 As Double, lat2 As Double) As Double
Dim zw As Double
zw = kugel_zw(lon1, lat1, lon2, lat2)
abstand_kugel = r_erde_kugel * zw
End Function
Beispiel: Berechnung der Distanz Paris-Wien
Die Koordinaten in geogr. Länge und Breite kann man grob aus Europa-Karten entnehmen oder im Internet finden.
Paris (2.333, 48.867)
Wien (16.374, 48.209)
Vor der Verwendung von Winkelfunktionen werden die Winkel in das Bogenmaß umgerechnet:
lon1 =  2.333 / 180 * π = 0.040716
lat1 = 48.867 / 180 * π = 0.852890
lon2 = 16.374 / 180 * π = 0.285780
lat2 = 48.209 / 180 * π = 0.841406
Dieser Algorithmus ist allgemeiner anwendbar als jener für ↑ UTM, jedoch wesentlich aufwändiger zu rechnen. 		Distanz Paris - Wien
zw = arccos(sin(0.853) * sin(0.841) + cos(0.853) * cos(0.841) * cos(0.286 - 0.041) )
   = arccos(0.753 * 0.745 + 0.658 * 0.666 * 0.970)
   = arccos(0.9868365) = 0.1624347
abstand_kugel = 6378000 * 0.1624347 = 1036008m = 1036km
Der Zentralwinkel zw=0.162 beträgt im Gradmaß zw=9.3°, die Distanz 1036km.

Richtung in der Ebene (UTM)
In der Kartografie verwendet man die N-Richtung (Y-Achse auf Landkarten) als Nullpunkt der Richtungs-Skala. Die Angabe der Richtung erfolgt wahlweise in Grad (0 ... 360), als Himmelsrichtung (N→O→S→W→N), gelegentlich nach der Stundenrichtung der Uhr (0 Uhr ... 12 Uhr), für die Rechnung mit Winkelfunktionen immer im Bogenmaß (0 ... 2*π). Die Richtung wird aus den Koordinaten von 2 Punkten (Standpunkt, Zielpunkt) berechnet.
▶ Die Reihenfolge der Punkte ist entscheidend, bei ihrer Umkehrung wendet die Richtung um 180°.
Achtung: In der Geometrie gelten teilweise andere Konventionen:
Die Richtungs-Skala wird an der X-Achse orientiert
Entgegengesetzte Richtungen werden oft nicht unterschieden.
Bei der Übertragung geometrischer Algorithmen auf die Kartografie muss man daher vorsichtig vorgehen und alle möglichen Fälle sorgfältig testen.
Richtung in der Ebene (UTM)
Funktion richtung berechnet die Richtung (Kurs, heading, Azimuth) von einem Standpunkt zu einem Zielpunkt auf einer Ebene. Sie erfordert, dass beide Punkte in der gleichen UTM-Zone liegen.

Punkt P1(x1,y1) ist der Standpunkt, Punkt P2(x2,y2) der Zielpunkt. Bei Vertauschung ergibt die Funktion die umgekehrte Richtung.

Die Berechnung des rohen Richtungswinkels w (im Bogenmaß) ist einfach:
w = arctan(dx / dy)
Die Tangens-Funktion hat eine Periode von π (180°), daher müssen für dy<0 Korrekturen angebracht werden.
Außerdem muss der Sonderfall dy=0 abgefangen werden, sonst ergibt sich ein Fehler wegen Division durch 0.

Nach der Berechnung des Winkels w mit allen Korrekturen wird dieser durch Addition oder Subtraktion eines vollen Kreises (2*π) in den Bereich 0..w..2*π gebracht.

Die abgeschaltete Zeile führt optional eine kleine Korrektur der Richtung durch, die proportional zum Abstand vom Zentral-Meridian ( an dem x1=500000) ist. 		Richtungswinkel von einm Standpunkt zu einem Zielpunkt mit VBA:
Function richtung(x1 As Double, y1 As Double, _
x2 As Double, y2 As Double) As Double
Dim dx, dy, w As Double
dx = x2 - x1
dy = y2 - y1
Select Case dy
Case Is > 0
w = Atn(dx / dy)
Case Is < 0
w = Atn(dx / dy)
Select Case dx
Case Is >= 0
w = w + pi
Case Is < 0
w = w - pi
End Select
Case 0
Select Case dx
Case Is > 0
w = pi / 2
Case Is < 0
w = -pi / 2
Case 0
w = 0
End Select
End Select
While w < 0
w = w + 2 * pi
Wend
While w >= (2 * pi)
w = w - 2 * pi
Wend
' w = w + (x1 - 500000) / r_erde_kugel
richtung = w
End Function
Beispiel:
Auf der österr.Karte 5320-Wien kann man die UTM-Koordinaten direkt in Meter ablesen. Alle Orte liegen in der gleichen UTM-Zone 33
Korneuburg (east=x1=598945m, north=y1=5355620m)
Großenzersdorf (east=x2=615240m, north=y2=5339985)

Für weitere Rechnungen verwendet man den Winkel im Bogenmaß. Für die Anwendung im Gelände wandelt man den Winkel in das Gradmaß um, das man auf jedem guten Kompass direkt einstellen kann.
Alternativ oder zusätzlich berechnet man am besten noch die ↓ Himmelsrichtung. 		Richtung von Korneuburg nach Großenzersdorf:
Die Daten werden von der Karte in Meter abgelesen und eingesetzt:
dx = 615240 - 598945 = 16295
dy = 5339985 - 5355615 = -15630
' case dy<0:
w = atn(16295 / -15630) = -0.806
' case dx>=0:
w = w + pi = 2.335
w_grad = w / pi * 180 = 133.8°
Ergebnis: 134° rechts von der Nord-Richtung.
Die optionale Korrektur beträgt
dw = (x1 - 500000) / r_erde_kugel
  =  598945 - 500000 / 6378000 = 0.0155
dw_grad = 0.0155 / pi * 180 = 0.889°
Ergebnis: 135°
Genauigkeit
Die Berechnung nach der gezeigten einfachen Methode ist nur eine Näherung. Bedingt durch die Regeln der kartographischen Abbildung laufen die UTM-Koordinaten nicht genau in Süd-Nord-Richtung.
Eine genaue Berechnung aus UTM-Daten erfordert aufwändige Korrekturen. In diesem Fall verwendet man besser die Berechnung auf der ↓ Kugel-Oberfläche oder überhaupt ein Ellipsoid-Modell.
Für die Anwendung im Gelände ist der vorgestellte Algorithmus völlig ausreichend. Daten nahe dem Zentral-Meridian (CM) einer UTM-Zone sind ohnehin nicht verzerrt.
Die Verzerrung (und notwendige Korrektur) ist verkehrt proportional zum Abstand vom Zentral-Meridian.
Die ↓ Himmelsrichtung ist zwar weniger genau als ein Winkel, hat jedoch meist einen größeren praktischen Wert. Sie lässt sich auf jedem Kompass oder GPS-Gerät ablesen und aus dem Sonnenstand sogar ohne Hilfsmittel schätzen.
Außerdem kann man die Richtung auf dem Lande kaum genauer einhalten als durch die 8 wichtigsten Himmels-Richtungen angegeben. 		Himmelsrichtung aus dem Richtungs-Winkel ( ↓ VBA-Funktion nosw ):
w_grad = 134
w = 134 / 180 * pi = 2.339
i = int(2.339 * 8 / pi + 0.5) = int(6.456) = 6
richtung(6) = "SO"
Gerade durch 2 Punkte
Für die Anwendung in der analytischen Geometrie kann man die Gleichung der Geraden durch 2 Punkte berechnen (↑ Abstand in der Ebene). 		Dabei ist zu beachten, dass eine Gerade im Gegensatz zur kartograph. Richtung keine Vorzugsrichtung kennt und sich daher sowohl in die angegebene Richtung als auch in die genau umgekehrte Richtung erstreckt.

Richtung auf einer Kugel-Oberfläche
Die Richtung auf einer Kugel-Oberfläche wird genauso angegeben wie in einer Ebene - Als Differenz-Winkel zur jeweiligen lokalen Nord-Richtung.
Der Weg auf der Kugel folgt einer Orthodrome: Dieser Kreisbogen ist Teil jenes Großkreises, welcher die beiden Punkte verbindet.

Die Richtungs-Unterschiede zur ebenen Rechnung sind am Äquator minimal und an den Polen maximal. Ein bekanntes Beispiel ist die Flug-Route, welche Europa und Nordamerika verbindet. Sie verläuft wesentlich weiter nördlich als der Breitenkreis, welcher die beiden Gebiete verbindet.

Zunächst werden die Sonderfälle richtung_kugel=Ost (pi/2) und West (3*pi/2)   sowie Identität P1=P2 abgefangen.
Im allgemeinen Fall wird zuerst der Zentralwinkel zw mit der ↑ oben vorgestellten Hilfs-Funktion kugel_zw berechnet.

Daraus wird der rohe Richtungswinkel w im Bogenmaß berechnet:
w = arctan( (cos(lat1)*cos(lat2)*sin(lon2-lon1) ) / ( sin(lat2) - sin(lat1)*cos(zw) ) )
Wenn Divisionen in Algorithmen auftreten, dann ist immer Vorsicht geboten: Division durch Null muss verhindert werden.
Der Funktionsteil
eh2 = ( sin(lat2) - sin(lat1)*cos(zw) )
wird daher getrennt berechnet und die Fälle eh2=0 abgefangen. Erst danach wird durch eh2 dividiert.

Alle Winkel sind im ↑ Bogenmaß angegeben. 		Richtungswinkel auf einer Kugel-Oberfläche:
Function richtung_kugel(lon1 As Double, lat1 As Double, lon2 As Double, lat2 As Double) As Double
Dim dlon, dlat, zw, eh1, eh2, eh As Double
dlon = lon2 - lon1
dlat = lat2 - lat1
If dlat = 0 Then
If dlon > 0 Then
eh = pi / 2
ElseIf dlon < 0 Then
eh = 3 * pi / 2
Else
eh = 0
End If
Else
zw = kugel_zw(lon1, lat1, lon2, lat2)
eh1 = Cos(lat1) * Cos(lat2) * Sin(dlon)
eh2 = (Sin(lat2) - Sin(lat1) * Cos(zw))
If eh2 = 0 Then
If dlon > 0 Then
eh = pi / 2
ElseIf dlon < 0 Then
eh = 3 * pi / 2
Else
eh = 0
End If
Else
eh = Atn(eh1 / eh2)
If dlat < 0 Then
eh = eh + pi
End If
End If
End If
richtung_kugel = eh
End Function
Beispiel: Richtungswinkel Paris→Wien
Die Koordinaten in geogr. Länge und Breite kann man grob aus Europa-Karten entnehmen oder im Internet finden.
Paris (2.333, 48.867)
Wien (16.374, 48.209)
Vor der Verwendung von Winkelfunktionen werden alle Winkel in das Bogenmaß umgerechnet:
lon1 =  2.333 / 180 * π = 0.040716
lat1 = 48.867 / 180 * π = 0.852890
lon2 = 16.374 / 180 * π = 0.285780
lat2 = 48.209 / 180 * π = 0.841406

Das Beispiel rechts zeigt gerundete Zahlenwerte. Verwenden sie zur Rechnung jedoch immer die volle Genauigkeit (Typ 'Double') ohne Rundung.
Der berechnete Zentralwinkel beträgt zw=0.1624, das ist im Gradmaß zw=9.3° 		Richtung Paris→Wien
dlon = 0.286 - 0.041 = 0.245
dlat = 0.841 - 0.853 = -0.0115
zw = arccos(sin(0.853) * sin(0.841) + cos(0.853) * cos(0.841) * cos(0.286 - 0.041) )
   = arccos(0.753 * 0.746 + 0.658 * 0.666 * 0.970)
   = arccos(0.9868365) = 0.1624
eh1 = cos(0.853) * cos(0.841) * sin(0.245)
  = 0.657 * 0.666 * 0,243 = 0.106
eh2 = sin(0.841) - sin(0.853) * cos(0.162)
= 0.746 - 0.753 * 0.987 = 0.0023
eh = atn(0.106 / 0.0023) = atn(46.03) = 1.549
' case dlat<0:
eh = eh + pi = 4.691
eh_grad = eh / pi * 180 = 269°

Himmelsrichtung
Himmelsrichtung
Die Himmelsrichtung ist eine besonders einfache und Praxis-freundliche Codierung des Richtungs-Winkels. Sie kann (im Gegensatz zur Angabe in Grad) nicht missgedeutet werden, erfordert jedoch zumindest minimale Kenntnisse der Orientierung.
Funktion nosw berechnet die Himmelsrichtung als Funktion des Richtungswinkels im Bogenmaß. Falls diese Funktion nicht direkt aus dem Kalkulationsblatt aufgerufen werden soll, kann sie Private definiert werden. Das Ergebnis ist eine Zeichenkette (String).

Zunächst wird der Winkel auf den Bereich 0 <= wbog < 2*pi normiert.

Die 'Berechnung' ist ein Lookup, d.h. Verweis auf eine Werte-Liste, hier auf das Array richtung.
Der Index i des Array wird berechnet, das entsprechende Element zurückgegeben.

Wenn die Funktion oft angewendet wird, dann sollten sie das Array richtung global definieren und einmalig initialisieren. Leider scheint VBA die Array-Funktion nur für den Typ Variant zu kennen - In anderen Programmiersprachen verwenden sie natürlich die Type String für das Array richtung. 		Berechnung der Himmelsrichtung aus dem Richtungswinkel:
Option Base 0
Const pi=3.14159265358979

Function nosw(ByVal wbog As Double) As String
Dim i As Integer
Dim richtung As Variant
richtung = Array("N", "NNO", "NO", "ONO", "O", "OSO", "SO", "SSO", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW", "N")
' 0...winkel_bog...2*pi
While wbog >= (2*pi)
wbog = wbog - 2*pi
Wend
While wbog < 0
wbog = wbog + 2*pi
Wend
' Berechnung
i = Int(wbog * 8 / pi + 0.5)
nosw = CStr(richtung(i))
End Function
Beispiel für den oben berechneten Richtungs-Winkel von 134°
Die Himmelsrichtung ist zwar weniger genau als der Winkel, hat jedoch meist einen größeren praktischen Wert. Sie lässt sich auf jedem Kompass oder GPS-Gerät ablesen und aus dem Sonnenstand sogar ohne Hilfsmittel schätzen. 		Himmelsrichtung aus dem Richtungs-Winkel:
w_grad = 123
w = 134 / 180 * pi = 2.339
i = int(2.339 * 8 / pi + 0.5) = int(6.456) = 6
richtung(6) = "SO"

Fortbewegung (Richtung, Strecke)
Bewegung in der Ebene (UTM)
So wird der Zielpunkt P2(x2,y2) berechnet, wenn man sich von einem Startpunkt P1(x1,y1) auf einer Ebene in der Richtung w um eine Strecke d bewegt (Strecken in Meter, Winkel im ↑ Bogenmaß, beginnend bei Nord=0). 		Bewegung von einem Punkt P1(x1,y1) in Richtung w um eine Strecke d:
x2 = x1 + d * sin(w)
y2 = y1 + d * cos(w)
Die Bewegungs-Funktion sollte beide Koordinaten des Zielpunkts zurückgeben. Leider kann eine VBA-Funktion (zumindest an Excel) nur einen Wert zurückgeben. Daher werden hier getrennte Funktionen zur Berechnung der x- und y-Koordinate des Zielpunkts vorgestellt.
X-Koordinate (Ost-Wert) des Zielpunkts bei Bewegung um d Meter von P1(x1,y1) in Richtung w:
Function move_utm_x(x1 As Double, w As Double, d As Double) As Double
move_utm_x = x1 + d * Sin(w)
End Function
Y-Koordinate (Nord-Wert) des Zielpunkts bei Bewegung um d Meter von P1(x1,y1) in Richtung w:
Function move_utm_y(y1 As Double, w As Double, d As Double) As Double
move_utm_y = y1 + d * Cos(w)
End Function
Beispiel: (Österr.Karte 5320-Wien)
Start in Großenzersdorf, Bewegung um 14 km in Richtung 29°(NNO)
Großenzersdorf (east=x1=615240m, north=y1=5339985)
 Berechnung:
w_grad = 29
w = 29 / 180 * pi = 0.506
x2 = 615240 + 14000 * sin(0.506) = 622027
y2 = 5339985 + 14000 * cos(0.506) = 5352235
Ergebnis: Das Ziel liegt in der Ortschaft Strasshof:
P2(622027,5352235)
Gerade durch Punkt + Richtung
Ein Vorteil von UTM ist die Möglichkeit, mit Methoden der analytischen Geometrie zu rechnen.
Zur Berechnung einer Geraden (Steigung, Ordinaten-Abschnitt) wird der Winkel in das geometrische System (Nullpunkt rechts) umgerechnet:
w_geom_deg = 90 - w_karto_deg
w_geom = w_geom_deg / 180 * pi
k = tan(w_geom)
d = y - k * x

Beispiel
Gerade durch Großenzersdorf:
w_geom_deg = 90 - 29 = 61°
w_geom = 61 / 180 * π = 1.065
k = tan(1.065) = 1.804
d = 5339985 - 1.804 * 615240 = 4230063
Im Gegensatz zur kartographischen Richtung ist eine Gerade nicht gerichtet, d.h. sie erstreckt sich beiderseits des Punktes unendlich weit.
Bewegung auf einer Kugel-Oberfläche
Zunächst wird der Zentralwinkel zw (Winkelabstand) zwischen Start P1(lon1,lat1) und Ziel P2(lon2,lat2) berechnet, diesmal jedoch aus dem Verhältnis der Wegstrecke d zum Erdumfang:
zw = d / erdradius
Alle Winkel (lon,lat,zw) werden Im Bodenmaß verwendet.
Die geogr.Breite des Zielpunkts lat2 wird berechnet:
sin(lat2) = cos(w)*cos(lat1)*sin(zw) + sin(lat1)*cos(zw)
Zur Berechnung der geogr.Länge des Zielpunkts lon2 braucht man die bereits berechneten Werte von Zentralwinkel zw und Breite lat2
lon2 = lon1 + arccos((cos(zw) - sin(lat1)*sin(lat2)) / (cos(lat1)*cos(lat2)))
Einige Sonderfälle dieser Funktion müssen abgefangen werden.
Function move_sphere_lat(lon1 As Double, lat1 As Double, w As Double, d As Double) As Double
Dim zw, lon2, lat2 As Double
zw = d / r_erde_kugel
lat2 = my_arcsin(Cos(w) * Cos(lat1) * Sin(zw) + Sin(lat1) * Cos(zw))
move_sphere_lat = lat2
End Function

Die Funktion benötigt den mittleren Erdradius (als globale Konstante) und (nur in VBA) die Hilfs-Funktion my_arcsin zur Umkehr der Sinus-Funktion.

Die Berechnung der geogr.Länge(rechts) ist wesentlich komplizierter, weil man einige Sonderfälle abfangen muss. Die Berechnung von Zentralwinkel und Breite ist darin enthalten, genauso wie zur Berechnung der Breite (oben) gezeigt.
Die Hilfs-Variable vz (Vorzeichen) dient zur Unterscheidung der beiden Lösungen, d.h. zwischen 2 entgegengesetzten Änderungen der Länge..
Function move_sphere_lon(lon1 As Double, lat1 As Double, w As Double, d As Double) As Double
Dim vz As Integer
Dim zw, dlon, lat2 As Double
While w >= (2 * pi)
w = w - 2 * pi
Wend
While w < 0
w = w + 2 * pi
Wend
If (w = 0 Or w = pi) Then
move_sphere_lon = lon1
Exit Function
End If
If w <= pi Then
vz = 1
Else
vz = -1
End If
zw = d / r_erde_kugel
lat2 = my_arcsin(Cos(w) * Cos(lat1) * Sin(zw) + Sin(lat1) * Cos(zw))
dlon = vz * my_arccos((Cos(zw) - Sin(lat1) * Sin(lat2)) / (Cos(lat1) * Cos(lat2)))
move_sphere_lon = lon1 + dlon
End Function
Beispiel:
(österr.Karte 5320-Wien)
Start in Großenzersdorf, Bewegung um 14 km in Richtung 29°(NNO)
Großenzersdorf (lon=16.5510°, lat=48.2027°)

Vor der Verwendung von Winkelfunktionen werden alle Winkel in das Bogenmaß umgerechnet:
lon1 = 16.5510 / 180 * π = 0.2889
lat1 = 48.2027 / 180 * π = 0.8413
w    = 29 / 180 * π = 0.5061

Zwecks kompakter Darstellung werden im Beispiel rechts nur gerundete Werte gezeigt. Die Rechnung muss jedoch mit voller Genauigkeit ausgeführt werden, typisch mit Variablen des Typs 'Double'. 		Berechnung:
zw = 14000 / 6378000 = 0.002195
lat2 = arcsin(Cos(0.506) * Cos(0.841) * Sin(0.0022) + Sin(0.841) * Cos(0.0022))
= arcsin(0.875 * 0.666 * 0.0022 + 0.746 * 1)
= arcsin(0.747) = 0.843
lat2_deg = 0.843 / pi * 180 = 48.313°

vz = 1
dlon = 1 * arccos((Cos(0.0022) - Sin(0.841) * Sin(0.843)) / (Cos(0.841) * Cos(0.843)))
= arccos( (1 - 0.746 * 0.747) / (0.666 * 0.665) )
= arccos(0.44326371 / 0.44326428) = arccos(0.9999987) = 0.0016
lon2 = 0.2889 + 0.0016 = 0.29047
lon2_deg = 0.29047 / pi * 180 = 16.643°
Ergebnis: Das Ziel liegt in der Ortschaft Strasshof:
P2(lon=16.643°, lat=48.313°)

Höhenwinkel (Altitude)
Der Höhenwinkel (Altitude, Elevation) beschreibt den vertikalen Winkel zu einem Zielpunkt, d.h. die Sichthöhe, gemessen in Grad oberhalb oder unterhalb der Horizontalen - Diese ist übrigens nur selten identisch mit dem sichtbaren Horizont !

In einfachen Fällen wird der Höhenwinkel aus der Höhendifferenz (h_ziel - h_standpunkt) und dem Abstand d zum Zielpunkt berechnet (beide in der gleichen Einheit, normalerweise Meter):
w_alt = tan((h_ziel - h_stand)/d)
Das Ergebnis wird im Bogenmaß -pi..w..+pi geliefert.

Ab ca. 25km Entfernung muss die Erdkrümmung berücksichtigt werden: Der Fußpunkt des Ziels sinkt mit zunehmender Entfernung immer rascher unter den Horizont. Die scheinbare Höhen-Änderung beträgt
dh = r_kugel * ( cos(d/r_erde_kugel)-1 )
für r_erde_kugel = 6378000m 		Dieser Effekt macht bei 100km Entfernung bereits ca. -800m aus. Diese Sichtweite wird im Gebirge bei gutem Wetter oft erreicht.

Der Großglockner (3770m) aus der Sicht von Wien (d=300km) versinkt mit dh=-7000m bereits weit unter dem Horizont.

Ein weiterer Effekt kann am Boden vernachlässigt werden: Das Lot entfernter Punkte ändert seine Richtung relativ zum Lot am Standpunkt mit dem ↑ Zentral-Winkel. Dadurch scheinen entfernte Körper geringfügig nach hinten gekippt und damit etwas kleiner.
dh2 = - h_ziel * cos(zw)
Dieser Effekt wirkt sich erst bei extremer Sichtweite (Flugzeug, Erdumlaufbahn) merkbar aus.
Zuerst wird die Entfernung (Distanz) d in Meter berechnet.

Nächster Schritt ist die Berechnung der scheinbaren Höhen-Änderung dh in Meter.

Zuletzt wird der Höhenwinkel alt berechnet, wobei dh als Korrektur der Ziel-Höhe verwendet wird. Das Ergebnis wird im Bogenmaß ausgegeben. 		Berechnung des Höhenwinkels (Altitude) mit VBA:
Function alt_utm(x1 As Double, y1 As Double, h1 As Double, _
x2 As Double, y2 As Double, h2 As Double) As Double
Dim d, dh, alt As Double
d = Sqr((x2 - x1) ^ 2 + (y2 - y1) ^ 2)
dh = r_erde_kugel * (Cos(d / r_erde_kugel) - 1)
alt = Tan((h2 - h1 + dh) / d)
alt_utm = alt
End Function
Beispiel
In welcher Sichthöhe erscheint der Schneeberg (Wiener 'Hausberg') von Wien ?
Beide Punkte liegen in UTM-Zone 33.
Koordinaten (Ost, Nord, Höhe) in Meter:
Wien-Laaerberg (602900, 5335140, 230)
Schneeberg (560695, 5291440, 2061)

Ergebnis:
Entfernung d=61km
Korrektur für Erdkrümmung dh=-289m
Korrektur für Zentralwinkel: dh2=-8cm
Höhenwinkel alt_deg=1.45°
Richtung (↑) azi_deg=225° = "SW"
Tatsächlich ist der Schneeberg bei klarem Wetter vom Süden Wiens und von vielen anderen hoch gelegenen Orten im Osten Österreichs gut sichtbar (und leicht erkennbar). 		Sichthöhe (Altitude) des Schneebergs von Wien:
d = sqr((560695-602900)^2 + (5291440-5335140)^2)
= sqr(42205^2 + 43700^2) = sqr(3690952025) = 60753m

dh = 6378000 * (cos(60753 / 6378000) -1)
= 6378000 * (cos(0.00952) -1)
= 6378000 * (-4.54E-5) = -289m

alt = tan((2061-230-289)/60753)
= tan(1542/60753) = tan(0.025376) = 0.025381
alt_deg = 0.0254 / pi * 180 = 1.45°
Panorama:
Mit Ortskoordinaten und Höhen-Daten lässt sich das sichtbare Panorama näherungsweise recht gut berechnen (aufklären).
• Sammeln sie die Daten (UTM-Koordinaten und Höhe) aller wichtigen Punkte (z.B. Gipfel, Türme, ...).
• Berechnen sie für jeden Punkt ↑ Entfernung, (horizontale) ↑ Richtung (Azimuth) und Höhenwinkel (Altitude).
• Sortieren sie die Zielpunkte (Gipfel ..) nach fallender Entfernung.
• Erzeugen sie ein XY-Diagramm der Daten:
Die Richtung 0...360° wird in X-Richtung aufgetragen, die Sichthöhe (meist -5°...+5°) in Y-Richtung.
• Zeichnen sie an den berechneten Punkten 'Berge' ein, z.B. als Dreiecke mit gefüllter Fläche.
Bei richtiger Sortierung werden entfernte Gipfel durch nähere überdeckt, falls sie in die Sichtlinie fallen.
Camera als Meßgerät
Eine Digital-Camera ist ein sehr genaues Gerät zur Winkel-Messung:

• Die Optik übersetzt Sichtwinkel mit guter Korrelation in Bildpunkte. Pixel in X-Richtung sind stets um den gleichen Richtungs-Winkel (Azimuth) voneinander entfernt, Pixel in Y-Richtung um den gleichen Höhenwinkel (Altitude).

• Am besten verwendet man eine reproduzierbare Brennweite, z.B. die End-Position eines Zoom-Objektivs. Die Verzerrung ist bei Weitwinkel-Einstellung (kurze Brennweite) am größten, bei Tele-Einstellung (lange Brennweite) am geringsten.
Beispiele:
• Eichung der Optik: Aus den Aufnahmen bekannter Objekte (Türme, Berge, ...) von einem bekannten Standpunkt lässt sich die Auflösung der Optik (Sichtwinkel pro Pixel) berechnen.

Mit geeichter Optik lassen sich alle einschlägigen Beispiele der Trigonometrie lösen, z.B.:
• Bestimmung des Standorts aus Aufnahmen bekannter Objekte.
• Bestimmung der Koordinaten unbekannter Objekte aus mehreren Aufnahmen von bekannten Standorten (nächster Absatz).
• Berechnung der Höhe von Gebäuden aus Aufnahmen in bekannter Entfernung...
Schnittpunkt von 2 Geraden
Für einfache trigonometrische Rechnungen muss oft der Schnittpunkt von 2 Geraden bestimmt werden.
Die Gleichung einer Geraden lässt sich u.a. ↑ aus 2 Punkten oder ↑ aus Punkt + Richtung berechnen.
Schnittpunkt (x,y) von 2 Geraden (k,d):
x = (d2 - d1) / (k1 - k2)
y = k1 + x * d1
y = k2 + x * d2

Beispiel
Ein unbekannter Berg wird von den Orten Amstetten und St.Pölten fotografiert. Aus den Aufnahmen wurde die Richtung berechnet. Um welchen Berg handelt es sich ?

UTM-Koordinaten der Standpunkte (Zone 33) / Sichtwinkel in Grad:
Amstetten (490230, 5330030) / 139.5° (SO)
St.Pölten (546590, 5339275) / 220° (SW)
Berechnung der beiden Sicht-Geraden (↑ Details):
w_geom_deg_1 = 90 - 139.5 = -49.5°
w_geom_1 = -49.5 / 180 * pi = -0.864
k_1 = tan(-0.864) = -1.171
d_1 = 5330030 + 1.171 * 490230 = 5904016

w_geom_deg_2 = 90 - 220 = -130°
w_geom_2 = -130 / 180 * pi = -2.269
k_2 = tan(-2.269) = 1.192
d_2 = 5339275 - 1.192 * 546590 = 4687874
Schnittpunkt der beiden Geraden:
x = (5904016 - 4687874) / (1.192 + 1.171) = 514746
y = 1.192 * 514746 + 4687874 = 5301325
Ein Blick auf die Karte zeigt, dass es sich beim gesuchten Berg um den Ötscher (515150, 5300970) in 44km Entfernung handelt. Die Genauigkeit von +/- 500m ist realistisch.

Koordinaten speichern, transportieren und transformieren
Datenbank
Nach der Theorie dürfen nur Basis-Daten in einer Datenbank gespeichert werden (z.B. nur geogr.Länge und Breite). Wenn es Daten gibt, die daraus berechnet werden können (z.B. UTM-Koordinaten einer oder mehrerer Zonen), dann dürfen diese nicht in der Datenbank gespeichert werden. Sie müssen aus den Basis-Daten stets Live berechnet werden. Damit wird garantiert, dass keine widersprüchlichen Daten gespeichert werden, und Änderungen der Basis-Daten auf alle abgeleiteten Daten übertragen werden.
In der Praxis durchbricht man (zumindest bei Orts-Koordinaten) meist die strenge Regel: Da sich die Basis-Daten normalerweise nicht ändern, speichert man auch die abgeleiteten Daten (meist UTM-Koordinaten). Damit wird die Anwendung um Größenordnungen beschleunigt. Man muss jedoch durch eigene Programme sicherstellen, dass bei jeder Änderung von Daten alle abgeleiteten Daten sofort neu berechnet und gespeichert werden.
Transport
Für den professionellen Austausch geografischer Daten haben sich einige Standards etabliert.
Im Amateur-Bereich dominieren Datenformate diverser GPS-Geräte und -Programme. Die meisten GPS-Programme können mehrere dieser Formate lesen / schreiben / umwandeln. 		NMEA
DEM (Digital Elevation Model) - Höhen-Daten in Form eines Gitter-Rasters variabler Maschenweite.
Algorithmen
für den eigenen Grbrauch werden auf einer eigenen Seite dieses Webs vorgestellt: → Koordinaten-Transformation 		Die Genauigkeit ist begrenzt, jedoch für den Amateur-Bereich meist ausreichend. Anwendbar in Kalkulations-Programmen (OpenOffice, Excel) und Programmiersprachen.
Live-Online Transformation im Internet
Wer mit Kartografie arbeitet, muss häufig Orts-Koordinaten umrechnen. Für wenige Daten kann man einen der vielen kostenlosen Online-Umrechner im Internet benutzen.
Beispiele für Algorithmen finden sie auf der UTM-Seite.
DMAP,
Programme zur Koordinaten-Transformation:
Wenn größere Datenmengen öfters zu verarbeiten sind, dann werden dazu eigene Programme benötigt.
Das Windows-Programm Transdat vom deutschen Hersteller Killetsoft ist in professioneller Qualität als Demo-Version verfügbar.

Erd-Modelle:   Kugel und Ellipsoid
Kugel
Einige Formeln bei Verwendung des Kugel-Modells, mit üblichen Zahlenwerten.
Radiusr = 6378000 m = 6378 km
Umfangu = 2 * r * π = 40074156 m = 40074 km
Fläche f = 4 * r^2 * π = 5.112E+14 m^2 = 5.112E+08 km^2
Volumen v = 4/3 * r^3 * π = 1.087E+21 m^3 = 1.087E+12 km^3
Masse m = 5.974E+24 kg
Das Kugel-Modell ist zwar eine 'ungenaue' Vereinfachung, jedoch wesentlich rascher und einfacher zu berechnen. Es wird zur Abschätzung und immer dann verwendet, wenn keine Daten oder Algorithmen für ein Ellipsen-Modell verfügbar sind.

Großkreis
Jeder Kreis auf der Kugel-Oberfläche, dessen Mittelpunkt sich im Kugel-Mittelpunkt befindet, ist ein Großkreis (z. B. der Äquator). Alle Großkreise sind gleich lang.
Der Ausdruck wird oft auch auf Kreisbögen mit dem gleichen Mittelpunkt angewendet (z.B. auf Meridiane = Kreisbögen von Pol zu Pol). Ein Großkreis-Bogen ist die kürzeste Verbindung zwischen 2 Punkten der Kugel-Oberfläche - Deshalb fliegt man zwischen Europa und Nordamerika nicht in Ost↔West Richtung, sondern auf der kürzeren 'Polroute'.
Auf einem Ellipsoid sind die Verhältnisse schwieriger zu beschreiben, dort hängt die Bogenlänge von der Richtung ab.
Kreisbogen:
Länge eines Großkreis-Bogens (z.B. Meridian) für einen Zentralwinkel zw in Grad (0..360):
b = 2 * r * π * zw / 360
ergibt z.B. für zw=1°   b = 111317 m = 111.3 km, unabhängig von der Orientierung des Bogens für alle Großkreise (nur auf einer Kugel).

Breitenkreis - ein Kreis parallel zum Äquator:
r(bk) = r * cos(lat)
u(bk) = 2 * r * cos(lat) * pi
für die geografische Breite lat = Zentralwinkel zwischen dem Äquator und einem Punkt am Kreis im Bogenmaß. Das ergibt z.B. für 48° Breite:
lat = 48 / 180 * π = 0.838
r(bk48°) = 4267715 m   (=67% des Erdradius)

Für Punkte der Oberfläche gilt
x = r * cos(lat) * cos(lon)
y = r * cos(lat) * sin(lon)
z = r * sin (lat)
für die geografische Breite lat bzw. geografische Länge lon im Bogenmaß. Die z-Achse durch die beiden Pole ist durch die Rotation natürlich festgelegt, die x-Achse willkürlich durch den Null-Meridian (Greenwich).
Ellipsoid
Hier wird nur der Spezialfall 'Rotationsellipsoid' angeführt, mit den 3 Achsen a > b = c
● Ein Rotations-Ellipsoid hat im Gegensatz zu einer Kugel eine Vorzugs-Richtung - Die Rotationsachse mit 2 Polen.
Die Ebene normal auf die Achse sowie durch den Mittelpunkt ist die Äquator-Ebene, ihre Schnittlinie mit der Oberfläche der Äquator-Kreis (= längster Umfang). Diese Elemente werden als Fixpunkte der Orientierung eingesetzt. - Nicht nur auf der Erde sondern für alle größeren (rotierenden) Himmelskörper.

Zahlenwerte des international am häufigsten verwendeten Erdmodells WGS84 (ohne Gewähr):
Große Halbachse
semi-major axis		 a = 6378137 m
a^2 = b^2 + e^2
Kleine Halbachse
semi-minor axis		 b = a * (1-f) = sqr (a^2 - e^2) = 6356752.31424518 m
Exzentrizität
eccentricity		 e = sqr(a^2 - b^2) = 521854 m
Abflachung
flattening		 f = 1 - b/a = 0.00335281066474746
Inv.Abflachung
inverse flattening		 1 / f = 298.257223563
In Tabellen findet man immer den Wert von a
Als zweiter Wert wird meist 1/f angeführt, seltener b oder e
Aus 2 dieser Angaben lassen sich alle anderen Werte berechnen, allenfalls durch Umformen der Gleichungen.

Zentrum
Unabhängig von der Form unterscheiden sich die Erdmodelle geringfügig durch die räumliche Lage des Zentrums.
Der Unterschied (dx, dy, dz) wird meist in Metern relativ zum Modell WGS84 angegeben.
Die beiden folgenden Werte werden in Berechnungen oft verwendet und daher meist als Konstanten vorgegeben oder einmalig vor Beginn anderer Rechnungen als globale Konstanten oder Variable definiert:
e1^2 e1^2 = (a^2 - b^2) / a^2 = 2*f - f^2 = 0,00669437999014132
e2^2 e2^2 = (a^2 - b^2) / b^2 = 0,00673949674227643
Die große Halbachse a ist der Abstand eines Punkts am Äquator vom Zentrum.
Die kleine Halbachse b ist der Abstand eines Pols vom Zentrum.
Die Exzentrizität e ist der Abstand zwischen Zentrum und Brennpunkt (focus) einer Ellipse. In einem Rotations-Ellipsoid gibt es keine Brennpunkte sondern einen Kreis mit den entsprechenden Eigenschaften in der Äquator-Ebene. e ist der Radius dieses Kreises.

Kreisbogen:
Die Berechnung ist vergleichsweise aufwändig und wird mit einer Reihen-Entwicklung gelöst (s. Koordinaten-Transformation).

Breitenkreis:
r(bk) = a * cos(lat)
u(bk) = 2 * a * cos(lat) * pi
für die geografische Breite lat im Bogenmaß.
Der Radius des Breitenkreises auf lat=48° ist um 92km länger als im Kugel-Modell.

Für Punkte der Oberfläche gilt
x = a * cos(lat) * cos(lon)
y = a * cos(lat) * sin(lon)
z = b * sin(lat)
für die geografische Breite lat bzw. geografische Länge lon im Bogenmaß.
Datum-Transformation
Im Laufe der Geschichte wurde die Vermessung der Erde immer genauer. Unabhängig davon wurden die Erd-Modelle für bestimmte Staaten und Regionen optimiert. So entstanden zahlreiche standardisierte Erd-Modelle. Heute verwendet man nur noch global optimierte Modelle.

Die Umrechnung (geodätische "Datum-Transformation") von Orts-Koordinaten in ein anderes Erd-Modell ist relativ aufwändig: Die Punkte der Oberfläche werden als fix angenommen, das beschreibende Koordinaten-System (Ellipsoid) wird geringfügig verschoben und gedreht. Die Orts-Daten eines Punktes unterscheiden sich je nach Erd-Modell um einige 10-100 Meter.

Zur Umrechnung benötigt man 3 Faktoren für die Lage des Zentrums (Ursprung des Koordinatensystems) in Meter.
3 weitere Faktoren für die Orientierung der Achsen im Raum werden meist in Bogensekunden angegeben. Für die Verwendung in Winkelfunktionen werden sie in das Bogenmaß umgerechnet: 1' → 2*π/360/60/60 = 4.848E-6

Ein weiterer Faktor 'Skalierung' wurde aus praktischen Gründen zusätzlich eingeführt, obwohl das eine Redundanz darstellt. Der Skalenfaktor s ist sehr nahe 1 und wird meist als Abweichung von 1 angegeben. (z.B. ds = 1ppm = 1 parts per million = 1E-6, dann ist s=1.000001)

Das ergibt insgesamt 7 Faktoren, die bei der Helmert-Transformation vollständig berücksichtigt werden. Genaue Angaben und Online-Umrechnungen findet man im Internet.
Für rasche und weniger genaue Transformationen werden meist nur 3 Faktoren (Molodensky-Transformation) verwendet und die übrigen Faktoren vernachlässigt.
Die verwendeten Methoden sollten sich an der Genauigkeit der verfügbaren Daten orientieren. Sowohl die terrestrische Vermessung mit (Laser)-Theodoliten als auch die Vermessung mit Satelliten (GPS) bieten nur begrenzte Genauigkeit. Daten älterer Herkunft sind um Größenordnungen weniger genau.

Beispiel: Transformation ITRS ↔ ETRS89:
tx=0.054m, ty=0.051m, tz=-0.048m
rx=0.000081*δt, ry=0.00049*δt, rz=-0.000792*δt
s=0
für δt = Zeit in Jahren seit Fixierung des ETRS (1989). Diese Transformation beschreibt die Kontinentaldrift von Europa auf der Erdoberfläche.
Ältere Landkarten verwenden u.a. das Modell 'Bessel 1841' (Österreich - auch 'Hermannskogel', 'MGI'), 'Potsdam' (Deutschland), 'CH1903' (Schweiz), 'Lambert' (Frankreich) usw.
In vielen Tabellen sind Erd-Modelle durch ihre Differenzen zum derzeitigen Standard-Modell WGS84 definiert:
Differenz der großen Halbachse:
da = a(WGS84) - a(Modell)
Differenz der Abplattung:
df = f(WGS84) - f(Modell)
Die Differenz df wird manchmal (Garmin-GPS) mit dem Faktor 10000 multipliziert.
Für das Ellipsoid Bessel 1841 (Österreich) ergibt sich (unverbindlich)
dx=653, dy=-212, dz=449
da=653,135, df=0.10037483
Wenn man die Drehung der Koordinatensysteme vernachlässigt, dann kann man (GPS-Geräte) die Erdmodelle mit den hier angeführten 5 Faktoren berechnen.

Ausgewählte Links zum Thema 'Kartografie'

DMAP - UTM zones of the world, UTM Calculator (Online), Mehrere US-Universitäten hatten ausgezeichnetes Material zu Kartografie, Orientierung usw. im Internet publiziert. Leider wurde das Meiste davon wegen der Terror-Paranoia der US-Regierung gelöscht. Das erinnert an die grotesken Versuche kommunistischer Regierungen, 'Republik-Flüchtlinge' durch falsche Landkarten und Ortstafeln zu verwirren und zu behindern. Man kann nur hoffen, dass die Zukunft wieder mehr Freiheit bringt. Europäische Universitäten sind aus Furcht vor eventuellen juristischen Folgen oder aus Bequemlichkeit leider nur selten bereit, Ihre Kurse oder Daten ins Internet zu stellen.
Univ.Texas - Map Collection, Map Sites (Links !) CIA - Maps & Publications

Geoid: USA,

DEM: USA
Wikipedia Wikipedia:
Erd-Modelle: WGS84, GRS80, ETRS89,
Bildformat GeoTIFF, 		MapRef (en): Sammlung geodätischer Daten: Theorie, Projektion & Transformation, Werkzeug, nationale referenzen (at, ch, de)

Kowoma - Erdmodelle,

Zum Inhaltsverzeichnis Homepage von PS-Trainer   >   Entwicklung   >   Algorithmen   >   Kartografie