lsa-planer

LSA-Planer Professional – Planungssoftware für Lichtsignalanlagen nach RiLSA 2015 und § 45 StVO. EUPL-1.2.

/ src domain geometrie projektion.ts

11,9 KB Rohdatei
src/domain/geometrie/projektion.ts — 334 Zeilen
1 import type { Meters } from '../units';
2
3 /**
4 * Umrechnung geografischer Koordinaten in ein metrisches Bezugssystem.
5 *
6 * WARUM NICHT WEB-MERCATOR
7 *
8 * Die naheliegende Wahl waere EPSG:3857 (Web-Mercator), das jeder Kartendienst
9 * spricht. Dort sind Karteneinheiten aber KEINE Meter: Bei der geografischen
10 * Breite phi entspricht eine Karteneinheit nur cos(phi) Metern am Boden. In
11 * Deutschland - etwa bei 52 Grad Nord - sind das 0,62 m. Ein daraus
12 * abgegriffener Raeumweg waere um mehr als 60 Prozent zu gross.
13 *
14 * Fuer die Vermassung wird deshalb ETRS89 / UTM Zone 32N (EPSG:25832) genutzt.
15 * Das ist zugleich das Bezugssystem der deutschen Landesvermessung; die
16 * Karteneinheiten sind Meter, der Massstabsfehler liegt bei hoechstens etwa
17 * 0,4 Promille - also unter einem Millimeter je Meter und damit weit unterhalb
18 * der Ablesegenauigkeit eines Lageplans.
19 *
20 * Zone 32N deckt 6 bis 12 Grad oestlicher Laenge ab und damit den groessten Teil
21 * Deutschlands. Fuer Gebiete weiter westlich oder oestlich stehen die Zonen 31
22 * und 33 zur Verfuegung.
23 */
24
25 /** Geografische Koordinate in Grad. */
26 export interface GeoPunkt {
27 readonly breite: number;
28 readonly laenge: number;
29 }
30
31 /** Koordinate in einem metrischen Bezugssystem. */
32 export interface UtmPunkt {
33 readonly ost: Meters;
34 readonly nord: Meters;
35 readonly zone: number;
36 }
37
38 // WGS84 bzw. GRS80 - die Unterschiede liegen im Zentimeterbereich und sind
39 // fuer die Vermassung eines Knotenpunkts ohne Belang.
40 const A = 6378137.0;
41 const F = 1 / 298.257223563;
42 const K0 = 0.9996;
43 const FALSE_EAST = 500000;
44
45 const E2 = F * (2 - F);
46 const EP2 = E2 / (1 - E2);
47
48 /** UTM-Zone zu einer geografischen Laenge. */
49 export function utmZone(laenge: number): number {
50 return Math.floor((laenge + 180) / 6) + 1;
51 }
52
53 /** EPSG-Code der ETRS89/UTM-Zone. */
54 export function epsgFuerZone(zone: number): number {
55 return 25800 + zone;
56 }
57
58 /**
59 * Rechnet eine geografische Koordinate in UTM um.
60 * Reihenentwicklung nach Krueger; fuer Abstaende bis einige Kilometer vom
61 * Mittelmeridian auf Millimeter genau.
62 */
63 export function nachUtm(punkt: GeoPunkt, zone = utmZone(punkt.laenge)): UtmPunkt {
64 const phi = (punkt.breite * Math.PI) / 180;
65 const lambda = (punkt.laenge * Math.PI) / 180;
66 const lambda0 = (((zone - 1) * 6 - 180 + 3) * Math.PI) / 180;
67
68 const sinPhi = Math.sin(phi);
69 const cosPhi = Math.cos(phi);
70 const tanPhi = Math.tan(phi);
71
72 const N = A / Math.sqrt(1 - E2 * sinPhi * sinPhi);
73 const T = tanPhi * tanPhi;
74 const C = EP2 * cosPhi * cosPhi;
75 const Aa = (lambda - lambda0) * cosPhi;
76
77 const M =
78 A *
79 ((1 - E2 / 4 - (3 * E2 * E2) / 64 - (5 * E2 * E2 * E2) / 256) * phi -
80 ((3 * E2) / 8 + (3 * E2 * E2) / 32 + (45 * E2 * E2 * E2) / 1024) * Math.sin(2 * phi) +
81 ((15 * E2 * E2) / 256 + (45 * E2 * E2 * E2) / 1024) * Math.sin(4 * phi) -
82 ((35 * E2 * E2 * E2) / 3072) * Math.sin(6 * phi));
83
84 const ost =
85 FALSE_EAST +
86 K0 *
87 N *
88 (Aa +
89 ((1 - T + C) * Aa ** 3) / 6 +
90 ((5 - 18 * T + T * T + 72 * C - 58 * EP2) * Aa ** 5) / 120);
91
92 const nord =
93 K0 *
94 (M +
95 N *
96 tanPhi *
97 ((Aa * Aa) / 2 +
98 ((5 - T + 9 * C + 4 * C * C) * Aa ** 4) / 24 +
99 ((61 - 58 * T + T * T + 600 * C - 330 * EP2) * Aa ** 6) / 720));
100
101 return { ost, nord, zone };
102 }
103
104 /** Umkehrung - fuer die Anzeige und zur Pruefung der Hinrechnung. */
105 export function nachGeo(punkt: UtmPunkt): GeoPunkt {
106 const x = punkt.ost - FALSE_EAST;
107 const y = punkt.nord;
108 const lambda0 = (((punkt.zone - 1) * 6 - 180 + 3) * Math.PI) / 180;
109
110 const e1 = (1 - Math.sqrt(1 - E2)) / (1 + Math.sqrt(1 - E2));
111 const M = y / K0;
112 const mu = M / (A * (1 - E2 / 4 - (3 * E2 * E2) / 64 - (5 * E2 * E2 * E2) / 256));
113
114 const phi1 =
115 mu +
116 ((3 * e1) / 2 - (27 * e1 ** 3) / 32) * Math.sin(2 * mu) +
117 ((21 * e1 * e1) / 16 - (55 * e1 ** 4) / 32) * Math.sin(4 * mu) +
118 ((151 * e1 ** 3) / 96) * Math.sin(6 * mu) +
119 ((1097 * e1 ** 4) / 512) * Math.sin(8 * mu);
120
121 const sinPhi1 = Math.sin(phi1);
122 const cosPhi1 = Math.cos(phi1);
123 const tanPhi1 = Math.tan(phi1);
124
125 const N1 = A / Math.sqrt(1 - E2 * sinPhi1 * sinPhi1);
126 const T1 = tanPhi1 * tanPhi1;
127 const C1 = EP2 * cosPhi1 * cosPhi1;
128 const R1 = (A * (1 - E2)) / (1 - E2 * sinPhi1 * sinPhi1) ** 1.5;
129 const D = x / (N1 * K0);
130
131 const phi =
132 phi1 -
133 ((N1 * tanPhi1) / R1) *
134 ((D * D) / 2 -
135 ((5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * EP2) * D ** 4) / 24 +
136 ((61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * EP2 - 3 * C1 * C1) * D ** 6) / 720);
137
138 const lambda =
139 lambda0 +
140 (D -
141 ((1 + 2 * T1 + C1) * D ** 3) / 6 +
142 ((5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * EP2 + 24 * T1 * T1) * D ** 5) / 120) /
143 cosPhi1;
144
145 return { breite: (phi * 180) / Math.PI, laenge: (lambda * 180) / Math.PI };
146 }
147
148 /** Ausschnitt in Kartenkoordinaten - Reihenfolge wie im WMS-Parameter BBOX. */
149 export interface Ausschnitt {
150 readonly minOst: Meters;
151 readonly minNord: Meters;
152 readonly maxOst: Meters;
153 readonly maxNord: Meters;
154 readonly zone: number;
155 }
156
157 /**
158 * Quadratischer Ausschnitt um einen Punkt mit vorgegebener Kantenlaenge.
159 * Weil in UTM gerechnet wird, ist die Kantenlaenge tatsaechlich in Metern am
160 * Boden zu verstehen.
161 */
162 export function ausschnittUm(mitte: GeoPunkt, kantenlaenge: Meters): Ausschnitt {
163 const utm = nachUtm(mitte);
164 const halb = kantenlaenge / 2;
165 return {
166 minOst: utm.ost - halb,
167 minNord: utm.nord - halb,
168 maxOst: utm.ost + halb,
169 maxNord: utm.nord + halb,
170 zone: utm.zone,
171 };
172 }
173
174 /**
175 * Teilausschnitt eines in Kacheln zerlegten Ausschnitts.
176 *
177 * WOZU
178 *
179 * Ein Kartendienst rechnet ein grosses Bild in einem Stueck aus, und der
180 * Aufwand waechst mit der Flaeche; jenseits einiger Millionen Bildpunkte
181 * bricht der Abruf ab oder laeuft in die Zeitgrenze. Ein grosser Ausschnitt
182 * wird deshalb kachelweise geholt und danach zusammengesetzt. Diese Funktion
183 * liefert den Kartenausschnitt EINER Kachel.
184 *
185 * Zerlegt wird in `spalten` mal `zeilen` gleich grosse Kacheln. Gezaehlt wird
186 * wie im Bild: `spalte` 0 liegt links, `zeile` 0 liegt OBEN.
187 *
188 * DIE STELLE, AN DER ES SCHIEFGEHT
189 *
190 * Bildzeile 0 liegt oben, also am NOERDLICHEN Rand des Ausschnitts. Die
191 * Nordwerte laufen der Bildzeile damit ENTGEGEN: Je groesser die Zeile, desto
192 * kleiner der Nordwert. Wird das vertauscht, ist das zusammengesetzte Mosaik
193 * senkrecht gespiegelt - Zeile fuer Zeile in sich richtig, in der Reihenfolge
194 * aber verkehrt.
195 *
196 * Auffallen wuerde das niemandem: Ein Luftbild hat keine Leserichtung. Der
197 * Anwender zeichnet seine Fahrlinien in ein Bild, in dem die noerdliche Zufahrt
198 * unten liegt, misst darin voellig plausible Laengen - und vermasst dabei die
199 * falschen Wege. Deshalb wird hier von `maxNord` nach unten gerechnet und nicht
200 * von `minNord` nach oben.
201 *
202 * Unbrauchbare Angaben werden begrenzt statt zurueckgewiesen: Der Rueckgabewert
203 * ist immer ein echter Teil des Gesamtausschnitts, nie einer daneben. Eine
204 * Kachel ausserhalb des Bildes waere ein Bildteil aus einer fremden Gegend -
205 * und die faellt im Luftbild ebenso wenig auf wie die Spiegelung.
206 */
207 export function unterAusschnitt(
208 gesamt: Ausschnitt,
209 spalte: number,
210 zeile: number,
211 spalten: number,
212 zeilen: number,
213 ): Ausschnitt {
214 const anzahlSpalten = anzahl(spalten);
215 const anzahlZeilen = anzahl(zeilen);
216 const s = index(spalte, anzahlSpalten);
217 const z = index(zeile, anzahlZeilen);
218
219 return {
220 // Ost laeuft mit der Spalte: links ist Westen.
221 minOst: teile(gesamt.minOst, gesamt.maxOst, s, anzahlSpalten),
222 maxOst: teile(gesamt.minOst, gesamt.maxOst, s + 1, anzahlSpalten),
223 // Nord laeuft der Zeile entgegen: Zeile 0 ist der NOERDLICHE Rand.
224 maxNord: teile(gesamt.maxNord, gesamt.minNord, z, anzahlZeilen),
225 minNord: teile(gesamt.maxNord, gesamt.minNord, z + 1, anzahlZeilen),
226 zone: gesamt.zone,
227 };
228 }
229
230 /**
231 * Kachelgrenze als Anteil zwischen zwei Raendern.
232 *
233 * Bewusst als Anteil des Ganzen gerechnet und nicht als "Anfang plus i mal
234 * Kachelbreite": So trifft die letzte Kante den Rand des Gesamtausschnitts
235 * genau. Aufsummierte Kachelbreiten liessen dort einen Spalt von der Groesse
236 * des aufgelaufenen Rundungsfehlers - im Mosaik ein Versatz zwischen den
237 * Kacheln, der sich als scheinbarer Zeichenfehler auf jeden Weg legt.
238 */
239 function teile(von: Meters, bis: Meters, i: number, kacheln: number): Meters {
240 return von + (bis - von) * (i / kacheln);
241 }
242
243 /** Kachelzahl je Richtung; weniger als eine Kachel gibt es nicht. */
244 function anzahl(wert: number): number {
245 if (!Number.isFinite(wert)) return 1;
246 return Math.max(1, Math.floor(wert));
247 }
248
249 /** Kachelnummer, auf das vorhandene Raster begrenzt. */
250 function index(wert: number, kacheln: number): number {
251 if (!Number.isFinite(wert)) return 0;
252 return Math.min(kacheln - 1, Math.max(0, Math.floor(wert)));
253 }
254
255 /** BBOX-Parameter fuer einen WMS-Aufruf. */
256 export function bboxParameter(ausschnitt: Ausschnitt): string {
257 return [ausschnitt.minOst, ausschnitt.minNord, ausschnitt.maxOst, ausschnitt.maxNord]
258 .map((wert) => wert.toFixed(3))
259 .join(',');
260 }
261
262 /** Kantenlaenge des Ausschnitts in Metern. */
263 export function kantenlaenge(ausschnitt: Ausschnitt): Meters {
264 return ausschnitt.maxOst - ausschnitt.minOst;
265 }
266
267 /** Eine Zahl mit oder ohne Nachkommastellen, Punkt oder Komma. */
268 const ZAHL = String.raw`-?\d+(?:[.,]\d+)?`;
269
270 /**
271 * Liest eine Koordinate aus einer eingefuegten Zeichenkette.
272 *
273 * Erkannt werden:
274 * - "50.9412784, 6.9582814"
275 * - "50,9412784 6,9582814" (deutsches Dezimalkomma)
276 * - Google Maps "@50.9412784,6.9582814"
277 * - Apple Karten "?ll=50.9412784,6.9582814"
278 * - OpenStreetMap "#map=19/50.9412784/6.9582814"
279 * - OSM mit Marke "?mlat=50.9412784&mlon=6.9582814"
280 * - Bing Maps "?cp=50.9412784~6.9582814"
281 * - "geo:50.9412784,6.9582814"
282 *
283 * Der Weg ueber den eingefuegten Kartenverweis ist bewusst vorgesehen: Wer eine
284 * Kreuzung sucht, findet sie meist zuerst in einem Kartendienst und kann den
285 * Verweis dann einfach uebernehmen.
286 *
287 * OpenStreetMap kam zunaechst nicht vor. Fuer eine Anwendung, deren Ergebnisse
288 * in behoerdliche Verfahren gehen, ist gerade das der naheliegende Dienst - und
289 * seine Adresse traegt die Koordinate im Fragment ("#map="), das keiner der
290 * uebrigen Muster erfasste. Die Zahlen dort haben nicht zwingend
291 * Nachkommastellen: "#map=5/51.5/11" ist eine gueltige Adresse.
292 */
293 export function leseKoordinate(text: string): GeoPunkt | null {
294 const roh = text.trim();
295 if (roh === '') return null;
296
297 /*
298 * Reihenfolge mit Absicht: Erst die Muster, die ein Schluesselwort tragen,
299 * dann die beiden blanken Zahlen. Ein Verweis enthaelt viele Zahlenpaare -
300 * Zoomstufe, Bildgroesse, Zeitstempel -, und ohne Schluesselwort waere nicht
301 * zu entscheiden, welches die Koordinate ist.
302 */
303 const muster: readonly RegExp[] = [
304 // Google Maps: der Teil nach dem At-Zeichen.
305 new RegExp(String.raw`@(${ZAHL}),\s*(${ZAHL})`),
306 // OpenStreetMap und Abkoemmlinge: Zoomstufe, Breite, Laenge im Fragment.
307 new RegExp(String.raw`#map=[\d.]+/(${ZAHL})/(${ZAHL})`, 'i'),
308 // OpenStreetMap mit gesetzter Marke.
309 new RegExp(String.raw`[?&]mlat=(${ZAHL})&mlon=(${ZAHL})`, 'i'),
310 // Bing Maps: Mittelpunkt mit Tilde getrennt.
311 new RegExp(String.raw`[?&]cp=(${ZAHL})~(${ZAHL})`, 'i'),
312 // Ausdrueckliche Parameter und geo:
313 new RegExp(String.raw`(?:^geo:|[?&](?:q|ll|center)=)(${ZAHL})[,;]\s*(${ZAHL})`, 'i'),
314 // Zwei blanke Zahlen - nur, wenn sonst nichts dasteht.
315 new RegExp(String.raw`^(${ZAHL})[\s,;]+(${ZAHL})$`),
316 ];
317
318 for (const regel of muster) {
319 const treffer = regel.exec(roh);
320 if (treffer) return pruefe(zahl(treffer[1]), zahl(treffer[2]));
321 }
322
323 return null;
324 }
325
326 function zahl(wert: string | undefined): number {
327 return Number((wert ?? '').replace(',', '.'));
328 }
329
330 function pruefe(breite: number, laenge: number): GeoPunkt | null {
331 if (!Number.isFinite(breite) || !Number.isFinite(laenge)) return null;
332 if (Math.abs(breite) > 90 || Math.abs(laenge) > 180) return null;
333 return { breite, laenge };
334 }