Große Koordinaten: wenn die Punktwolke auf halbe Meter rundet.

Das Problem in einer Zahl

Computer speichern Kommazahlen als Gleitkommazahl nach IEEE 754: ein Vorzeichen, ein Exponent und eine Mantisse mit fester Zahl an Stellen. Wie fein eine Zahl aufgelöst wird, hängt deshalb von ihrer Größe ab. Für eine Zahl $x$ mit dem Exponenten $e=\lfloor {\log }_2|x|\rfloor$ ist der Abstand zur nächstgrößeren darstellbaren Zahl

$$\Delta (x)=2^{e-p}$$

mit $p=23$ gespeicherten Mantissenbits bei einfacher Genauigkeit (float32, 4 Byte) und $p=52$ bei doppelter Genauigkeit (float64, 8 Byte). Die Formel gilt für die normalen Zahlen, also für alles, womit ein Vermesser je zu tun hat. Für den Hochwert $N=\mathrm{5.542.928}\text{m}$ gilt $2^{22}\le N<2^{23}$, also $e=22$:

$${\Delta }_{\text{float32}}=2^{22-23}\text{m}=\mathrm{0,5}\text{m}{\Delta }_{\text{float64}}=2^{22-52}\text{m}\approx \mathrm{0,9}\text{nm}$$

In float32 kann dieser Hochwert also nur in Schritten von einem halben Meter gespeichert werden. Jeder Punkt wird auf den nächsten Schritt gerundet, der Fehler beträgt bis zu $\Delta /2=25\text{cm}$. In float64 ist dieselbe Zahl auf unter ein Nanometer genau, dort gibt es kein Problem.

Eine handliche Faustregel: float32 hat 24 geltende Binärstellen. Jede ganze Zahl bis $2^{24}=\mathrm{16.777.216}$ ist darin exakt darstellbar, darüber nicht mehr jede. Das sind gut sieben geltende Dezimalstellen. Ein UTM-Hochwert verbraucht sieben Stellen schon vor dem Komma, für die Nachkommastellen bleibt nichts übrig.

Warum die Schrittweite in Stufen wächst

Zwischen zwei Zweierpotenzen ist die Schrittweite konstant. Alle Werte von 4.194.304 m $(2^{22})$ bis 8.388.608 m $(2^{23})$ liegen in float32 genau 0,5 m auseinander, unterhalb von 4.194.304 m sind es 0,25 m. An der Zweierpotenz selbst ist der Abstand nach unten deshalb nur halb so groß wie der nach oben: Von 4.194.304 m aus ist der nächste kleinere float32-Wert 0,25 m entfernt, der nächste größere 0,5 m. Für die Abschätzung des Rundungsfehlers genügt die Schrittweite nach oben, sie ist die größere der beiden.

Wie grob es wird

Die folgende Tabelle zeigt die Schrittweite in float32 für typische Koordinatengrößen, gerechnet mit der Funktion numpy.spacing, die für positive Zahlen den Abstand zum nächstgrößeren darstellbaren Wert ausgibt. Das ist die Schrittweite im jeweiligen Binärintervall (die ULP, englisch unit in the last place); nur genau an einer Zweierpotenz ist der Abstand nach unten kleiner, siehe oben. Der Rechtswert ist einmal ohne und einmal mit vorangestellter Zonennummer aufgeführt, weil beide Schreibweisen im Umlauf sind. Ein Gauß-Krüger-Rechtswert trägt die Kennziffer des Meridianstreifens als Millionenstelle, im Streifen 3 also zum Beispiel 3.413.000 m.

Drei Dinge fallen auf. Höhen sind unkritisch: Selbst 3.000 m über NHN ergeben in float32 noch eine Schrittweite von 0,24 mm. Gauß-Krüger ist nicht besser als UTM, denn der Hochwert ist in beiden Systemen gleich groß und der Rechtswert mit Kennziffer sogar größer als der UTM-Rechtswert ohne Zone. Alte Gauß-Krüger-Daten mit vorangestellter Kennziffer (3.413.000 m) runden in float32 auf 25 cm, kaum besser als der UTM-Hochwert mit 50 cm. Und die Zonennummer vor dem UTM-Rechtswert kostet noch einmal den Faktor 64: Aus 31 mm werden 2 m.

Wichtig ist: Die Schrittweite gilt für die ganze Zahl, nicht nur für die letzten Stellen. Zwei Punkte, die 11 cm auseinanderliegen, können nach dem Runden denselben Wert haben (Abstand null) oder 50 cm auseinander liegen. Das betrifft also nicht nur die absolute Lage, sondern jede Strecke, jede Fläche und jede Normale, die aus den gerundeten Werten berechnet wird.

Mittlerer und größter Fehler

Die folgenden Formeln setzen voraus, dass die Koordinaten zufällig und ohne Bezug zum float32-Raster liegen. Dann sind die Rundungsfehler näherungsweise gleichmäßig zwischen $-\Delta /2$ und $+\Delta /2$ verteilt, und es gilt

$$\overline{|\epsilon |}=\frac{\Delta }{4}{\epsilon }_{\text{RMS}}=\frac{\Delta }{\sqrt{12}}{\epsilon }_{\text{max}}=\frac{\Delta }{2}$$

Eine Simulation mit einer Million zufälliger Punkte in einem 300 m breiten Streifen bestätigt das: Beim Hochwert um 5.542.928 m beträgt der mittlere Betrag des Fehlers 125 mm, der quadratische Mittelwert (RMS) 144 mm und der größte Fehler 250 mm. Beim Rechtswert um 413.000 m sind es 7,8 mm, 9,0 mm und 15,6 mm. Wer nur den Mittelwert prüft, unterschätzt das Problem also um die Hälfte.

Für die Punkte einer Messung mit vielen Millimetern Streuung über eine große Fläche ist die Annahme gut erfüllt. Sie gilt aber nicht von selbst für jede Punktmenge. Liegen die Werte selbst auf einem Raster, entstehen nur wenige, feste Fehlerwerte. Nachgerechnet für Hochwerte in einem 200 m breiten Streifen ab einem vollen Meter: Ein Raster von 1 m ist in float32 fehlerfrei, weil jedes Vielfache von 0,5 m exakt darstellbar ist. Ein Raster von 10 cm ergibt nur drei verschiedene Fehlerbeträge (0, 10 cm und 20 cm), im Mittel 120 mm, höchstens 200 mm. Erst ein 1-cm-Raster verhält sich wieder wie die Gleichverteilung (125 mm und 250 mm). Und eine kleine Geometrie, die nur ein oder zwei Rasterstufen überdeckt, hat keinen mittleren Fehler im Sinne der Formel: Wie groß ihr Fehler ist, hängt davon ab, wo sie zum Raster liegt. Der nächste Abschnitt zeigt das an einer Wand.

Was aus den gerundeten Punkten wird

Am deutlichsten zeigt sich der Schaden bei allem, was aus einer kleinen Nachbarschaft berechnet wird. Ein Beispiel, nachgerechnet mit simulierten Punkten. Der Aufbau im Einzelnen:

• Wand. Senkrecht, die wahre Normale liegt waagerecht und um 45° zwischen Rechts- und Hochwert gedreht. Die Mitte des Ausschnitts liegt bei Rechtswert 413.000,3 m, Hochwert 5.542.928,3 m und Höhe 150 m.
• Punkte. 2.000 Punkte, in der Wandfläche gleichmäßig zufällig über einen Ausschnitt von 20 cm mal 20 cm verteilt (waagerecht entlang der Wand und senkrecht).
• Messrauschen. Jeder Punkt ist senkrecht zur Wand um einen normalverteilten Betrag mit 1 mm Standardabweichung verschoben.
• Zufallszahlen. Generator numpy.random.default_rng mit dem festen Startwert (Seed) 1, die Rechnung ist damit wiederholbar.
• Normale. Vom Ausschnitt wird zuerst der Schwerpunkt abgezogen, aus den zentrierten Koordinaten wird die 3 × 3-Streumatrix gebildet. Die Normale ist der Eigenvektor zum kleinsten Eigenwert, gerechnet in float64. Die Rundung auf float32 betrifft also nur die Eingangskoordinaten.
• Kennzahlen. Fehler der Normale ist der Winkel zwischen berechneter und wahrer Normale, unabhängig vom Vorzeichen. Streuung um die Ebene ist die Wurzel aus dem kleinsten Eigenwert geteilt durch die Punktzahl, also der quadratische Mittelwert der Abstände zur angepassten Ebene.

In absoluten float32-Werten bleiben von 2.000 Hochwerten genau zwei übrig. Die Wand zerfällt in zwei Scheiben im Abstand von 0,5 m, die Normale kippt um mehr als ein Drittel eines rechten Winkels, und die Streuung um die angepasste Ebene ist dreißigmal so groß wie das Messrauschen. Jede Auswertung, die auf Normalen aufbaut (Ebenen suchen, Punkt-zu-Ebene-Registrierung, Bauteilerkennung, Schattierung in der Ansicht), arbeitet dann mit falschen Werten. Relativ zu einem Ursprung gerechnet sind die Ergebnisse von float64 nicht zu unterscheiden.

Die Zahlen sind das Ergebnis dieser einen Simulation, keine Grenzwerte für jede Wand. Mit 100 verschiedenen Startwerten bei derselben Lage schwankt der Fehler der Normale in absoluten float32-Werten zwischen 36,2° und 36,9° und die Streuung zwischen 32,3 mm und 34,6 mm; in float64 und relativ gerechnet bleibt der Winkel unter 0,07° und die Streuung zwischen 0,97 mm und 1,04 mm. Viel stärker als vom Zufall hängt das Ergebnis von der Lage der Wand zum 0,5-m-Raster ab. Der Ausschnitt überdeckt im Hochwert nur rund 14 cm. Liegt seine Mitte bei 5.542.928,25 m, verteilen sich die Punkte auf zwei Rasterwerte (36,8° Fehler, 21,6 mm Streuung). Liegt sie bei 5.542.928,05 m, 5.542.928,15 m oder 5.542.928,40 m, landen alle Punkte auf demselben Hochwert: Die berechnete Normale zeigt dann genau in Hochwertrichtung und ist um 45° falsch, die Streuung um die Ebene ist aber null. Die Wand sieht in der Auswertung vollkommen eben aus, obwohl sie falsch steht. Eine kleine Streuung ist in absoluten float32-Werten also kein Beleg für eine richtige Ebene.

Wo float32 im Spiel ist

Am Scanner selbst entsteht das Problem nicht. Ein Laserscanner misst im eigenen Gerätesystem, die Koordinaten sind dort so groß wie die Entfernung vom Gerät, also einige Meter bis einige hundert Meter. Selbst wenn ein Scan in diesem System in float32 gespeichert ist, liegt die Schrittweite bei 300 m Entfernung bei 0,03 mm. Kritisch wird es erst, wenn die Punkte ins Landessystem gedreht und verschoben werden und das Ergebnis an einer der folgenden Stellen in float32 landet:

• Grafikkarte. Die Grafikschnittstellen im Browser kennen nur einfache Genauigkeit (siehe nächster Abschnitt). Eine Punktwolke, die mit absoluten UTM-Werten an die Grafikkarte geht, wird dort gerundet, auch wenn sie im Hauptspeicher richtig vorliegt.
• Dateiformate. E57 erlaubt für Koordinaten Gleitkommazahlen in einfacher oder doppelter Genauigkeit oder skalierte Ganzzahlen (Festkomma). PLY schreibt den Datentyp jeder Eigenschaft in den Dateikopf, erlaubt sind unter anderem float mit 4 Byte und double mit 8 Byte. Wer eine Punktwolke im Landessystem als E57 mit einfacher Genauigkeit oder als PLY mit float schreibt, rundet sie beim Speichern. Ob das passiert ist, verrät bei PLY ein Blick in den Textkopf der Datei.
• Zwischenspeicher. Übersichtswolken, Kacheln und Octrees halten Millionen Punkte bereit. Wer dort aus Platzgründen float32 wählt, bekommt den Fehler in jede schnelle Vorschau und jede Messung, die darauf aufsetzt.

LAS arbeitet nach einem anderen Prinzip, das einen eigenen Abschnitt verdient.

LAS: Festkomma statt Gleitkomma

LAS speichert jede Koordinate als vorzeichenbehaftete 32-Bit-Ganzzahl. Im Dateikopf stehen je Achse ein Maßstab (scale) und ein Versatz (offset), beide in doppelter Genauigkeit. Die Spezifikation schreibt die Rückrechnung so vor:

$$X_{\text{Koordinate}}=X_{\text{Datensatz}}·X_{\text{Maßstab}}+X_{\text{Versatz}}$$

Damit ist die Auflösung überall gleich, nämlich der Maßstab. Bei einem Maßstab von 0,001 liegen alle darstellbaren Werte genau 1 mm auseinander, egal ob die Koordinate 100 m oder 5.542.928 m groß ist. Der Rundungsfehler beträgt höchstens 0,5 mm, im quadratischen Mittel 0,29 mm. Das ist ein grundsätzlich anderes Verhalten als bei float32, wo die Schrittweite mit der Zahl wächst.

Der Versatz hat in LAS eine andere Aufgabe als der Rechenursprung bei Gleitkommazahlen: Er sorgt dafür, dass die Ganzzahlen in den Wertebereich von int32 passen. Dieser reicht bis $2^{31}-1=\mathrm{2.147.483.647}$. Bei einem Maßstab von 1 mm deckt eine LAS-Datei damit 2.147.483,647 m beiderseits des Versatzes ab. Ein Hochwert von 5.542.928,184 m würde ohne Versatz die Ganzzahl 5.542.928.184 ergeben und passt nicht hinein. Mit einem Versatz von 5.542.928 m wird daraus die Ganzzahl 184. Die Spezifikation nennt genau diesen Fall: Im Allgemeinen seien die Versatzwerte null, für bestimmte Fälle reiche die Auflösung der Ganzzahlen aber nicht aus.

Für die Praxis heißt das: Bei LAS ist der Maßstab die entscheidende Angabe. Ein Maßstab von 0,01 rundet jede Koordinate um bis zu 5 mm, ein Maßstab von 0,001 um bis zu 0,5 mm; welcher genügt, entscheidet die Genauigkeitsanforderung des Auftrags. Das Format selbst prüft das nicht: Auch eine LAS-Datei mit Maßstab 0,1 ist gültig und rundet dann auf 5 cm. Wer eine LAS-Datei weitergibt, sollte Maßstab und Versatz aus dem Kopf ablesen und im Bericht nennen.

Textformate wie PTS oder XYZ sind von beidem nicht betroffen, dort bestimmt die Zahl der geschriebenen Nachkommastellen die Auflösung. Vier Nachkommastellen ergeben einen größten Rundungsfehler von 0,05 mm. Erst beim Einlesen entscheidet das lesende Programm, in welcher Zahlenart es die Werte ablegt; legt es sie in float32 ab, ist die Rundung aus der Tabelle oben wieder da, auch wenn die Datei auf 0,1 mm genau war.

Warum nicht einfach überall float64?

float64 ist genauso Teil von IEEE 754 wie float32 und auf dem Prozessor die übliche Zahlenart. Programmiersprachen und Rechenbibliotheken arbeiten standardmäßig damit. Dass trotzdem float32 verwendet wird, hat drei Gründe:

• Die Grafikschnittstellen im Browser kennen es nicht. Die Spezifikationen von WebGL 1.0 und 2.0 sehen für Eckpunktdaten und Texturen nur Gleitkommaformate bis 32 Bit vor, einen Typ für doppelte Genauigkeit gibt es nicht. Auch die Shadersprache von WebGPU (WGSL) kennt als Gleitkommatypen nur f32 (IEEE 754 binary32) und, als Erweiterung, f16. Eine 3D-Ansicht im Browser übergibt ihre Punkte deshalb als float32.
• Grafikkarten rechnen float64 langsam. Wie stark, steht für NVIDIA-Karten im CUDA C++ Programming Guide, Fassung 12.4, Abschnitt 5.4.1, in einer Tabelle des Durchsatzes nativer Rechenbefehle (Zahl der Ergebnisse je Takt und Multiprozessor). Die folgende Tabelle zeigt daraus das Verhältnis der Zeilen für Addition, Multiplikation und Multiplizieren-Addieren in 64 Bit zu denselben Befehlen in 32 Bit.
• Der Speicher verdoppelt sich. 100 Millionen Punkte belegen mit x, y und z in float32 1,2 GB, in float64 2,4 GB. Das gilt für Arbeitsspeicher, Grafikspeicher, Zwischenspeicher auf der Platte und jede Übertragung. Bei großen Projekten entscheidet das, ob die Ansicht flüssig läuft oder der Speicher ausgeht.

Die Rechenfähigkeit ist NVIDIAs Kennung der Chipgeneration, welche Generation eine bestimmte Karte hat, steht in deren Datenblatt. Die Spanne reicht also von halbem Durchsatz bis zu einem Vierundsechzigstel. Das sind theoretische Werte für einzelne Befehle, keine Laufzeiten ganzer Programme: Wie schnell eine Anwendung tatsächlich ist, hängt auch von Speicherbandbreite, Datenzugriff und Aufbau des Rechenkerns ab, und neuere Chipgenerationen als 9.0 stehen in dieser Fassung der Tabelle nicht. Für rechenintensive GPU-Kernel kann FP64 auf diesen Architekturen deutlich langsamer sein; für die reine Koordinatenverwaltung ist eine Relativierung meist die effizientere Lösung.

Die übliche Lösung ist deshalb eine Arbeitsteilung: float64 dort, wo gerechnet und dauerhaft gespeichert wird, float32 dort, wo es schnell und klein sein muss, dann aber nur für Werte relativ zu einem Ursprung. Wie das funktioniert, zeigt der Abschnitt zum Rechenursprung.

Auch Rechnungen verlieren Stellen

Nicht nur das Speichern rundet. Jede Rechenoperation rundet ihr Ergebnis auf die Zahlenart, in der gerechnet wird, und bei großen Zahlen summieren sich diese Rundungen. Zwei Fälle kommen bei Punktwolken ständig vor: der Schwerpunkt einer Punktmenge und ihre Streuung.

Summen und Schwerpunkte

Der Schwerpunkt ist die Summe aller Koordinaten geteilt durch ihre Anzahl. Bei einer Million Hochwerten um 5,5 Millionen Meter wächst die Zwischensumme gegen Ende auf rund $\mathrm{5,5}·{10}^{12}$ m. In float32 beträgt die Schrittweite bei dieser Größe 524.288 m, jeder weitere Summand wird beim Addieren grob gerundet. Nachgerechnet mit einer Million Hochwerten in einem 100 m breiten Streifen:

Die einfache Schleife liegt um 25 km daneben. Die Rechenbibliothek numpy summiert deshalb nicht der Reihe nach, sondern paarweise in Teilsummen, ihre Dokumentation weist ausdrücklich darauf hin und empfiehlt für float32-Daten, die Summe in float64 zu bilden. Goldberg leitet in seinem Standardaufsatz die Fehlerschranken her: Beim Aufaddieren der Reihe nach kann der Fehler eines Summanden mit der Zahl der folgenden Additionen wachsen, beim Summationsverfahren nach Kahan bleibt er je Summand unabhängig von der Anzahl auf das Doppelte der Maschinengenauigkeit beschränkt, bis auf einen kleinen Zusatzterm zweiter Ordnung. Am einfachsten ist aber auch hier der Ursprung: Relativ gerechnet sind die Summanden klein, und schon die einfache Schleife ist genau.

Streuungen

Selbst float64 schützt nicht vor jeder Falle. Die Streuung einer Punktmenge lässt sich auf zwei Wegen berechnen:

$$s^2=\frac{1}{n}\sum _{i=1}^n{(x_i-\overline{x})}^2\text{oder}{s}^2=\frac{1}{n}\sum _{i=1}^n{x}_i^2-{\overline{x}}^2$$

Mathematisch sind beide gleich, die zweite Form kommt mit einem einzigen Durchgang durch die Daten aus und ist deshalb beliebt. Bei Landeskoordinaten ist sie unbrauchbar: Das Quadrat eines Hochwerts liegt bei $3·{10}^{13}{\text{m}}^2$, die Schrittweite von float64 beträgt dort 0,0039 m². Die gesuchte Varianz einer Wand mit 1 mm Rauschen ist aber nur ${10}^{-6}{\text{m}}^2$ groß und geht in der Differenz zweier fast gleicher Riesenzahlen unter. Nachgerechnet mit 100.000 Werten ergibt die zweite Form eine Standardabweichung von 88 mm statt 1,0 mm. Die erste Form und die zweite Form mit relativen Koordinaten liefern beide 1,0 mm. Dieselbe Falle steckt in jeder Kovarianzmatrix, aus der Normalen, Ebenen oder Hauptachsen berechnet werden, etwa in der Form $\overline{xy}-\overline{x}\overline{y}$.

Daraus folgt eine Regel, die über die Zahlenart hinausgeht: Auch eine Rechnung, die durchgehend in float64 läuft, ist nicht von selbst stabil. Die Differenz zweier großer, fast gleicher Summen löscht in jeder Zahlenart die führenden Stellen aus, float64 hat nur mehr Stellen zu verlieren. Für Varianzen, Kovarianzen und Hauptachsen nimmt man deshalb Verfahren, die vor dem Quadrieren zentrieren: entweder zwei Durchgänge (erst Schwerpunkt, dann Abweichungen vom Schwerpunkt) oder, wenn die Daten nur einmal durchlaufen werden können, das Verfahren nach Welford. Es führt Mittelwert und Quadratsumme der Abweichungen laufend mit: Für jeden neuen Wert $x_k$ wird $d=x_k-m_{k-1}$, $m_k=m_{k-1}+d/k$ und $M_k=M_{k-1}+d·(x_k-m_k)$ gerechnet, am Ende ist $s^2=M_n/n$. Nachgerechnet mit einer weiteren Stichprobe von 100.000 Hochwerten mit 1 mm Rauschen, alles in float64: Die Ein-Durchgang-Formel liefert 108 mm, das Verfahren nach Welford und die Rechnung in zwei Durchgängen liefern beide 1,0 mm. Am einfachsten bleibt der Rechenursprung, nach dem Abziehen sind die Zahlen klein und die Auslöschung harmlos.

Woran man den Fehler erkennt

• Die Punktwolke zerfällt in Streifen oder Treppen mit festem Abstand, besonders deutlich in der Richtung mit der größeren Koordinate (meist der Hochwert). Der Abstand der Streifen ist die Schrittweite aus der Tabelle oben, beim Hochwert also 0,5 m.
• Gemessene Strecken ändern sich in Sprüngen von einigen Zentimetern bis Dezimetern, je nachdem, welche Punkte man anklickt.
• Ebenen, die an solche Punkte angepasst werden, haben eine unplausibel große Streuung, Normalen kippen.
• Beim Drehen oder Heranzoomen zittern die Punkte oder springen in festen Abständen.

Das Zittern hat eine eigene Ursache. Es entsteht nicht durch die gespeicherten Punkte allein, sondern in der Kamerarechnung: Die Grafikkarte zieht von jedem Punkt die Lage des Augpunkts ab und dreht das Ergebnis in die Blickrichtung. Steht der Augpunkt in absoluten Koordinaten in float32, kann er sich nur in Schritten von 0,5 m bewegen. Nachgerechnet: Fünf Augpunkte mit Hochwerten von 5.542.931,81 m bis 5.542.932,21 m im Abstand von 10 cm landen in float32 alle auf 5.542.932,00 m. Die Kamera bleibt stehen und springt dann um einen halben Meter, die ganze Szene ruckt mit. Deshalb genügt es nicht, nur die Punkte relativ zu speichern: Augpunkt, Drehpunkt und jede Verschiebung, die an die Grafikkarte geht, müssen ebenfalls relativ zum selben Ursprung übergeben werden.

Ein schneller Selbsttest: Dieselbe Szene einmal in lokalen Koordinaten und einmal in Landeskoordinaten öffnen. Sieht nur die zweite Fassung zerrissen aus, ist es die Zahlendarstellung und nicht der Scan. Ein zweiter Test ist das Abfragen zweier benachbarter Punkte: Haben viele Punkte auf einer Fläche exakt denselben Hochwert, oder unterscheiden sich die Hochwerte nur in Vielfachen von 0,5 m, ist die Wolke irgendwo durch float32 gegangen.

Die Lösung: ein Rechenursprung

Für Abstände kommt es nur auf Differenzen an. Zieht man von allen Punkten denselben Ursprung ${\mathrm{𝐱}}_0$ ab, ändert sich keine Strecke:

$${\mathrm{𝐱}}_i^{\prime }={\mathrm{𝐱}}_i-{\mathrm{𝐱}}_0\Vert {\mathrm{𝐱}}_i^{\prime }-{\mathrm{𝐱}}_j^{\prime }\Vert =\Vert {\mathrm{𝐱}}_i-{\mathrm{𝐱}}_j\Vert$$

Die relativen Werte sind aber klein, und damit ist auch ihre float32-Schrittweite klein. Wie klein, hängt nur noch von der Ausdehnung des Projekts ab. Maßgeblich ist dabei jede Achse für sich: float32 rundet jede Koordinate einzeln, die Schrittweite folgt aus dem Betrag der jeweiligen Komponente, also aus dem relativen Rechts-, Hoch- oder Höhenwert, nicht aus dem räumlichen Abstand zum Ursprung. Die folgende Tabelle ist mit einer Million zufälliger Hochwerte um 5.542.928 m gerechnet, der Ursprung ist jeweils die auf 100 m abgerundete Mitte. Die mittlere Spalte ist der größte vorkommende Betrag der relativen Koordinate und damit eine obere Schranke für jede Komponente in diesem Projekt:

Selbst bei einem Projekt von 5 km Ausdehnung bleibt der Fehler bei gut einem Zehntel Millimeter, also eine Größenordnung unter einem Millimeter, der Einheit, in der das Messrauschen terrestrischer Laserscanner üblicherweise angegeben wird. Bei Bauwerksaufnahmen von einigen hundert Metern liegt er bei Tausendsteln eines Millimeters. Für einen ganzen Landkreis in einem einzigen Ursprung reicht float32 dagegen nicht mehr ohne Weiteres: Hat eine Koordinatenkomponente den Betrag 50 km, ist ihre Schrittweite wieder 3,9 mm. Ein Punkt, der räumlich 50 km entfernt liegt, aber schräg zu den Achsen, hat kleinere Komponenten (bei 45° je rund 35 km) und damit je Achse höchstens dieselbe, oft eine kleinere Schrittweite.

Die Reihenfolge entscheidet

Entscheidend ist: erst in doppelter Genauigkeit den Ursprung abziehen, dann auf float32 verkleinern. Wer erst verkleinert und dann abzieht, hat die Information schon verloren. Im Rechenbeispiel unten werden die Hochwerte 5.542.928,18 m und 5.542.928,29 m in dieser Reihenfolge zu 28,180000 m und 28,290000 m. In der falschen Reihenfolge werden sie zuerst zu 5.542.928,0 m und 5.542.928,5 m gerundet, nach dem Abziehen des Ursprungs stehen dort 28,0 m und 28,5 m: kleine Zahlen, aber mit dem Fehler der großen.

Wahl und Ablage des Ursprungs

Der Ursprung selbst wird einmal in doppelter Genauigkeit neben den Daten abgelegt. Sinnvoll ist ein runder Wert in der Nähe der Projektmitte, etwa die auf 100 m abgerundete Mitte. Dann lässt sich jede relative Koordinate im Kopf zurückrechnen, und Exporte, Koordinatenanzeigen und Berichte können wieder volle Landeskoordinaten ausgeben. Die Mitte ist besser als eine Ecke, weil der größte Abstand zum Ursprung dann nur die halbe Ausdehnung beträgt.

Ein Rechenursprung ist eine reine Verschiebung der Zahlen. Er ändert kein Bezugssystem, keinen Maßstab und keine Höhenreduktion: Strecken aus UTM-Koordinaten bleiben auch nach dem Abziehen des Ursprungs Strecken in der Abbildung und müssen bei Bedarf weiterhin wie im Beitrag zum UTM-Maßstab reduziert werden.

Ein Risiko gibt es: zwei Datensätze mit verschiedenen Ursprüngen. Werden zwei relativ gespeicherte Wolken zusammengeführt, ohne dass ein Programm die Ursprünge berücksichtigt, liegen sie genau um die Differenz der Ursprünge auseinander. Liegen die Ursprünge auf vollen 100 m, ist der Versatz ein glattes Vielfaches davon und fällt meist sofort auf. Wer Daten mit Versatz weitergibt, schreibt den Versatz deshalb in eine Begleitdatei und in den Bericht.

Das Prinzip ist nicht neu. Bietet ein Programm beim Import einen „globalen Versatz“ oder „Shift“ an, steckt derselbe Gedanke dahinter, nur muss man dort selbst daran denken, ihn beim Export wieder zu addieren. Der Versatz im LAS-Dateikopf ist dagegen, wie oben beschrieben, Teil einer Festkommadarstellung und löst eine andere Aufgabe.

Rechenbeispiel

Zwei Punkte auf einer Wand, 11 cm auseinander, Hochwerte 5.542.928,18 m und 5.542.928,29 m. Ursprung $N_0=\mathrm{5.542.900}\text{m}$.

Bei 28 m beträgt die float32-Schrittweite 0,0019 mm, die relative Fassung ist also auf etwa 0,002 mm genau. Genauer gerechnet speichert float32 die beiden Werte als 28,180000305 m und 28,290000916 m, der Abstand ergibt sich zu 0,110000610 m, also 0,0006 mm zu lang. Die absolute Fassung macht aus 11 cm einen halben Meter.

Für die Praxis

• Ausgabeformat prüfen. Bei E57 und PLY nachsehen, ob Koordinaten im Landessystem in einfacher Genauigkeit geschrieben werden. Bei LAS Maßstab und Versatz aus dem Dateikopf ablesen. Bei Textformaten die Zahl der Nachkommastellen prüfen.
• Versatz dokumentieren. Wird für ein CAD-Programm ein Versatz abgezogen, gehört er in eine Begleitdatei und in den Bericht, mit Achse und Vorzeichen.
• Kontrolle mit bekannten Punkten. Nach jedem Austausch einige Zielmarken oder Passpunkte abfragen und mit den Sollkoordinaten vergleichen. Abweichungen in Vielfachen von 0,25 m oder 0,5 m deuten auf float32, glatte Abweichungen von 100 m oder 1.000 m auf einen vergessenen Versatz.
• Selbst rechnen mit Ursprung. Wer eigene Auswertungen schreibt, zieht zuerst einen Ursprung ab, rechnet Schwerpunkte und Kovarianzen relativ und addiert den Ursprung erst bei der Ausgabe wieder.