====== std.geo — Geokoordinaten ======
→ [[lyx_-_programmiersprache:units|Zurück zur Unit-Übersicht]]
Geografische Koordinaten in **Mikrograd**: 1 Grad = 1 000 000 Einheiten, gespeichert als ''int64''. Damit lässt sich ohne Fließkomma auf etwa 11 cm genau rechnen. Die Unit liest und schreibt Koordinaten, prüft ihre Gültigkeit, bildet Mittelpunkte und umschließende Rechtecke und schätzt Entfernungen.
Einsatzbereiche: Standortdaten, Umkreissuche, Geo-Fencing, Kartenausschnitte.
**Autor:** Andreas Röne\\
**Copyright:** 2024–2025 Andreas Röne\\
**Quelle:** ''std/geo.lyx''
**Seit lyxc 1.0.20F stimmen die Entfernungen** (#1481). Dieselbe Strecke Berlin–Hamburg wie in der früheren Fehlermessung:
^ Funktion ^ vorher ^ jetzt ^
| ''DistanceM'' | 437 128 m (+71 %) | **255 775 m** |
| ''HaversineDistanceM'' | 390 522 m (+53 %) | **255 775 m** |
Die Luftlinie beträgt rund 255 km — beide Verfahren liefern damit denselben, plausiblen Wert. Ursache war eine um Faktor 10 verrechnete Kosinustabelle.
===== Import =====
import std.geo;
----
===== Typen =====
==== GeoPoint ====
^ Feld ^ Typ ^ Bedeutung ^
| ''x'' | ''int64'' | **Längengrad** in Mikrograd |
| ''y'' | ''int64'' | **Breitengrad** in Mikrograd |
Die Reihenfolge ist die eines Koordinatensystems, nicht die der üblichen Schreibweise: ''x'' ist die Länge (Ost/West), ''y'' die Breite (Nord/Süd). Auch ''GeoPointNew'' nimmt **zuerst die Länge**:
var berlin: GeoPoint := GeoPointNew(13404954, 52520008); // lon, lat
Wer die vertraute Reihenfolge „52.52 N, 13.40 E" direkt überträgt, vertauscht die Felder — der Punkt landet dann im Indischen Ozean statt in Berlin. Ein Struct-Literal gibt es nicht, ''GeoPointNew'' ist der einzige Weg.
==== GeoBBox ====
^ Feld ^ Typ ^ Bedeutung ^
| ''minX'' | ''int64'' | westlichster Längengrad |
| ''minY'' | ''int64'' | südlichster Breitengrad |
| ''maxX'' | ''int64'' | östlichster Längengrad |
| ''maxY'' | ''int64'' | nördlichster Breitengrad |
Die Felder liegen flach nebeneinander statt als zwei ''GeoPoint''. Der Quelltextkommentar begründet das mit fehlender Unterstützung für verschachtelte Structs — das trifft **nicht mehr zu**: ''Circle'' in [[lyx_-_programmiersprache:units:circle|std.circle]] enthält einen ''Vec2'', und ''Range2D'' greift mit ''r2d.x.min'' über zwei Ebenen zu. Die flache Form ist damit heute nur noch eine Altlast. ''GeoBBoxGetMin'' und ''GeoBBoxGetMax'' bauen bei Bedarf einen ''GeoPoint'' daraus.
----
===== Funktionen =====
==== Lesen und Schreiben ====
^ Signatur ^ Beschreibung ^ Status ^
| ''ParseLat(s: pchar): int64'' | Dezimalgrad als Text → Mikrograd | ok |
| ''ParseLon(s: pchar): int64'' | dito | ok |
| ''FormatDecimal(coord: int64): pchar'' | Mikrograd → Text mit sechs Nachkommastellen | ok |
| ''ParseDMS(degrees, minutes, seconds, direction: int64): int64'' | Grad/Minuten/Sekunden → Mikrograd; ''direction'' ist der Zeichencode von ''N''/''S''/''E''/''W'' | ok |
| ''FormatDMS(coord: int64, isLatitude: bool): pchar'' | Mikrograd → Grad/Minuten/Sekunden | **feste Zeichenkette** (#1483) |
''ParseDMS'' rundet ganzzahlig; die Abweichung liegt bei etwa 13 Mikrograd, also gut einem Zentimeter.
==== Prüfen ====
^ Signatur ^ Beschreibung ^ Status ^
| ''IsValidLat(lat: int64): bool'' | −90…+90 Grad | ok |
| ''IsValidLon(lon: int64): bool'' | −180…+180 Grad | ok |
| ''IsValidGeoPoint(p: GeoPoint): bool'' | beide Achsen | ok |
| ''GeoPointNew(lon, lat: int64): GeoPoint'' | Punkt anlegen | ok |
==== Entfernungen ====
^ Signatur ^ Beschreibung ^ Status ^
| ''DistanceM(p1, p2: GeoPoint): int64'' | Abstand in Metern | **+71 %** (#1482) |
| ''DistanceKm(p1, p2: GeoPoint): int64'' | Abstand in Kilometern | erbt #1482 |
| ''DistanceMLegacy(lat1, lon1, lat2, lon2: int64): int64'' | Abstand aus Einzelwerten — **Breite zuerst** | erbt #1482 |
| ''DistanceMCorrected(p1, p2: GeoPoint): int64'' | mit Breitenkorrektur | **+69 %** (#1481) |
| ''HaversineDistanceM(p1, p2: GeoPoint): int64'' | Kugelentfernung | **+53 %** (#1481) |
| ''IsWithinDistanceM(p1, p2: GeoPoint, thresholdM: int64): bool'' | Abstand unter Schwelle? | ok |
| ''DistanceSq(p1, p2: GeoPoint): int64'' | Quadrat der **Mikrograd**-Differenz, nicht der Meter | ok, aber missverständlich |
| ''CorrectLongitudeForLatitude(dLon, lat: int64): int64'' | Längendifferenz auf die Breite umrechnen | ok |
''DistanceSq'' ist als reiner Vergleichswert gedacht — es zieht keine Wurzel und rechnet nicht in Meter um. Innerhalb einer engen Breitenzone taugt es zum Sortieren nach Nähe; über größere Nord-Süd-Spannen hinweg nicht, weil die Längenachse dort anders skaliert.
==== Mittelpunkte und Rechtecke ====
^ Signatur ^ Beschreibung ^ Status ^
| ''Midpoint(p1, p2: GeoPoint): GeoPoint'' | Mittelpunkt zweier Punkte | ok |
| ''MidpointLat(lat1, lat2: int64): int64'' · ''MidpointLon'' | einzelne Achse | ok |
| ''BoundingBoxFromPoints(p1, p2: GeoPoint): GeoBBox'' | kleinstes Rechteck um zwei Punkte | ok |
| ''GeoBBoxContains(bb: GeoBBox, p: GeoPoint): bool'' | Punkt im Rechteck? | ok |
| ''GeoBBoxOverlaps(a, b: GeoBBox): bool'' | Überschneidung zweier Rechtecke | ok |
| ''GeoBBoxCenter(bb: GeoBBox): GeoPoint'' | Mittelpunkt | ok |
| ''GeoBBoxGetMin(bb: GeoBBox): GeoPoint'' · ''GeoBBoxGetMax'' | Ecke als Punkt | ok |
| ''IsPointInRect(p, min, max: GeoPoint): bool'' | Punkt gegen zwei Eckpunkte | ok |
| ''IsPointInRect2(p: GeoPoint, minLat, minLon, maxLat, maxLon: int64): bool'' | dito mit Einzelwerten | ok |
| ''BoundingBoxMin''/''Max''/''Center'', ''DoBoundingBoxesOverlap'' | ältere Fassungen ohne ''GeoBBox'' | ok |
| ''CalculateBoundingBox(center: GeoPoint, radiusM: int64): GeoBBox'' | Umkreis-Rechteck | ok |
| ''CalculateBoundingBoxMin''/''Max'' | dito als Einzelpunkte | dito |
Für neuen Code die ''GeoBBox''-Fassungen verwenden — die älteren Helfer mit losen ''GeoPoint''-Paaren machen dasselbe, nur umständlicher.
==== Richtung und Verschiebung ====
^ Signatur ^ Beschreibung ^ Status ^
| ''Bearing(p1, p2: GeoPoint): int64'' | Peilung in Mikrograd, 0 = Nord | ok |
| ''AddOffsetM(center: GeoPoint, bearing, distanceM: int64): GeoPoint'' | Punkt verschieben; ''bearing'' in Mikrograd | ok |
----
===== Beispiel =====
import std.geo;
fn main(): int64 {
// Koordinaten einlesen: Text -> Mikrograd
var lat: int64 := ParseLat("52.520008");
var lon: int64 := ParseLon("13.404954");
PrintLn("Berlin: lat ", IntToStr(lat), " lon ", IntToStr(lon), " (Mikrograd)");
PrintLn("zurueck: ", FormatDecimal(lat), " / ", FormatDecimal(lon));
var berlin: GeoPoint := GeoPointNew(lon, lat);
var hamburg: GeoPoint := GeoPointNew(ParseLon("9.993682"), ParseLat("53.551086"));
PrintLn("gueltig? ", IntToStr(IsValidGeoPoint(berlin) as int64),
" Suedpol+5 Grad gueltig? ", IntToStr(IsValidLat(95000000) as int64));
// Mittelpunkt und umschliessendes Rechteck
var m: GeoPoint := Midpoint(berlin, hamburg);
PrintLn("Mittelpunkt: ", FormatDecimal(m.y), " / ", FormatDecimal(m.x));
var bb: GeoBBox := BoundingBoxFromPoints(berlin, hamburg);
PrintLn("BBox: ", FormatDecimal(bb.minY), " bis ", FormatDecimal(bb.maxY), " Nord");
PrintLn(" ", FormatDecimal(bb.minX), " bis ", FormatDecimal(bb.maxX), " Ost");
PrintLn("Berlin darin? ", IntToStr(GeoBBoxContains(bb, berlin) as int64));
PrintLn("Muenchen darin? ", IntToStr(GeoBBoxContains(bb, GeoPointNew(11576124, 48137154)) as int64));
// DMS einlesen
PrintLn("ParseDMS(52, 31, 12, 'N') = ", FormatDecimal(ParseDMS(52, 31, 12, 78)));
return 0;
}
Ausgabe:
Berlin: lat 52520008 lon 13404954 (Mikrograd)
zurueck: 52.520008 / 13.404954
gueltig? 1 Suedpol+5 Grad gueltig? 0
Mittelpunkt: 53.035547 / 11.699318
BBox: 52.520008 bis 53.551086 Nord
9.993682 bis 13.404954 Ost
Berlin darin? 1
Muenchen darin? 0
ParseDMS(52, 31, 12, 'N') = 52.520013
''ParseDMS'' bekommt die Himmelsrichtung als Zeichencode — 78 ist ''N'', 83 ''S'', 69 ''E'', 87 ''W''.
----
===== Was nicht stimmt =====
Geprüft mit ''lyxc 1.0.21A'', Referenzstrecken aus bekannten Städtekoordinaten.
==== Die Kosinustabelle ist um Faktor 10 zu schwach (#1481) ====
Ein Längengrad wird zum Pol hin kürzer — bei 50° Nord misst er 71,7 km statt 111,3 km. ''CosLatTable'' soll diesen Faktor liefern, rechnet aber ein Zehntel der vorgesehenen Korrektur:
^ Breite ^ geliefert ^ cos korrekt ^
| 0° | 1 000 000 | 1 000 000 |
| 45° | 987 728 | 707 107 |
| 60° | 978 182 | 500 000 |
| 90° | 950 910 | **0** |
Am Pol müsste die Korrektur gegen null gehen; sie liegt bei 0,95 und ist damit praktisch abgeschaltet. Der Quelltextkommentar beschreibt die richtige Formel — der Code teilt durch 165 statt durch 16 500.
Betroffen sind alle Funktionen, die auf dieser Tabelle aufsetzen: ''HaversineDistanceM'', ''DistanceMCorrected'', ''CorrectLongitudeForLatitude'', ''CalculateBoundingBox'' und ''AddOffsetM''.
Rechnet man die Haversine-Formel mit dem korrekten cos(53°) nach, ergibt sie für Berlin–Hamburg **255,7 km** — also genau den richtigen Wert. Die Formel selbst ist in Ordnung, nur die Tabelle nicht.
==== DistanceM rechnet doppelt daneben (#1482) ====
Berlin -> Hamburg = 437 128 m echte Luftlinie 255 km
Berlin -> Muenchen = 589 685 m echte Luftlinie 504 km
1 Grad Ost-West bei 50N = 111 319 m korrekt 71 700 m
1 Grad Nord-Sued = 111 319 m korrekt — hier stimmt es
Zwei Ursachen: Längendifferenzen werden wie Breitendifferenzen behandelt (keine Korrektur), und statt der Wurzel aus der Quadratsumme wird ''größerer + halber kleinerer'' gerechnet — eine Achteck-Näherung, die die Diagonale um bis zu 12 % überschätzt.
**Reine Nord-Süd-Strecken stimmen exakt.** Genau deshalb fällt der Fehler bei einem schnellen Test leicht durch.
Wer Entfernungen braucht, rechnet vorerst selbst — mit ''IntSqrt'' aus [[lyx_-_programmiersprache:units:math|std.math]], das korrekt arbeitet:
import std.math;
fn AbstandM(p1: GeoPoint, p2: GeoPoint): int64 {
var dLat: int64 := p2.y - p1.y;
var dLon: int64 := p2.x - p1.x;
// cos der mittleren Breite selbst bestimmen (Tabelle oder Naeherung)
var cosLat: int64 := 602000; // Beispiel: 53 Grad Nord
var dLonCorr: int64 := (dLon * cosLat) / 1000000;
var y: int64 := (dLat * 111319) / 1000000;
var x: int64 := (dLonCorr * 111319) / 1000000;
return IntSqrt(x * x + y * y);
}
==== Umkreis-Rechteck und Verschiebung (#1483) ====
Beide multiplizieren den Längenoffset mit dem Kosinus, statt durch ihn zu teilen — und lassen eine Skalierung aus:
CalculateBoundingBox(Berlin, 10 km):
min = (-74943165, 52430188) max = (101753073, 52609828)
Die Breitenspanne stimmt (0,179640° ≈ 20 km), die Längenspanne reicht über **176 Grad** statt der erwarteten 0,29°. Eine Umkreissuche über dieses Rechteck umfasst den halben Globus.
AddOffsetM(Berlin, 0, 10000) // 10 km nach Norden
-> lon = -88334714706, lat = 52609828
Die Breite stimmt exakt, die Länge ist zerstört — bei einer Bewegung, die gar keine Längenänderung enthält.
==== Bearing kennt nur vier Richtungen (#1483) ====
nach Norden = 0 ok
nach Osten = 90000000 ok
nach Sueden = 180000000 ok
nach Westen = 270000000 ok
Nordost = 0 erwartet 45000000
Die Achsenrichtungen sind fest verdrahtet; dazwischen liefert die Funktion 0.
==== FormatDMS gibt eine Platzhalterzeichenkette (#1483) ====
FormatDMS(52520000, true) = 000° 00' 00" X
Grad, Minuten und Sekunden werden berechnet und anschließend verworfen; zurück kommt ein festes Literal. Die Gegenrichtung ''ParseDMS'' arbeitet korrekt.
----
===== Verwandte Units =====
* [[lyx_-_programmiersprache:units:math|std.math]] — ''IntSqrt'' für eigene Entfernungsrechnungen
* [[lyx_-_programmiersprache:units:rect|std.rect]] — achsenparallele Rechtecke in ebenen Koordinaten
* [[lyx_-_programmiersprache:units:validate:iso6346|std.validate]] — Prüfziffern und Codeformate