Mercurial > hg > egg-tcls

changeset 436:3c816fdc302f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
weather: Implement the "great circle distance" algorithm for calculating distances between two latitude/longitude coordinate pairs, which should be more accurate than simply using Euclidean distance.
author Matti Hamalainen <ccr@tnsp.org>
date Mon, 09 Jan 2017 05:45:40 +0200
parents 5edaebbbb7f2
children e28c3347a948
files weather.tcl
diffstat 1 files changed, 25 insertions(+), 7 deletions(-) [+]
line wrap: on
line diff
--- a/weather.tcl	Mon Jan 09 04:19:10 2017 +0200
+++ b/weather.tcl	Mon Jan 09 05:45:40 2017 +0200
@@ -327,9 +327,9 @@
     }
 
     # Check argument types
-    set d_lat [lindex $nlist 0]
-    set d_lng [lindex $nlist 1]
-    if {![string is double -strict $d_lat] || ![string is double -strict $d_lng]} {
+    set pt1_lat [lindex $nlist 0]
+    set pt1_lng [lindex $nlist 1]
+    if {![string is double -strict $pt1_lat] || ![string is double -strict $pt1_lng]} {
       weather_msg $upublic $unick $uchan $weather_msg_usage_nearest_invalid
       return 0
     }
@@ -338,9 +338,27 @@
     set result {}
     foreach {ukey uvalue} [array get weather_data] {
       if {![string match "w_*" $ukey]} {
-        set delta_lat [expr {$d_lat - [lindex $uvalue 2]}]
-        set delta_lng [expr {$d_lng - [lindex $uvalue 3]}]
-        set dist [expr { sqrt($delta_lat * $delta_lat + $delta_lng * $delta_lng) }]
+        set pt2_lat [lindex $uvalue 2]
+        set pt2_lng [lindex $uvalue 3]
+
+        # Calculate "great circle angle" of two lat/long coordinate pairs
+        # https://en.wikipedia.org/wiki/Great-circle_distance
+        #
+        # The result then can be used to calculate approximate distance
+        # between the points when radius of the sphere "R" is known.
+        set d_lng [expr { $pt2_lng - $pt1_lng }]
+        set c_lat1 [expr { cos($pt1_lat) }]
+        set c_lat2 [expr { cos($pt2_lat) }]
+        set s_lat1 [expr { sin($pt1_lat) }]
+        set s_lat2 [expr { sin($pt2_lat) }]
+
+        set f1 [expr { $c_lat2 * sin($d_lng) }]
+        set f2 [expr { $c_lat2 * cos($d_lng) }]
+        set f3 [expr { $c_lat1 * $s_lat2 - $s_lat1 * $f2 }]
+        set c_ang [expr { atan2(sqrt($f1 * $f1 + $f3 * $f3), $s_lat1 * $s_lat2 + $c_lat1 * $f2) }]
+
+        # Calculate distance based on Earth radius approximation in km * 1000m
+        set dist [expr { 6371.439 * 1000.0 * $c_ang }]
         lappend result [list $ukey $dist]
       }
     }
@@ -356,7 +374,7 @@
 
     # Print out the result
     set res [join $uresult " ; "]
-    weather_msg $upublic $unick $uchan $weather_msg_nearest_stations [list $d_lat $d_lng $res]
+    weather_msg $upublic $unick $uchan $weather_msg_nearest_stations [list $pt1_lat $pt1_lng $res]
     return 0
   } elseif {$rarg == "vakio" || $rarg == "default" || $rarg == "vakiot" || $rarg == "defaults"} {
     # List or set the default weather station name patterns for this user