Diese Seite enthält eine Vielzahl von Berechnungen für Längen- und Breitengrade mit den Formeln und Codefragmenten für deren Implementierung.

Alle diese Formeln dienen zur Berechnung auf der Grundlage einer kugelförmigen Erde (ohne Berücksichtigung von Ellipsoideffekten) - was für die meisten Zwecke genau genug ist ^* [[Tatsächlich ist die Erde sehr leicht ellipsoid; Die Verwendung eines sphärischen Modells führt zu Fehlern von typischerweise bis zu 0,3% ¹ - siehe Anmerkungen für weitere Einzelheiten].

Entfernung

Hierbei wird die ' Haversine' -Formel verwendet, um die Großkreisentfernung zwischen zwei Punkten zu berechnen , dh die kürzeste Entfernung über der Erdoberfläche. Dabei wird eine Luftlinie zwischen den Punkten angegeben (wobei alle Hügel, die sie fliegen, ignoriert werden) natürlich vorbei!).

var R = 6371e3 ; // Meter var ? 1 = lat1 . toRadians (); var ? 2 = lat2 . toRadians (); var ?? = ( lat2 - lat1 ). toRadians (); var ?? = ( lon2 - lon1 ). toRadians ();  
  
  
   
   

var a = Math . sin (?? / 2 ) * Math . sin (?? / 2 ) + Math . cos (? 1 ) * Math . cos (? 2 ) * Math . sin (?? / 2 ) * Math . sin (?? / 2 ); var c = 2 * Math . atan2 ( Math . sqrt    
           
          
   ( a ), Math . sqrt ( 1 - a )); 

var d = R * c ;

Beachten Sie, dass ich in diesen Skripten im Allgemeinen Lat / Lon für Breitengrad / Längengrad in Grad und ? / ? für Breiten- / Längengrad im Bogenmaß verwende - nachdem ich festgestellt habe, dass das Mischen von Grad und Bogenmaß oft der einfachste Weg ist, um Fehler am Kopf zu kratzen ...

Der Haversine Formel ¹ ‚vor besonders gut konditionierte für numerische Berechnung auch bei kleinen Abständen‘ - im Gegensatz zu Berechnungen auf der Grundlage des sphärischen Kosinussatz . Der '(re) versierte Sinus' ist 1 - cos? , und der 'halbversierte Sinus' ist (1 - cos?) / 2 oder sin² (? / 2), wie oben verwendet. Früher von Navigatoren weit verbreitet, wurde es 1984 von Roger Sinnott in der Zeitschrift Sky & Telescope („Tugenden der Haversine“) beschrieben: Sinnott erklärte, dass die Winkeltrennung zwischen Mizar und Alcor in Ursa Major - 0 ° 11'49.69 ? - möglich sei auf einem TRS-80 genau berechnet werden mit dem Haversine.

Für Neugierige ist c der Winkelabstand im Bogenmaß und a das Quadrat der halben Sehnenlänge zwischen den Punkten.

Wenn atan2nicht verfügbar, könnte c aus 2 ? asin (min (1, ? a )) berechnet werden (einschließlich Schutz vor Rundungsfehlern).

Bei Verwendung von Chrome auf einem mittelgroßen Core i5-PC dauert eine Entfernungsberechnung etwa 2 bis 5 Mikrosekunden ( daher etwa 200.000 bis 500.000 pro Sekunde). Wenig bis gar kein Nutzen wird erzielt, wenn gemeinsame Begriffe herausgerechnet werden. wahrscheinlich optimiert der JIT-Compiler sie aus.

Sphärisches Gesetz des Kosinus

Tatsächlich verwenden JavaScript (und die meisten modernen Computer und Sprachen) 64-Bit-Gleitkommazahlen 'IEEE 754', die 15 signifikante Genauigkeitsangaben liefern. Nach meiner Schätzung liefert das einfache sphärische Gesetz der Kosinusformel (cos c = cos a cos b + sin a sin b cos C ) mit dieser Genauigkeit gut konditionierte Ergebnisse bis zu Entfernungen von nur wenigen Metern auf der Erdoberfläche. ( Beachten Sie, dass die geodätische Form des Kosinusgesetzes von der kanonischen geändert wird, sodass der Breitengrad direkt und nicht der Kolatitude verwendet werden kann. )

Dies macht das einfachere Kosinusgesetz für viele geodätische Zwecke (wenn nicht für die Astronomie) zu einer vernünftigen 1-Zeilen-Alternative zur Haversinformel. Die Auswahl kann durch Programmiersprache, Prozessor, Codierungskontext, verfügbare Triggerfunktionen (in verschiedenen Sprachen) usw. bestimmt werden - und für sehr kleine Entfernungen kann eine gleichwinklige Näherung besser geeignet sein.

var ? 1 = lat1 . toRadians (), ? 2 = lat2 . toRadians (), ?? = ( lon2 - lon1 ). toRadians (), R = 6371e3 ; // gibt d in Metern var d = Math . acos ( Math . sin (? 1 ) * Math . sin (? 2 ) + Math . cos (? 1)         
    ) * Math . cos (? 2 ) * Math . cos (??) ) * R ;    

Während es einfacher ist, ist das Kosinusgesetz in meinen Tests etwas langsamer als das Haversinus.

Gleichwinklige Approximation

Wenn die Leistung ist ein Problem , und die Genauigkeit weniger wichtig, für kleine Entfernungen Satz des Pythagoras auf verwendet werden kann equirectangular Projektion : ^*

var x = (? 2 -? 1 ) * Math . cos ((? 1 + ? 2 ) / 2 ); var y = (? 2 -? 1 ); var d = Math . sqrt ( x * x + y * y ) * R ;   
 
  

Dies verwendet nur eine Trigger- und eine sqrt-Funktion - im Gegensatz zu einem halben Dutzend Triggerfunktionen für das cos-Gesetz und 7 trigs + 2 sqrts für haversine. Die Genauigkeit ist etwas komplex: Entlang der Meridiane gibt es keine Fehler, ansonsten hängen sie von Entfernung, Peilung und Breitengrad ab, sind jedoch für viele Zwecke klein genug ^* (und im Vergleich zur sphärischen Näherung selbst oft trivial).

Alternativ kann die Polarkoordinaten-Flacherdformel verwendet werden: Verwenden der Ko-Breiten ? ₁ = ? / 2 - ? ₁ und ? ₂ = ? / 2 - ? ₂ , dann d = R ? ? ? ₁ ² + ? ₂ ² - 2 ? ? ₁ ? ? ₂ ? cos ?? . Ich habe die Genauigkeit nicht verglichen.

Lager

Im Allgemeinen ändert sich Ihre aktuelle Überschrift, wenn Sie einem Großkreispfad (Orthodrom) folgen. Die endgültige Überschrift unterscheidet sich von der ursprünglichen Überschrift je nach Entfernung und Breitengrad um unterschiedliche Grade (wenn Sie von beispielsweise 35 ° N, 45 ° O (? Bagdad) zu 35 ° N, 135 ° O (? Osaka) gehen, Sie würde mit einer Überschrift von 60 ° beginnen und mit einer Überschrift von 120 ° enden!).

Diese Formel gilt für die anfängliche Peilung (manchmal auch als Vorwärtsazimut bezeichnet). Wenn Sie in einer geraden Linie entlang eines Großkreisbogens folgen, gelangen Sie vom Startpunkt zum Endpunkt: ¹

var y = Math . sin (? 2 -? 1 ) * Math . cos (? 2 ); var x = Math . cos (? 1 ) * Math . sin (? 2 ) - Math . sin (? 1 ) * Math . cos (? 2 ) * Math . cos (? 2 -? 1 ); var brng   
  
        
= Math . atan2 ( y , x ). toDegrees (); 

Da atan2 Werte im Bereich -? ... + ? (dh -180 ° ... + 180 °) zurückgibt, um das Ergebnis auf eine Kompasspeilung (im Bereich 0 ° ... 360 °, mit zu normalisieren) ?ve Werte in den Bereich 180 ° ... 360 ° transformiert), in Grad konvertieren und dann (? + 360)% 360 verwenden, wobei% (Gleitkomma-) Modulo ist.

Zur endgültigen Lagerung nimmt einfach den Anfang Lager von dem Ende Punkt zum Startpunkt und Rückwärts es (mit ? = (? + 180)% 360).

Mittelpunkt

Dies ist der halbe Weg entlang eines Großkreiswegs zwischen den beiden Punkten. ¹

var Bx = Math . cos (? 2 ) * Math . cos (? 2 -? 1 ); var By = Math . cos (? 2 ) * Math . sin (? 2 -? 1 ); var ? 3 = Math . atan2 ( Math . sin (? 1 ) + Math . sin (? 2 ),     
     
     
                    Mathe . sqrt ( ( Math . cos (? 1 ) + Bx ) * ( Math . cos (? 1 ) + Bx ) + By * By ) ); var ? 3 = ? 1 + Math . atan2 ( By , Math . cos (? 1 ) + Bx );     
        

Ebenso wie die anfängliche Peilung von der endgültigen Peilung abweichen kann, befindet sich der Mittelpunkt möglicherweise nicht auf halber Strecke zwischen Breiten- und Längengraden. Der Mittelpunkt zwischen 35 ° N, 45 ° O und 35 ° N, 135 ° E liegt bei 45 ° N, 90 ° E.

Zwischenpunkt

Ein Zwischenpunkt an einem beliebigen Bruchteil entlang des Großkreiswegs zwischen zwei Punkten kann ebenfalls berechnet werden. ¹

Zielpunkt bei gegebener Entfernung und Peilung vom Startpunkt

Bei gegebenem Startpunkt, anfänglicher Peilung und Entfernung werden der Zielpunkt und die endgültige Peilung berechnet, die sich entlang eines (kürzesten Entfernungs-) Großkreisbogens bewegen.

var ? 2 = Math . asin ( Math . sin (? 1 ) * Math . cos ( d / R ) + Math . cos (? 1 ) * Math . sin ( d / R ) * Math . cos ( brng ) ); var ? 2 = ? 1 + Math . atan2 ( Math .     
                     
     sin ( brng ) * Math . sin ( d / R ) * Math . cos (? 1 ), Math . cos ( d / R ) - Math . sin (? 1 ) * Math . sin (? 2 ));
                         

Für die endgültigen Lager, nehmen Sie einfach die Anfangs Peilung vom Ende Punkt zum Startpunkt und umgekehrt mit (brng+180)%360.

Schnittpunkt zweier Pfade mit Startpunkten und Peilungen

Dies ist eine etwas komplexere Berechnung als die meisten anderen auf dieser Seite, aber ich wurde mehrmals danach gefragt. Dies stammt aus Ed Williams Luftfahrtformel . Siehe unten für das JavaScript.

Dies ist viel einfacher, wenn Vektoren anstelle der sphärischen Trigonometrie verwendet werden: siehe latlong-vectors.html .

Streckenübergreifende Entfernung

Hier ist eine neue: Ich wurde manchmal nach der Entfernung eines Punktes von einem Großkreispfad gefragt (manchmal als Cross-Track-Fehler bezeichnet).

var ? 13 = d13 / R ; var dXt = Math . asin ( Math . sin (? 13 ) * Math . sin (? 13 -? 12 )) * R ;  
  

Hier wird der Großkreispfad durch einen Startpunkt und einen Endpunkt identifiziert. Abhängig von den Anfangsdaten, von denen aus Sie arbeiten, können Sie die obigen Formeln verwenden, um den relevanten Abstand und die Peilungen zu erhalten. Das Vorzeichen von d _xt gibt an, auf welcher Seite des Pfades sich der dritte Punkt befindet.

Die Entfernung entlang der Strecke vom Startpunkt zum nächstgelegenen Punkt auf dem Weg zum dritten Punkt beträgt

var ? 13 = d13 / R ; var dAt = Math . acos ( Math . cos (? 13 ) / Math . cos ( dXt / R )) * R ;  
  

Nächster Punkt zu den Polen

Und: 'Clairauts Formel' gibt Ihnen den maximalen Breitengrad eines Großkreiswegs bei gegebener Peilung ? und Breitengrad ? auf dem Großkreis:

var ? Max = Math . acos ( Math . abs ( Math . sin (?) * Math . cos (?)));   

Rhumb Lines

Eine "Loxodrome" (oder Loxodrom) ist ein Pfad konstanter Peilung, der alle Meridiane im gleichen Winkel kreuzt.

Segler pflegten (und manchmal immer noch) entlang der Loxodrome zu navigieren, da es einfacher ist, einer konstanten Kompasspeilung zu folgen, als die Peilung kontinuierlich anzupassen, wie dies erforderlich ist, um einem großen Kreis zu folgen. Rhumb-Linien sind gerade Linien auf einer Mercator-Projektionskarte (auch hilfreich für die Navigation).

Rhumbuslinien sind im Allgemeinen länger als Großkreisrouten (Orthodromrouten). Zum Beispiel ist London nach New York entlang einer Loxodrome 4% länger als entlang eines großen Kreises - wichtig für Flugbenzin, aber nicht besonders für Segelschiffe. New York nach Peking - nahe am extremsten Beispiel (obwohl nicht segelbar!) - ist entlang einer Loxodrome 30% länger.

Schlüssel zur Berechnung von rhumb Linien ist die inverse Gudermannfunktion ¹ , die die Höhe an einer Mercator - Projektion Karte eines bestimmten Breitengrad ergibt: ln (tan & PHgr; + sec?) oder ln (tan (? / 4 + ? / 2)) . Dies neigt natürlich zu Unendlichkeit an den Polen (in Übereinstimmung mit der Mercator-Projektion). Für Obsessive gibt es sogar eine ellipsoide Version, den 'isometrischen Breitengrad': ? = ln (tan (? / 4 + ? / 2) / [(1 - e?sin?) / (1 + e?sin?)] ^{e / 2} ) oder sein besser konditioniertes Äquivalent ? = atanh (sin?) - e?atanh (e?sin?) .

Die Formeln zur Ableitung der Mercator-Projektions-Ost- und Nordkoordinaten aus dem sphärischen Breiten- und Längengrad sind dann ¹

Die folgenden Formeln stammen aus der Luftfahrtformel von Ed Williams. ¹

Entfernung

Da eine Loxodrome eine gerade Linie auf einer Mercator-Projektion ist, entspricht der Abstand zwischen zwei Punkten entlang einer Loxodrome der Länge dieser Linie (nach Pythagoras). Die Verzerrung der Projektion muss jedoch ausgeglichen werden.

Auf einem konstanten Breitengradkurs (von Ost nach West) ist diese Kompensation einfach cos? ; Im allgemeinen Fall ist es ?? / ?? wo ?? = ln (tan (? / 4 + ? ₂ /2) / tan (? / 4 + ? ₁ /2)) (die 'projizierte' Breitengrad - Differenz)

var ?? = Math . log ( Math . tan ( Math . PI / 4 + ? 2 / 2 ) / Math . tan ( Math . PI / 4 + ? 1 / 2 )); var q = Math . abs (??) > 10e-12 ? ?? / ?? : Math . cos (? 1 ); // EW-Kurs wird mit 0/0 schlecht konditioniert   
        

// Wenn dLon über 180 ° ist, nehmen Sie eine kürzere Loxodrome über den Anti-Meridian: wenn ( Math . abs (??) > Math . PI ) ?? = ??> 0 ? - ( 2 * Math . PI- ??) : ( 2 * Math . PI + ??);
          

var dist = Math . sqrt (?? * ?? + q * q * ?? * ??) * R ;   

Lager

Eine Loxodrome ist eine gerade Linie auf einer Mercator-Projektion mit einem Winkel auf der Projektion, der der Kompasspeilung entspricht.

var ?? = Math . log ( Math . tan ( Math . PI / 4 + ? 2 / 2 ) / Math . tan ( Math . PI / 4 + ? 1 / 2 ));   

// Wenn dLon über 180 ° ist, nehmen Sie eine kürzere Loxodrome über den Anti-Meridian: wenn ( Math . abs (??) > Math . PI ) ?? = ??> 0 ? - ( 2 * Math . PI- ??) : ( 2 * Math . PI + ??);
          

var brng = Math . atan2 (??, ??). toDegrees ();  

Ziel

Bei einem Startpunkt und einem Abstand d entlang der konstanten Peilung ? wird der Zielpunkt berechnet. Wenn Sie eine konstante Peilung entlang einer Loxodrome beibehalten, werden Sie sich allmählich in Richtung eines der Pole drehen.

var ? = d / R ; var ?? = ? * Math . cos (?); var ? 2 = ? 1 + ??;  
     
     

var ?? = Math . log ( Math . tan (? 2 / 2 + Math . PI / 4 ) / Math . tan (? 1 / 2 + Math . PI / 4 )); var q = Math . abs (??) > 10e-12 ? ?? / ?? : Math . cos (? 1 );   
          // EW-Kurs wird mit 0/0 schlecht konditioniert

var ?? = ? * Math . sin (?) / q ; var ? 2 = ? 1 + ??;   
     

// Überprüfen Sie, ob ein dummer Mistkerl am Pol vorbeigeht, und normalisieren Sie den Breitengrad, wenn ja, wenn ( Math . abs (? 2 ) > Math . PI / 2 ) ? 2 = ? 2 > 0 ? Mathe . PI -? 2 : - Math . PI -? 2 ;
          

Mittelpunkt

Diese Formel zur Berechnung des 'loxodromen Mittelpunkts', des Punktes auf halber Strecke entlang einer Loxodrome zwischen zwei Punkten, ist Robert Hill und Clive Tooth ¹ (thx Axel!) Zu verdanken.

if ( Math . abs (? 2 -? 1 ) > Math . PI ) ? 1 + = 2 * Math . PI ; // Anti-Meridian kreuzen       

var ? 3 = (? 1 + ? 2 ) / 2 ; var f1 = Math . tan ( Math . PI / 4 + ? 1 / 2 ); var f2 = Math . tan ( Math . PI / 4 + ? 2 / 2 ); var f3 = Math . tan ( Math . PI   
   
   
 / 4 + ? 3 / 2 ); var ? 3 = ( (? 2 -? 1 ) * Math . log ( f3 ) + ? 1 * Math . log ( f2 ) - ? 2 * Math . log ( f1 ) ) / Math . log ( f2 / f1 );  
           

wenn (! isFinite (? 3 )) ? 3 = (? 1 + ? 2 ) / 2 ; // Breitengrad parallel     

Verwenden der Skripte auf Webseiten

Die Verwendung dieser Skripte auf Webseiten wäre ungefähr so:

<! doctype html> <html lang = "en" > <head> <title> Verwenden der Skripte auf Webseiten </ title> <meta charset = "utf-8" > <script type = "module" > importieren Sie LatLon aus 'https://cdn.jsdelivr.net/npm/geodesy@2/latlon-spherical.min.js' ; 
        Dokument . addEventListener ( 'DOMContentLoaded' , function () { 
            document . querySelector ( '
 

    
     
     
           onclick = function () { 
                berechneDistanz (); } }); Funktion berechneDistanz () { const p1 = LatLon . parse ( Dokument . querySelector ( '# point1' ). Wert ); const p2 = LatLon . parse ( Dokument . querySelector ( '# point2' ). Wert ); const dist = parseFloat  
            
        
         
             
             
            ( p1 . distanceTo ( p2 ). toPrecision ( 4 )); 
            Dokument . querySelector ( '# result-distance' ). textContent = dist + ' meter ' ; } </ script> </ head> <body> <form> 
    Punkt 1: <input type = "text" name = "point1" id = "point1" placeholder = "lat1, lon1" > 
    Punkt 2: < 
        
    


     = "text" name = "point2" id = "point2" placeholder = "lat2, lon2" > <button type = "button" id = "calc-dist" > Abstand berechnen </ button> <output id = "result- Entfernung " > </ output> </ form> </ body> </ html>   
      
     

Konvertieren Sie zwischen Grad-Minuten-Sekunden und Dezimalgrad

Anmerkungen:

• Genauigkeit: Da die Erde keine Kugel ist, gibt es kleine Fehler bei der Verwendung der sphärischen Geometrie. Die Erde ist tatsächlich ungefähr ellipsoid (oder genauer gesagt, abgeflachtes Sphäroid) mit einem Radius zwischen etwa 6.378 km (äquatorial) und 6.357 km (polar) und einem lokalen Krümmungsradius zwischen 6.336 km (äquatorialer Meridian) und 6.399 km (polar) ). 6.371 km ist der allgemein akzeptierte Wert für den mittleren Radius der Erde . Dies bedeutet, dass die Fehler bei der Annahme einer sphärischen Geometrie bis zu 0,55% betragen können, wenn sie den Äquator überqueren, jedoch im Allgemeinen unter 0,3%, abhängig von Breitengrad und Fahrtrichtung ( Whuber untersucht dies beim Stapelaustausch ausführlich). Eine Genauigkeit von besser als 3 m in 1 km ist für mich meistens gut genug, aber wenn Sie eine höhere Genauigkeit wünschen, können Sie die Vincenty- Formel zur Berechnung geodätischer Entfernungen auf Ellipsoiden verwenden, die Ergebnisse mit einer Genauigkeit von bis zu 1 mm liefert. (Aus purer Perversität - ich habe noch nie eine solche Genauigkeit benötigt - habe ich diese Formel nachgeschlagen und festgestellt, dass die JavaScript-Implementierung einfacher war als erwartet).
• Triggerfunktionen verwenden Argumente im Bogenmaß , daher müssen Breite, Länge und Peilung in Grad (entweder Dezimal oder Grad / Minuten / Sekunden) in Bogenmaß konvertiert werden, rad = ?.deg / 180. Bei der Umrechnung des Bogenmaßes in Grad (Grad = 180.rad / ?) ist West negativ, wenn vorzeichenbehaftete Dezimalgrade verwendet werden. Für Lager müssen Werte im Bereich von -? bis + ? [-180 ° bis + 180 °] in 0 bis + 2? [0 ° –360 °] umgewandelt werden. Dies kann durch (brng + 2.?)% 2.? [oder brng + 360)% 360] erfolgen, wobei% der (Gleitkomma-) Modulooperator ist (beachten Sie, dass verschiedene Sprachen die Modulooperation auf unterschiedliche Weise implementieren ).
• Alle Lager beziehen sich auf den wahren Norden , 0 ° = N, 90 ° = E usw.; Wenn Sie mit einem Kompass arbeiten, variiert der magnetische Norden auf komplexe Weise um die Erde vom wahren Norden, und der Unterschied muss durch Abweichungen ausgeglichen werden, die auf lokalen Karten angegeben sind.
• Die hier weit verbreitete Funktion atan2 () verwendet zwei Argumente, atan2 (y, x), und berechnet den Bogentangens des Verhältnisses y / x. Es ist flexibler als atan (y / x), da es x = 0 verarbeitet und auch Werte in allen 4 Quadranten -? bis + ? zurückgibt (die atan-Funktion gibt Werte im Bereich -? / 2 bis + ? / zurück) 2).
• Wenn Sie eine Formel mit atan2 in einer Tabelle implementieren (Microsoft Excel , LibreOffice Calc, Google Sheets, Apple Numbers), müssen Sie die Argumente umkehren, da Excel usw. sie umgekehrt zu JavaScript verwendet - die herkömmliche Reihenfolge ist atan2 (y) , x), aber Excel verwendet atan2 (x, y). Um atan2 in einem (VBA) -Makro zu verwenden, können Sie WorksheetFunction.Atan2 () verwenden.
• Wenn Sie Google Maps verwenden , werden jetzt mehrere dieser Funktionen in der "sphärischen" Bibliothek der Google Maps API V3 (computeDistanceBetween (), computeHeading (), computeOffset (), interpolate () usw. bereitgestellt. Beachten Sie, dass sie einen Standard-Erdradius verwenden von 6.378.137 Metern).
• Wenn Sie UK Ordnance Survey Grid-Referenzen verwenden, habe ich ein Skript zum Konvertieren zwischen Lat / Long- und OS-Grid-Referenzen implementiert .
• Wenn Sie UTM-Koordinaten oder MGRS-Rasterreferenzen verwenden, habe ich Skripts für die Konvertierung zwischen Lat / Long, UTM und MGRS implementiert .
• Ich habe viel aus den GIS-FAQ des US Census Bureau gelernt, die nicht mehr verfügbar sind, also habe ich eine Kopie erstellt.
• Vielen Dank an Ed Williams ' Aviation Formulary für viele der Formeln.
• Teilen Sie km für Meilen durch 1.609344
• Teilen Sie für Seemeilen km durch 1,852

Weiter unten finden Sie den JavaScript-Quellcode, der auch auf GitHub verfügbar ist . Eine vollständige Dokumentation sowie eine Testsuite sind verfügbar .

Hinweis Ich verwende griechische Buchstaben in Variablen, die mathematische Symbole darstellen, die üblicherweise als griechische Buchstaben dargestellt werden: Ich schätze den großen Vorteil der Lesbarkeit gegenüber den geringfügigen Unannehmlichkeiten beim Tippen (wenn Sie auf Probleme stoßen, stellen Sie sicher, dass Ihre <head>Includes enthalten sind <meta charset="utf-8">) und verwende beim Speichern von Dateien die UTF-8-Codierung ).

Mit seiner untypisierten C-Syntax liest sich JavaScript bemerkenswert ähnlich wie Pseudocode: Die Algorithmen werden mit einem Minimum an syntaktischen Ablenkungen verfügbar gemacht. Diese Funktionen sollten bei Bedarf einfach in andere Sprachen zu übersetzen sein, können jedoch auch unverändert in Browsern und Node.js verwendet werden.

Ich habe das Basis-JavaScript-Nummernobjekt mit den Methoden toRadians () und toDegrees () erweitert: Ich sehe keine große Wahrscheinlichkeit von Konflikten, da dies allgegenwärtige Operationen sind.

Ich habe auch eine Seite, die die Verwendung des sphärischen Kosinusgesetzes zur Auswahl von Punkten aus einer Datenbank innerhalb eines bestimmten Begrenzungskreises veranschaulicht. Das Beispiel basiert auf MySQL + PDO, sollte jedoch auf andere DBMS-Plattformen erweiterbar sein.

Mehrere Personen haben nach Beispiel- Excel- Tabellen gefragt , daher habe ich die Formeln für Entfernung und Peilung sowie den Zielpunkt als Tabellen in einer Form implementiert , in der alle beteiligten Phasen zur Veranschaulichung des Vorgangs aufgeschlüsselt sind.

Februar 2019 : Ich habe die Bibliothek überarbeitet, um ES-Module zu verwenden und ihren Umfang zu erweitern. Weitere Informationen finden Sie in der GitHub README und im CHANGELOG.

Leistung: Wie oben erwähnt, dauert die Berechnung des Haversinusabstands etwa 2 bis 5 Mikrosekunden ( daher etwa 200.000 bis 500.000 pro Sekunde). Ich habe noch keine Timing-Tests für andere Berechnungen durchgeführt.

Andere Sprachen : Ich kann keine Übersetzungen in andere Sprachen unterstützen. Wenn Sie den Code jedoch in eine andere Sprache portiert haben, kann ich Ihnen hier gerne Links zur Verfügung stellen.

• Brian Lambert hat eine Objective-C- Version erstellt.
• Jean Brouwers hat eine Python- Version erstellt.
• Vahan Aghajanyan hat eine C ++ - Version erstellt.

   * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * / / * Sphärische Geodäsie-Tools für Breite / Länge (c) Chris Veness 2002-2019 * / / * MIT-Lizenz * / / * www.movable-type.co.uk/scripts/latlong.html * / / * www.movable- type.co.uk/scripts/geodesy-library.html#latlon-spherical * / / * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * /






import Dms von ‘./dms.js' ;   

const ? = Math . PI ;   


/ **
 * Bibliothek von Geodäsiefunktionen für Operationen an einem sphärischen Erdmodell.
 * *
 * Beinhaltet Entfernungen, Peilungen, Ziele usw. sowohl für Großkreispfade als auch für Loxodrome.
 * und andere verwandte Funktionen.
 * *
 * Alle Berechnungen werden mit einfachen sphärischen trigonometrischen Formeln durchgeführt.
 * *
 * @module latlon-sphärisch
 * /

// Beachten Sie, dass griechische Buchstaben (z. B. ?, ?, ?) für Winkel im Bogenmaß verwendet werden, um von Winkeln in // Grad (z. B. lat, lon, brng) zu unterscheiden.



/ * LatLonSpherical - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * /


/ **
 * Breiten- / Längengrade auf einer sphärischen Modellerde und Methoden zur Berechnung von Entfernungen,
 * Peilungen, Ziele usw. auf (orthodromen) Großkreispfaden und (loxodromen) Loxodromen.
 * / class LatLonSpherical {
  

    / **
     * Erstellt einen Breiten- / Längengradpunkt auf der Erdoberfläche unter Verwendung eines sphärischen Modells Erde.
     * *
     * @param {number} lat - Breitengrad (in Grad).
     * @param {number} lon - Länge (in Grad).
     * @throws {TypeError} Ungültiges lat / lon.
     * *
     * @Beispiel
     * LatLon aus '/js/geodesy/latlon-spherical.js' importieren;
     * const p = neues LatLon (52,205, 0,119);
     * / 
    Konstruktor ( lat , lon ) { if ( isNaN ( lat )) wirft neuen TypeError ( `ungültiges lat '$ {lat}'` ); if ( isNaN ( lon )) wirft einen neuen TypeError ( `ungültiges lon '$ {lon}'` ); 
            
            

        das . _lat = Dms . wrap90 ( lat ); das . _lon = Dms . wrap180 ( lon ); }} 
         
    


    / **
     * Breitengrad in Grad nördlich vom Äquator (einschließlich Aliase Lat, Breitengrad): kann eingestellt werden als
     * numerisch oder hexagesimal (Grad-Min-Sek); als numerisch zurückgegeben.
     */
    get lat()       { return this._lat; }
    get latitude()  { return this._lat; }
    set lat(lat) {
        this._lat = isNaN(lat) ? Dms.wrap90(Dms.parse(lat)) : Dms.wrap90(lat);
        if (isNaN(this._lat)) throw new TypeError(`invalid lat ‘${lat}’`);
    }
    set latitude(lat) {
        this._lat = isNaN(lat) ? Dms.wrap90(Dms.parse(lat)) : Dms.wrap90(lat);
        if (isNaN(this._lat)) throw new TypeError(`invalid latitude ‘${lat}’`);
    }

    /**
     * Longitude in degrees east from international reference meridian (including aliases lon, lng,
     * longitude): can be set as numeric or hexagesimal (deg-min-sec); returned as numeric.
     */
    get lon()       { return this._lon; }
    get lng()       { return this._lon; }
    get longitude() { return this._lon; }
    set lon(lon) {
        this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(lon);
        if (isNaN(this._lon)) throw new TypeError(`invalid lon ‘${lon}’`);
    }
    set lng(lon) {
        this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(lon);
        if (isNaN(this._lon)) throw new TypeError(`invalid lng ‘${lon}’`);
    }
    set longitude(lon) {
        this._lon = isNaN(lon) ? Dms.wrap180(Dms.parse(lon)) : Dms.wrap180(lon);
        if (isNaN(this._lon)) throw new TypeError(`invalid longitude ‘${lon}’`);
    }


    /** Conversion factors; 1000 * LatLon.metresToKm gives 1. */
    static get metresToKm()            { return 1/1000; }
    /** Conversion factors; 1000 * LatLon.metresToMiles gives 0.621371192237334. */
    static get metresToMiles()         { return 1/1609.344; }
    /** Conversion factors; 1000 * LatLon.metresToMiles gives 0.5399568034557236. */
    static get metresToNauticalMiles() { return 1/1852; }


    /**
     * Parses a latitude/longitude point from a variety of formats.
     *
     * Latitude & longitude (in degrees) can be supplied as two separate parameters, as a single
     * comma-separated lat/lon string, or as a single object with { lat, lon } or GeoJSON properties.
     *
     * The latitude/longitude values may be numeric or strings; they may be signed decimal or
     * deg-min-sec (hexagesimal) suffixed by compass direction (NSEW); a variety of separators are
     * accepted. Examples -3.62, '3 37 12W', '3°37?12?W'.
     *
     * Thousands/decimal separators must be comma/dot; use Dms.fromLocale to convert locale-specific
     * thousands/decimal separators.
     *
     * @param   {number|string|Object} lat|latlon - Latitude (in degrees) or comma-separated lat/lon or lat/lon object.
     * @param   {number|string}        [lon]      - Longitude (in degrees).
     * @returns {LatLon} Latitude/longitude point.
     * @throws  {TypeError} Invalid point.
     *
     * @example
     *   const p1 = LatLon.parse(52.205, 0.119);                                    // numeric pair (? new LatLon)
     *   const p2 = LatLon.parse('52.205', '0.119');                                // numeric string pair (? new LatLon)
     *   const p3 = LatLon.parse('52.205, 0.119');                                  // single string numerics
     *   const p4 = LatLon.parse('52°12?18.0?N', '000°07?08.4?E');                  // DMS pair
     *   const p5 = LatLon.parse('52°12?18.0?N, 000°07?08.4?E');                    // single string DMS
     *   const p6 = LatLon.parse({ lat: 52.205, lon: 0.119 });                      // { lat, lon } object numeric
     *   const p7 = LatLon.parse({ lat: '52°12?18.0?N', lng: '000°07?08.4?E' });    // { lat, lng } object DMS
     *   const p8 = LatLon.parse({ type: 'Point', coordinates: [ 0.119, 52.205] }); // GeoJSON
     */
    static parse(...args) {
        if (args.length == 0) throw new TypeError('invalid (empty) point');
        if (args[0]===null || args[1]===null) throw new TypeError('invalid (null) point');

        let lat=undefined, lon=undefined;

        if (args.length == 2) { // regular (lat, lon) arguments
            [ lat, lon ] = args;
            lat = Dms.wrap90(Dms.parse(lat));
            lon = Dms.wrap180(Dms.parse(lon));
            if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${args.toString()}’`);
        }

        if (args.length == 1 && typeof args[0] == 'string') { // single comma-separated lat,lon string
            [ lat, lon ] = args[0].split(',');
            lat = Dms.wrap90(Dms.parse(lat));
            lon = Dms.wrap180(Dms.parse(lon));
            if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${args[0]}’`);
        }

        if (args.length == 1 && typeof args[0] == 'object') { // single { lat, lon } object
            const ll = args[0];
            if (ll.type == 'Point' && Array.isArray(ll.coordinates)) { // GeoJSON
                [ lon, lat ] = ll.coordinates;
            } else { // regular { lat, lon } object
                if (ll.latitude  != undefined) lat = ll.latitude;
                if (ll.lat       != undefined) lat = ll.lat;
                if (ll.longitude != undefined) lon = ll.longitude;
                if (ll.lng       != undefined) lon = ll.lng;
                if (ll.lon       != undefined) lon = ll.lon;
                lat = Dms.wrap90(Dms.parse(lat));
                lon = Dms.wrap180(Dms.parse(lon));
            }
            if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${JSON.stringify(args[0])}’`);
        }

        if (isNaN(lat) || isNaN(lon)) throw new TypeError(`invalid point ‘${args.toString()}’`);

        return new LatLonSpherical(lat, lon);
    }


    /**
     * Returns the distance along the surface of the earth from ‘this’ point to destination point.
     *
     * Uses haversine formula: a = sin²(??/2) + cos?1·cos?2 · sin²(??/2); d = 2 · atan2(?a, ?(a-1)).
     *
     * @param   {LatLon} point - Latitude/longitude of destination point.
     * @param   {number} [radius=6371e3] - Radius of earth (defaults to mean radius in metres).
     * @returns {number} Distance between this point and destination point, in same units as radius.
     * @throws  {TypeError} Invalid radius.
     *
     * @example
     *   const p1 = new LatLon(52.205, 0.119);
     *   const p2 = new LatLon(48.857, 2.351);
     *   const d = p1.distanceTo(p2);       // 404.3×10³ m
     *   const m = p1.distanceTo(p2, 3959); // 251.2 miles
     */
    distanceTo(point, radius=6371e3) {
        if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms
        if (isNaN(radius)) throw new TypeError(`invalid radius ‘${radius}’`);

        // a = sin²(??/2) + cos(?1)?cos(?2)?sin²(??/2)
        // ? = 2·atan2(?(a), ?(1?a))
        // see mathforum.org/library/drmath/view/51879.html for derivation

        const R = radius;
        const ?1 = this.lat.toRadians(),  ?1 = this.lon.toRadians();
        const ?2 = point.lat.toRadians(), ?2 = point.lon.toRadians();
        const ?? = ?2 - ?1;
        const ?? = ?2 - ?1;

        const a = Math.sin(??/2)*Math.sin(??/2) + Math.cos(?1)*Math.cos(?2) * Math.sin(??/2)*Math.sin(??/2);
        const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
        const d = R * c;

        return d;
    }


    /**
     * Returns the initial bearing from ‘this’ point to destination point.
     *
     * @param   {LatLon} point - Latitude/longitude of destination point.
     * @returns {number} Initial bearing in degrees from north (0°..360°).
     *
     * @example
     *   const p1 = new LatLon(52.205, 0.119);
     *   const p2 = new LatLon(48.857, 2.351);
     *   const b1 = p1.initialBearingTo(p2); // 156.2°
     */
    initialBearingTo(point) {
        if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms
        if (this.equals(point)) return NaN; // coincident points

        // tan? = sin???cos?2 / cos?1?sin?2 ? sin?1?cos?2?cos??
        // see mathforum.org/library/drmath/view/55417.html for derivation

        const ?1 = this.lat.toRadians();
        const ?2 = point.lat.toRadians();
        const ?? = (point.lon - this.lon).toRadians();

        const x = Math.cos(?1) * Math.sin(?2) - Math.sin(?1) * Math.cos(?2) * Math.cos(??);
        const y = Math.sin(??) * Math.cos(?2);
        const ? = Math.atan2(y, x);

        const bearing = ?.toDegrees();

        return Dms.wrap360(bearing);
    }


    /**
     * Returns final bearing arriving at destination point from ‘this’ point; the final bearing will
     * differ from the initial bearing by varying degrees according to distance and latitude.
     *
     * @param   {LatLon} point - Latitude/longitude of destination point.
     * @returns {number} Final bearing in degrees from north (0°..360°).
     *
     * @example
     *   const p1 = new LatLon(52.205, 0.119);
     *   const p2 = new LatLon(48.857, 2.351);
     *   const b2 = p1.finalBearingTo(p2); // 157.9°
     */
    finalBearingTo(point) {
        if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms

        // get initial bearing from destination point to this point & reverse it by adding 180°

        const bearing = point.initialBearingTo(this) + 180;

        return Dms.wrap360(bearing);
    }


    /**
     * Returns the midpoint between ‘this’ point and destination point.
     *
     * @param   {LatLon} point - Latitude/longitude of destination point.
     * @returns {LatLon} Midpoint between this point and destination point.
     *
     * @example
     *   const p1 = new LatLon(52.205, 0.119);
     *   const p2 = new LatLon(48.857, 2.351);
     *   const pMid = p1.midpointTo(p2); // 50.5363°N, 001.2746°E
     */
    midpointTo(point) {
        if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms

        // ?m = atan2( sin?1 + sin?2, ?( (cos?1 + cos?2?cos??)² + cos²?2?sin²?? ) )
        // ?m = ?1 + atan2(cos?2?sin??, cos?1 + cos?2?cos??)
        // midpoint is sum of vectors to two points: mathforum.org/library/drmath/view/51822.html

        const ?1 = this.lat.toRadians();
        const ?1 = this.lon.toRadians();
        const ?2 = point.lat.toRadians();
        const ?? = (point.lon - this.lon).toRadians();

        // get cartesian coordinates for the two points
        const A = { x: Math.cos(?1), y: 0, z: Math.sin(?1) }; // place point A on prime meridian y=0
        const B = { x: Math.cos(?2)*Math.cos(??), y: Math.cos(?2)*Math.sin(??), z: Math.sin(?2) };

        // vector to midpoint is sum of vectors to two points (no need to normalise)
        const C = { x: A.x + B.x, y: A.y + B.y, z: A.z + B.z };

        const ?m = Math.atan2(C.z, Math.sqrt(C.x*C.x + C.y*C.y));
        const ?m = ?1 + Math.atan2(C.y, C.x);

        const lat = ?m.toDegrees();
        const lon = ?m.toDegrees();

        return new LatLonSpherical(lat, lon);
    }


    /**
     * Returns the point at given fraction between ‘this’ point and given point.
     *
     * @param   {LatLon} point - Latitude/longitude of destination point.
     * @param   {number} fraction - Fraction between the two points (0 = this point, 1 = specified point).
     * @returns {LatLon} Intermediate point between this point and destination point.
     *
     * @example
     *   const p1 = new LatLon(52.205, 0.119);
     *   const p2 = new LatLon(48.857, 2.351);
     *   const pInt = p1.intermediatePointTo(p2, 0.25); // 51.3721°N, 000.7073°E
     */
    intermediatePointTo(point, fraction) {
        if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms
        if (this.equals(point)) return new LatLonSpherical(this.lat, this.lon); // coincident points

        const ?1 = this.lat.toRadians(), ?1 = this.lon.toRadians();
        const ?2 = point.lat.toRadians(), ?2 = point.lon.toRadians();

        // distance between points
        const ?? = ?2 - ?1;
        const ?? = ?2 - ?1;
        const a = Math.sin(??/2) * Math.sin(??/2)
            + Math.cos(?1) * Math.cos(?2) * Math.sin(??/2) * Math.sin(??/2);
        const ? = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

        const A = Math.sin((1-fraction)*?) / Math.sin(?);
        const B = Math.sin(fraction*?) / Math.sin(?);

        const x = A * Math.cos(?1) * Math.cos(?1) + B * Math.cos(?2) * Math.cos(?2);
        const y = A * Math.cos(?1) * Math.sin(?1) + B * Math.cos(?2) * Math.sin(?2);
        const z = A * Math.sin(?1) + B * Math.sin(?2);

        const ?3 = Math.atan2(z, Math.sqrt(x*x + y*y));
        const ?3 = Math.atan2(y, x);

        const lat = ?3.toDegrees();
        const lon = ?3.toDegrees();

        return new LatLonSpherical(lat, lon);
    }


    /**
     * Returns the destination point from ‘this’ point having travelled the given distance on the
     * given initial bearing (bearing normally varies around path followed).
     *
     * @param   {number} distance - Distance travelled, in same units as earth radius (default: metres).
     * @param   {number} bearing - Initial bearing in degrees from north.
     * @param   {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
     * @returns {LatLon} Destination point.
     *
     * @example
     *   const p1 = new LatLon(51.47788, -0.00147);
     *   const p2 = p1.destinationPoint(7794, 300.7); // 51.5136°N, 000.0983°W
     */
    destinationPoint(distance, bearing, radius=6371e3) {
        // sin?2 = sin?1?cos? + cos?1?sin??cos?
        // tan?? = sin??sin??cos?1 / cos??sin?1?sin?2
        // see mathforum.org/library/drmath/view/52049.html for derivation

        const ? = distance / radius; // angular distance in radians
        const ? = Number(bearing).toRadians();

        const ?1 = this.lat.toRadians(), ?1 = this.lon.toRadians();

        const sin?2 = Math.sin(?1) * Math.cos(?) + Math.cos(?1) * Math.sin(?) * Math.cos(?);
        const ?2 = Math.asin(sin?2);
        const y = Math.sin(?) * Math.sin(?) * Math.cos(?1);
        const x = Math.cos(?) - Math.sin(?1) * sin?2;
        const ?2 = ?1 + Math.atan2(y, x);

        const lat = ?2.toDegrees();
        const lon = ?2.toDegrees();

        return new LatLonSpherical(lat, lon);
    }


    /**
     * Returns the point of intersection of two paths defined by point and bearing.
     *
     * @param   {LatLon}      p1 - First point.
     * @param   {number}      brng1 - Initial bearing from first point.
     * @param   {LatLon}      p2 - Second point.
     * @param   {number}      brng2 - Initial bearing from second point.
     * @returns {LatLon|null} Destination point (null if no unique intersection defined).
     *
     * @example
     *   const p1 = new LatLon(51.8853, 0.2545), brng1 = 108.547;
     *   const p2 = new LatLon(49.0034, 2.5735), brng2 =  32.435;
     *   const pInt = LatLon.intersection(p1, brng1, p2, brng2); // 50.9078°N, 004.5084°E
     */
    static intersection(p1, brng1, p2, brng2) {
        if (!(p1 instanceof LatLonSpherical)) p1 = LatLonSpherical.parse(p1); // allow literal forms
        if (!(p2 instanceof LatLonSpherical)) p2 = LatLonSpherical.parse(p2); // allow literal forms
        if (isNaN(brng1)) throw new TypeError(`invalid brng1 ‘${brng1}’`);
        if (isNaN(brng2)) throw new TypeError(`invalid brng2 ‘${brng2}’`);

        // see www.edwilliams.org/avform.htm#Intersection

        const ?1 = p1.lat.toRadians(), ?1 = p1.lon.toRadians();
        const ?2 = p2.lat.toRadians(), ?2 = p2.lon.toRadians();
        const ?13 = Number(brng1).toRadians(), ?23 = Number(brng2).toRadians();
        const ?? = ?2 - ?1, ?? = ?2 - ?1;

        // angular distance p1-p2
        const ?12 = 2 * Math.asin(Math.sqrt(Math.sin(??/2) * Math.sin(??/2)
            + Math.cos(?1) * Math.cos(?2) * Math.sin(??/2) * Math.sin(??/2)));
        if (Math.abs(?12) < Number.EPSILON) return new LatLonSpherical(p1.lat, p1.lon); // coincident points

        // initial/final bearings between points
        const cos?a = (Math.sin(?2) - Math.sin(?1)*Math.cos(?12)) / (Math.sin(?12)*Math.cos(?1));
        const cos?b = (Math.sin(?1) - Math.sin(?2)*Math.cos(?12)) / (Math.sin(?12)*Math.cos(?2));
        const ?a = Math.acos(Math.min(Math.max(cos?a, -1), 1)); // protect against rounding errors
        const ?b = Math.acos(Math.min(Math.max(cos?b, -1), 1)); // protect against rounding errors

        const ?12 = Math.sin(?2-?1)>0 ? ?a : 2*?-?a;
        const ?21 = Math.sin(?2-?1)>0 ? 2*?-?b : ?b;

        const ?1 = ?13 - ?12; // angle 2-1-3
        const ?2 = ?21 - ?23; // angle 1-2-3

        if (Math.sin(?1) == 0 && Math.sin(?2) == 0) return null; // infinite intersections
        if (Math.sin(?1) * Math.sin(?2) < 0) return null;        // ambiguous intersection (antipodal?)

        const cos?3 = -Math.cos(?1)*Math.cos(?2) + Math.sin(?1)*Math.sin(?2)*Math.cos(?12);

        const ?13 = Math.atan2(Math.sin(?12)*Math.sin(?1)*Math.sin(?2), Math.cos(?2) + Math.cos(?1)*cos?3);

        const ?3 = Math.asin(Math.sin(?1)*Math.cos(?13) + Math.cos(?1)*Math.sin(?13)*Math.cos(?13));

        const ??13 = Math.atan2(Math.sin(?13)*Math.sin(?13)*Math.cos(?1), Math.cos(?13) - Math.sin(?1)*Math.sin(?3));
        const ?3 = ?1 + ??13;

        const lat = ?3.toDegrees();
        const lon = ?3.toDegrees();

        return new LatLonSpherical(lat, lon);
    }


    /**
     * Returns (signed) distance from ‘this’ point to great circle defined by start-point and
     * end-point.
     *
     * @param   {LatLon} pathStart - Start point of great circle path.
     * @param   {LatLon} pathEnd - End point of great circle path.
     * @param   {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
     * @returns {number} Distance to great circle (-ve if to left, +ve if to right of path).
     *
     * @example
     *   const pCurrent = new LatLon(53.2611, -0.7972);
     *   const p1 = new LatLon(53.3206, -1.7297);
     *   const p2 = new LatLon(53.1887, 0.1334);
     *   const d = pCurrent.crossTrackDistanceTo(p1, p2);  // -307.5 m
     */
    crossTrackDistanceTo(pathStart, pathEnd, radius=6371e3) {
        if (!(pathStart instanceof LatLonSpherical)) pathStart = LatLonSpherical.parse(pathStart); // allow literal forms
        if (!(pathEnd instanceof LatLonSpherical)) pathEnd = LatLonSpherical.parse(pathEnd);       // allow literal forms
        const R = radius;

        const ?13 = pathStart.distanceTo(this, R) / R;
        const ?13 = pathStart.initialBearingTo(this).toRadians();
        const ?12 = pathStart.initialBearingTo(pathEnd).toRadians();

        const ?xt = Math.asin(Math.sin(?13) * Math.sin(?13 - ?12));

        return ?xt * R;
    }


    /**
     * Returns how far ‘this’ point is along a path from from start-point, heading towards end-point.
     * That is, if a perpendicular is drawn from ‘this’ point to the (great circle) path, the
     * along-track distance is the distance from the start point to where the perpendicular crosses
     * the path.
     *
     * @param   {LatLon} pathStart - Start point of great circle path.
     * @param   {LatLon} pathEnd - End point of great circle path.
     * @param   {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
     * @returns {number} Distance along great circle to point nearest ‘this’ point.
     *
     * @example
     *   const pCurrent = new LatLon(53.2611, -0.7972);
     *   const p1 = new LatLon(53.3206, -1.7297);
     *   const p2 = new LatLon(53.1887,  0.1334);
     *   const d = pCurrent.alongTrackDistanceTo(p1, p2);  // 62.331 km
     */
    alongTrackDistanceTo(pathStart, pathEnd, radius=6371e3) {
        if (!(pathStart instanceof LatLonSpherical)) pathStart = LatLonSpherical.parse(pathStart); // allow literal forms
        if (!(pathEnd instanceof LatLonSpherical)) pathEnd = LatLonSpherical.parse(pathEnd);       // allow literal forms
        const R = radius;

        const ?13 = pathStart.distanceTo(this, R) / R;
        const ?13 = pathStart.initialBearingTo(this).toRadians();
        const ?12 = pathStart.initialBearingTo(pathEnd).toRadians();

        const ?xt = Math.asin(Math.sin(?13) * Math.sin(?13-?12));

        const ?at = Math.acos(Math.cos(?13) / Math.abs(Math.cos(?xt)));

        return ?at*Math.sign(Math.cos(?12-?13)) * R;
    }


    /**
     * Returns maximum latitude reached when travelling on a great circle on given bearing from
     * ‘this’ point (‘Clairaut’s formula’). Negate the result for the minimum latitude (in the
     * southern hemisphere).
     *
     * The maximum latitude is independent of longitude; it will be the same for all points on a
     * given latitude.
     *
     * @param   {number} bearing - Initial bearing.
     * @returns {number} Maximum latitude reached.
     */
    maxLatitude(bearing) {
        const ? = Number(bearing).toRadians();

        const ? = this.lat.toRadians();

        const ?Max = Math.acos(Math.abs(Math.sin(?) * Math.cos(?)));

        return ?Max.toDegrees();
    }


    /**
     * Returns the pair of meridians at which a great circle defined by two points crosses the given
     * latitude. If the great circle doesn't reach the given latitude, null is returned.
     *
     * @param   {LatLon}      point1 - First point defining great circle.
     * @param   {LatLon}      point2 - Second point defining great circle.
     * @param   {number}      latitude - Latitude crossings are to be determined for.
     * @returns {Object|null} Object containing { lon1, lon2 } or null if given latitude not reached.
     */
    static crossingParallels(point1, point2, latitude) {
        if (point1.equals(point2)) return null; // coincident points

        const ? = Number(latitude).toRadians();

        const ?1 = point1.lat.toRadians();
        const ?1 = point1.lon.toRadians();
        const ?2 = point2.lat.toRadians();
        const ?2 = point2.lon.toRadians();

        const ?? = ?2 - ?1;

        const x = Math.sin(?1) * Math.cos(?2) * Math.cos(?) * Math.sin(??);
        const y = Math.sin(?1) * Math.cos(?2) * Math.cos(?) * Math.cos(??) - Math.cos(?1) * Math.sin(?2) * Math.cos(?);
        const z = Math.cos(?1) * Math.cos(?2) * Math.sin(?) * Math.sin(??);

        if (z * z > x * x + y * y) return null; // great circle doesn't reach latitude

        const ?m = Math.atan2(-y, x);               // longitude at max latitude
        const ??i = Math.acos(z / Math.sqrt(x*x + y*y)); // ?? from ?m to intersection points

        const ?i1 = ?1 + ?m - ??i;
        const ?i2 = ?1 + ?m + ??i;

        const lon1 = ?i1.toDegrees();
        const lon2 = ?i2.toDegrees();

        return {
            lon1: Dms.wrap180(lon1),
            lon2: Dms.wrap180(lon2),
        };
    }


    /* Rhumb - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */


    /**
     * Returns the distance travelling from ‘this’ point to destination point along a rhumb line.
     *
     * @param   {LatLon} point - Latitude/longitude of destination point.
     * @param   {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
     * @returns {number} Distance in km between this point and destination point (same units as radius).
     *
     * @example
     *   const p1 = new LatLon(51.127, 1.338);
     *   const p2 = new LatLon(50.964, 1.853);
     *   const d = p1.distanceTo(p2); //  40.31 km
     */
    rhumbDistanceTo(point, radius=6371e3) {
        if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms

        // see www.edwilliams.org/avform.htm#Rhumb

        const R = radius;
        const ?1 = this.lat.toRadians();
        const ?2 = point.lat.toRadians();
        const ?? = ?2 - ?1;
        let ?? = Math.abs(point.lon - this.lon).toRadians();
        // if dLon over 180° take shorter rhumb line across the anti-meridian:
        if (Math.abs(??) > ?) ?? = ?? > 0 ? -(2 * ? - ??) : (2 * ? + ??);

        // on Mercator projection, longitude distances shrink by latitude; q is the 'stretch factor'
        // q becomes ill-conditioned along E-W line (0/0); use empirical tolerance to avoid it
        const ?? = Math.log(Math.tan(?2 / 2 + ? / 4) / Math.tan(?1 / 2 + ? / 4));
        const q = Math.abs(??) > 10e-12 ? ?? / ?? : Math.cos(?1);

        // distance is pythagoras on 'stretched' Mercator projection, ?(??² + q²·??²)
        const ? = Math.sqrt(??*?? + q*q * ??*??); // angular distance in radians
        const d = ? * R;

        return d;
    }


    /**
     * Returns the bearing from ‘this’ point to destination point along a rhumb line.
     *
     * @param   {LatLon}    point - Latitude/longitude of destination point.
     * @returns {number}    Bearing in degrees from north.
     *
     * @example
     *   const p1 = new LatLon(51.127, 1.338);
     *   const p2 = new LatLon(50.964, 1.853);
     *   const d = p1.rhumbBearingTo(p2); // 116.7°
     */
    rhumbBearingTo(point) {
        if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms
        if (this.equals(point)) return NaN; // coincident points

        const ?1 = this.lat.toRadians();
        const ?2 = point.lat.toRadians();
        let ?? = (point.lon - this.lon).toRadians();
        // if dLon over 180° take shorter rhumb line across the anti-meridian:
        if (Math.abs(??) > ?) ?? = ?? > 0 ? -(2 * ? - ??) : (2 * ? + ??);

        const ?? = Math.log(Math.tan(?2 / 2 + ? / 4) / Math.tan(?1 / 2 + ? / 4));

        const ? = Math.atan2(??, ??);

        const bearing = ?.toDegrees();

        return Dms.wrap360(bearing);
    }


    /**
     * Returns the destination point having travelled along a rhumb line from ‘this’ point the given
     * distance on the given bearing.
     *
     * @param   {number} distance - Distance travelled, in same units as earth radius (default: metres).
     * @param   {number} bearing - Bearing in degrees from north.
     * @param   {number} [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
     * @returns {LatLon} Destination point.
     *
     * @example
     *   const p1 = new LatLon(51.127, 1.338);
     *   const p2 = p1.rhumbDestinationPoint(40300, 116.7); // 50.9642°N, 001.8530°E
     */
    rhumbDestinationPoint(distance, bearing, radius=6371e3) {
        const ?1 = this.lat.toRadians(), ?1 = this.lon.toRadians();
        const ? = Number(bearing).toRadians();

        const ? = distance / radius; // angular distance in radians

        const ?? = ? * Math.cos(?);
        let ?2 = ?1 + ??;

        // check for some daft bugger going past the pole, normalise latitude if so
        if (Math.abs(?2) > ? / 2) ?2 = ?2 > 0 ? ? - ?2 : -? - ?2;

        const ?? = Math.log(Math.tan(?2 / 2 + ? / 4) / Math.tan(?1 / 2 + ? / 4));
        const q = Math.abs(??) > 10e-12 ? ?? / ?? : Math.cos(?1); // E-W course becomes ill-conditioned with 0/0

        const ?? = ? * Math.sin(?) / q;
        const ?2 = ?1 + ??;

        const lat = ?2.toDegrees();
        const lon = ?2.toDegrees();

        return new LatLonSpherical(lat, lon);
    }


    /**
     * Returns the loxodromic midpoint (along a rhumb line) between ‘this’ point and second point.
     *
     * @param   {LatLon} point - Latitude/longitude of second point.
     * @returns {LatLon} Midpoint between this point and second point.
     *
     * @example
     *   const p1 = new LatLon(51.127, 1.338);
     *   const p2 = new LatLon(50.964, 1.853);
     *   const pMid = p1.rhumbMidpointTo(p2); // 51.0455°N, 001.5957°E
     */
    rhumbMidpointTo(point) {
        if (!(point instanceof LatLonSpherical)) point = LatLonSpherical.parse(point); // allow literal forms

        // see mathforum.org/kb/message.jspa?messageID=148837

        const ?1 = this.lat.toRadians(); let ?1 = this.lon.toRadians();
        const ?2 = point.lat.toRadians(), ?2 = point.lon.toRadians();

        if (Math.abs(?2 - ?1) > ?) ?1 += 2 * ?; // crossing anti-meridian

        const ?3 = (?1 + ?2) / 2;
        const f1 = Math.tan(? / 4 + ?1 / 2);
        const f2 = Math.tan(? / 4 + ?2 / 2);
        const f3 = Math.tan(? / 4 + ?3 / 2);
        let ?3 = ((?2 - ?1) * Math.log(f3) + ?1 * Math.log(f2) - ?2 * Math.log(f1)) / Math.log(f2 / f1);

        if (!isFinite(?3)) ?3 = (?1 + ?2) / 2; // parallel of latitude

        const lat = ?3.toDegrees();
        const lon = ?3.toDegrees();

        return new LatLonSpherical(lat, lon);
    }


    /* Area - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


    /**
     * Calculates the area of a spherical polygon where the sides of the polygon are great circle
     * arcs joining the vertices.
     *
     * @param   {LatLon[]} polygon - Array of points defining vertices of the polygon.
     * @param   {number}   [radius=6371e3] - (Mean) radius of earth (defaults to radius in metres).
     * @returns {number}   The area of the polygon in the same units as radius.
     *
     * @example
     *   const polygon = [new LatLon(0,0), new LatLon(1,0), new LatLon(0,1)];
     *   const area = LatLon.areaOf(polygon); // 6.18e9 m²
     */
    static areaOf(polygon, radius=6371e3) {
        // uses method due to Karney: osgeo-org.1560.x6.nabble.com/Area-of-a-spherical-polygon-td3841625.html;
        // for each edge of the polygon, tan(E/2) = tan(??/2)·(tan(??/2)+tan(??/2)) / (1+tan(??/2)·tan(??/2))
        // where E is the spherical excess of the trapezium obtained by extending the edge to the equator
        // (Karney's method is probably more efficient than the more widely known L’Huilier’s Theorem)

        const R = radius;

        // close polygon so that last point equals first point
        const closed = polygon[0].equals(polygon[polygon.length-1]);
        if (!closed) polygon.push(polygon[0]);

        const nVertices = polygon.length - 1;

        let S = 0; // spherical excess in steradians
        for (let v=0; v Number.EPSILON) return false;
        if (Math.abs(this.lon - point.lon) > Number.EPSILON) return false;

        return true;
    }


    /**
     * Converts ‘this’ point to a GeoJSON object.
     *
     * @returns {Object} this point as a GeoJSON ‘Point’ object.
     */
    toGeoJSON() {
        return { type: 'Point', coordinates: [ this.lon, this.lat ] };
    }


    /**
     * Returns a string representation of ‘this’ point, formatted as degrees, degrees+minutes, or
     * degrees+minutes+seconds.
     *
     * @param   {string} [format=d] - Format point as 'd', 'dm', 'dms', or 'n' for signed numeric.
     * @param   {number} [dp=4|2|0] - Number of decimal places to use: default 4 for d, 2 for dm, 0 for dms.
     * @returns {string} Comma-separated formatted latitude/longitude.
     * @throws  {RangeError} Invalid format.
     *
     * @example
     *   const greenwich = new LatLon(51.47788, -0.00147);
     *   const d = greenwich.toString();                        // 51.4779°N, 000.0015°W
     *   const dms = greenwich.toString('dms', 2);              // 51°28?40.37?N, 000°00?05.29?W
     *   const [lat, lon] = greenwich.toString('n').split(','); // 51.4779, -0.0015
     */
    toString(format='d', dp=undefined) {
        // note: explicitly set dp to undefined for passing through to toLat/toLon
        if (![ 'd', 'dm', 'dms', 'n' ].includes(format)) throw new RangeError(`invalid format ‘${format}’`);

        if (format == 'n') { // signed numeric degrees
            if (dp == undefined) dp = 4;
            return `${this.lat.toFixed(dp)},${this.lon.toFixed(dp)}`;
        }
        const lat = Dms.toLat(this.lat, format, dp);
        const lon = Dms.toLon(this.lon, format, dp);
        return `${lat}, ${lon}`;
    }

}


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

export { LatLonSpherical as default, Dms };
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/* Geodesy representation conversion functions                        (c) Chris Veness 2002-2019  */
/*                                                                                   MIT Licence  */
/* www.movable-type.co.uk/scripts/latlong.html                                                    */
/* www.movable-type.co.uk/scripts/js/geodesy/geodesy-library.html#dms                             */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

/* eslint no-irregular-whitespace: [2, { skipComments: true }] */


/**
 * Latitude/longitude points may be represented as decimal degrees, or subdivided into sexagesimal
 * minutes and seconds. This module provides methods for parsing and representing degrees / minutes
 * / seconds.
 *
 * @module dms
 */


/* Degree-minutes-seconds (& cardinal directions) separator character */
let dmsSeparator = '\u202f'; // U+202F = 'narrow no-break space'


/**
 * Functions for parsing and representing degrees / minutes / seconds.
 */
class Dms {

    // note Unicode Degree = U+00B0. Prime = U+2032, Double prime = U+2033

    /**
     * Separator character to be used to separate degrees, minutes, seconds, and cardinal directions.
     *
     * Default separator is U+202F ‘narrow no-break space’.
     *
     * To change this (e.g. to empty string or full space), set Dms.separator prior to invoking
     * formatting.
     *
     * @example
     *   import LatLon, { Dms } from '/js/geodesy/latlon-spherical.js';
     *   const p = new LatLon(51.2, 0.33).toString('dms');  // 51°?12??00??N, 000°?19??48??E
     *   Dms.separator = '';                                // no separator
     *   const p? = new LatLon(51.2, 0.33).toString('dms'); // 51°12?00?N, 000°19?48?E
     */
    static get separator()     { return dmsSeparator; }
    static set separator(char) { dmsSeparator = char; }


    /**
     * Parses string representing degrees/minutes/seconds into numeric degrees.
     *
     * This is very flexible on formats, allowing signed decimal degrees, or deg-min-sec optionally
     * suffixed by compass direction (NSEW); a variety of separators are accepted. Examples -3.62,
     * '3 37 12W', '3°37?12?W'.
     *
     * Thousands/decimal separators must be comma/dot; use Dms.fromLocale to convert locale-specific
     * thousands/decimal separators.
     *
     * @param   {string|number} dms - Degrees or deg/min/sec in variety of formats.
     * @returns {number}        Degrees as decimal number.
     *
     * @example
     *   const lat = Dms.parse('51° 28? 40.37? N');
     *   const lon = Dms.parse('000° 00? 05.29? W');
     *   const p1 = new LatLon(lat, lon); // 51.4779°N, 000.0015°W
     */
    static parse(dms) {
        // check for signed decimal degrees without NSEW, if so return it directly
        if (!isNaN(parseFloat(dms)) && isFinite(dms)) return Number(dms);

        // strip off any sign or compass dir'n & split out separate d/m/s
        const dmsParts = String(dms).trim().replace(/^-/, '').replace(/[NSEW]$/i, '').split(/[^0-9.,]+/);
        if (dmsParts[dmsParts.length-1]=='') dmsParts.splice(dmsParts.length-1);  // from trailing symbol

        if (dmsParts == '') return NaN;

        // and convert to decimal degrees...
        let deg = null;
        switch (dmsParts.length) {
            case 3:  // interpret 3-part result as d/m/s
                deg = dmsParts[0]/1 + dmsParts[1]/60 + dmsParts[2]/3600;
                break;
            case 2:  // interpret 2-part result as d/m
                deg = dmsParts[0]/1 + dmsParts[1]/60;
                break;
            case 1:  // just d (possibly decimal) or non-separated dddmmss
                deg = dmsParts[0];
                // check for fixed-width unseparated format eg 0033709W
                //if (/[NS]/i.test(dmsParts)) deg = '0' + deg;  // - normalise N/S to 3-digit degrees
                //if (/[0-9]{7}/.test(deg)) deg = deg.slice(0,3)/1 + deg.slice(3,5)/60 + deg.slice(5)/3600;
                break;
            default:
                return NaN;
        }
        if (/^-|[WS]$/i.test(dms.trim())) deg = -deg; // take '-', west and south as -ve

        return Number(deg);
    }


    /**
     * Converts decimal degrees to deg/min/sec format
     *  - degree, prime, double-prime symbols are added, but sign is discarded, though no compass
     *    direction is added.
     *  - degrees are zero-padded to 3 digits; for degrees latitude, use .slice(1) to remove leading
     *    zero.
     *
     * @private
     * @param   {number} deg - Degrees to be formatted as specified.
     * @param   {string} [format=d] - Return value as 'd', 'dm', 'dms' for deg, deg+min, deg+min+sec.
     * @param   {number} [dp=4|2|0] - Number of decimal places to use – default 4 for d, 2 for dm, 0 for dms.
     * @returns {string} Degrees formatted as deg/min/secs according to specified format.
     */
    static toDms(deg, format='d', dp=undefined) {
        if (isNaN(deg)) return null;  // give up here if we can't make a number from deg
        if (typeof deg == 'string' && deg.trim() == '') return null;
        if (typeof deg == 'boolean') return null;
        if (deg == Infinity) return null;
        if (deg == null) return null;

        // default values
        if (dp === undefined) {
            switch (format) {
                case 'd':   case 'deg':         dp = 4; break;
                case 'dm':  case 'deg+min':     dp = 2; break;
                case 'dms': case 'deg+min+sec': dp = 0; break;
                default:          format = 'd'; dp = 4; break; // be forgiving on invalid format
            }
        }

        deg = Math.abs(deg);  // (unsigned result ready for appending compass dir'n)

        let dms = null, d = null, m = null, s = null;
        switch (format) {
            default: // invalid format spec!
            case 'd': case 'deg':
                d = deg.toFixed(dp);                       // round/right-pad degrees
                if (d<100) d = '0' + d;                    // left-pad with leading zeros (note may include decimals)
                if (d<10) d = '0' + d;
                dms = d + '°';
                break;
            case 'dm': case 'deg+min':
                d = Math.floor(deg);                       // get component deg
                m = ((deg*60) % 60).toFixed(dp);           // get component min & round/right-pad
                if (m == 60) { m = (0).toFixed(dp); d++; } // check for rounding up
                d = ('000'+d).slice(-3);                   // left-pad with leading zeros
                if (m<10) m = '0' + m;                     // left-pad with leading zeros (note may include decimals)
                dms = d + '°'+Dms.separator + m + '?';
                break;
            case 'dms': case 'deg+min+sec':
                d = Math.floor(deg);                       // get component deg
                m = Math.floor((deg*3600)/60) % 60;        // get component min
                s = (deg*3600 % 60).toFixed(dp);           // get component sec & round/right-pad
                if (s == 60) { s = (0).toFixed(dp); m++; } // check for rounding up
                if (m == 60) { m = 0; d++; }               // check for rounding up
                d = ('000'+d).slice(-3);                   // left-pad with leading zeros
                m = ('00'+m).slice(-2);                    // left-pad with leading zeros
                if (s<10) s = '0' + s;                     // left-pad with leading zeros (note may include decimals)
                dms = d + '°'+Dms.separator + m + '?'+Dms.separator + s + '?';
                break;
        }

        return dms;
    }


    /**
     * Converts numeric degrees to deg/min/sec latitude (2-digit degrees, suffixed with N/S).
     *
     * @param   {number} deg - Degrees to be formatted as specified.
     * @param   {string} [format=d] - Return value as 'd', 'dm', 'dms' for deg, deg+min, deg+min+sec.
     * @param   {number} [dp=4|2|0] - Number of decimal places to use – default 4 for d, 2 for dm, 0 for dms.
     * @returns {string} Degrees formatted as deg/min/secs according to specified format.
     *
     * @example
     *   const lat = Dms.toLat(-3.62, 'dms'); // 3°37?12?S
     */
    static toLat(deg, format, dp) {
        const lat = Dms.toDms(Dms.wrap90(deg), format, dp);
        return lat===null ? '–' : lat.slice(1) + Dms.separator + (deg<0 ? 'S' : 'N');  // knock off initial '0' for lat!
    }


    /**
     * Convert numeric degrees to deg/min/sec longitude (3-digit degrees, suffixed with E/W).
     *
     * @param   {number} deg - Degrees to be formatted as specified.
     * @param   {string} [format=d] - Return value as 'd', 'dm', 'dms' for deg, deg+min, deg+min+sec.
     * @param   {number} [dp=4|2|0] - Number of decimal places to use – default 4 for d, 2 for dm, 0 for dms.
     * @returns {string} Degrees formatted as deg/min/secs according to specified format.
     *
     * @example
     *   const lon = Dms.toLon(-3.62, 'dms'); // 3°37?12?W
     */
    static toLon(deg, format, dp) {
        const lon = Dms.toDms(Dms.wrap180(deg), format, dp);
        return lon===null ? '–' : lon + Dms.separator + (deg<0 ? 'W' : 'E');
    }


    /**
     * Converts numeric degrees to deg/min/sec as a bearing (0°..360°).
     *
     * @param   {number} deg - Degrees to be formatted as specified.
     * @param   {string} [format=d] - Return value as 'd', 'dm', 'dms' for deg, deg+min, deg+min+sec.
     * @param   {number} [dp=4|2|0] - Number of decimal places to use – default 4 for d, 2 for dm, 0 for dms.
     * @returns {string} Degrees formatted as deg/min/secs according to specified format.
     *
     * @example
     *   const lon = Dms.toBrng(-3.62, 'dms'); // 356°22?48?
     */
    static toBrng(deg, format, dp) {
        const brng =  Dms.toDms(Dms.wrap360(deg), format, dp);
        return brng===null ? '–' : brng.replace('360', '0');  // just in case rounding took us up to 360°!
    }


    /**
     * Converts DMS string from locale thousands/decimal separators to JavaScript comma/dot separators
     * for subsequent parsing.
     *
     * Both thousands and decimal separators must be followed by a numeric character, to facilitate
     * parsing of single lat/long string (in which whitespace must be left after the comma separator).
     *
     * @param   {string} str - Degrees/minutes/seconds formatted with locale separators.
     * @returns {string} Degrees/minutes/seconds formatted with standard Javascript separators.
     *
     * @example
     *   const lat = Dms.fromLocale('51°28?40,12?N');                          // '51°28?40.12?N' in France
     *   const p = new LatLon(Dms.fromLocale('51°28?40,37?N, 000°00?05,29?W'); // '51.4779°N, 000.0015°W' in France
     */
    static fromLocale(str) {
        const locale = (123456.789).toLocaleString();
        const separator = { thousands: locale.slice(3, 4), decimal: locale.slice(7, 8) };
        return str.replace(separator.thousands, '?').replace(separator.decimal, '.').replace('?', ',');
    }


    /**
     * Converts DMS string from JavaScript comma/dot thousands/decimal separators to locale separators.
     *
     * Can also be used to format standard numbers such as distances.
     *
     * @param   {string} str - Degrees/minutes/seconds formatted with standard Javascript separators.
     * @returns {string} Degrees/minutes/seconds formatted with locale separators.
     *
     * @example
     *   const Dms.toLocale('123,456.789');                   // '123.456,789' in France
     *   const Dms.toLocale('51°28?40.12?N, 000°00?05.31?W'); // '51°28?40,12?N, 000°00?05,31?W' in France
     */
    static toLocale(str) {
        const locale = (123456.789).toLocaleString();
        const separator = { thousands: locale.slice(3, 4), decimal: locale.slice(7, 8) };
        return str.replace(/,([0-9])/, '?$1').replace('.', separator.decimal).replace('?', separator.thousands);
    }


    /**
     * Returns compass point (to given precision) for supplied bearing.
     *
     * @param   {number} bearing - Bearing in degrees from north.
     * @param   {number} [precision=3] - Precision (1:cardinal / 2:intercardinal / 3:secondary-intercardinal).
     * @returns {string} Compass point for supplied bearing.
     *
     * @example
     *   const point = Dms.compassPoint(24);    // point = 'NNE'
     *   const point = Dms.compassPoint(24, 1); // point = 'N'
     */
    static compassPoint(bearing, precision=3) {
        if (![ 1, 2, 3 ].includes(Number(precision))) throw new RangeError(`invalid precision ‘${precision}’`);
        // note precision could be extended to 4 for quarter-winds (eg NbNW), but I think they are little used

        bearing = Dms.wrap360(bearing); // normalise to range 0..360°

        const cardinals = [
            'N', 'NNE', 'NE', 'ENE',
            'E', 'ESE', 'SE', 'SSE',
            'S', 'SSW', 'SW', 'WSW',
            'W', 'WNW', 'NW', 'NNW' ];
        const n = 4 * 2**(precision-1); // no of compass points at req’d precision (1=>4, 2=>8, 3=>16)
        const cardinal = cardinals[Math.round(bearing*n/360)%n * 16/n];

        return cardinal;
    }


    /**
     * Constrain degrees to range 0..360 (e.g. for bearings); -1 => 359, 361 => 1.
     *
     * @private
     * @param {number} degrees
     * @returns degrees within range 0..360.
     */
    static wrap360(degrees) {
        if (0<=degrees && degrees<360) return degrees; // avoid rounding due to arithmetic ops if within range
        return (degrees%360+360) % 360; // sawtooth wave p:360, a:360
    }

    /**
     * Constrain degrees to range -180..+180 (e.g. for longitude); -181 => 179, 181 => -179.
     *
     * @private
     * @param {number} degrees
     * @returns degrees within range -180..+180.
     */
    static wrap180(degrees) {
        if (-180 -89, 91 => 89.
     *
     * @private
     * @param {number} degrees
     * @returns degrees within range -90..+90.
     */
    static wrap90(degrees) {
        if (-90<=degrees && degrees<=90) return degrees; // avoid rounding due to arithmetic ops if within range
        return Math.abs((degrees%360 + 270)%360 - 180) - 90; // triangle wave p:360 a:±90 TODO: fix e.g. -315°
    }

}


// Extend Number object with methods to convert between degrees & radians
Number.prototype.toRadians = function() { return this * Math.PI / 180; };
Number.prototype.toDegrees = function() { return this * 180 / Math.PI; };

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

export default Dms;