std.geo — Geokoordinaten

→ 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

<WRAP info> Seit lyxc 1.0.20F stimmen die Entfernungen (#1481). Dieselbe Strecke Berlin–Hamburg wie in der früheren Fehlermessung:

Die Luftlinie beträgt rund 255 km — beide Verfahren liefern damit denselben, plausiblen Wert. Ursache war eine um Faktor 10 verrechnete Kosinustabelle. </WRAP>

import std.geo;

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.

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 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.

ParseDMS rundet ganzzahlig; die Abweichung liegt bei etwa 13 Mikrograd, also gut einem Zentimeter.

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.

Für neuen Code die GeoBBox-Fassungen verwenden — die älteren Helfer mit losen GeoPoint-Paaren machen dasselbe, nur umständlicher.

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.

Geprüft mit lyxc 1.0.21A, Referenzstrecken aus bekannten Städtekoordinaten.

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:

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.

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 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);
}

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.

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(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.

• std.math — IntSqrt für eigene Entfernungsrechnungen
• std.rect — achsenparallele Rechtecke in ebenen Koordinaten
• std.validate — Prüfziffern und Codeformate