comparison weather.tcl @ 436:3c816fdc302f

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 8efbb045d44d
children e28c3347a948
comparison
equal deleted inserted replaced
435:5edaebbbb7f2 436:3c816fdc302f
325 weather_usage $upublic $unick $uchan $weather_msg_usage_nearest 325 weather_usage $upublic $unick $uchan $weather_msg_usage_nearest
326 return 0 326 return 0
327 } 327 }
328 328
329 # Check argument types 329 # Check argument types
330 set d_lat [lindex $nlist 0] 330 set pt1_lat [lindex $nlist 0]
331 set d_lng [lindex $nlist 1] 331 set pt1_lng [lindex $nlist 1]
332 if {![string is double -strict $d_lat] || ![string is double -strict $d_lng]} { 332 if {![string is double -strict $pt1_lat] || ![string is double -strict $pt1_lng]} {
333 weather_msg $upublic $unick $uchan $weather_msg_usage_nearest_invalid 333 weather_msg $upublic $unick $uchan $weather_msg_usage_nearest_invalid
334 return 0 334 return 0
335 } 335 }
336 336
337 # Calculate distances between given coordinates for each location 337 # Calculate distances between given coordinates for each location
338 set result {} 338 set result {}
339 foreach {ukey uvalue} [array get weather_data] { 339 foreach {ukey uvalue} [array get weather_data] {
340 if {![string match "w_*" $ukey]} { 340 if {![string match "w_*" $ukey]} {
341 set delta_lat [expr {$d_lat - [lindex $uvalue 2]}] 341 set pt2_lat [lindex $uvalue 2]
342 set delta_lng [expr {$d_lng - [lindex $uvalue 3]}] 342 set pt2_lng [lindex $uvalue 3]
343 set dist [expr { sqrt($delta_lat * $delta_lat + $delta_lng * $delta_lng) }] 343
344 # Calculate "great circle angle" of two lat/long coordinate pairs
345 # https://en.wikipedia.org/wiki/Great-circle_distance
346 #
347 # The result then can be used to calculate approximate distance
348 # between the points when radius of the sphere "R" is known.
349 set d_lng [expr { $pt2_lng - $pt1_lng }]
350 set c_lat1 [expr { cos($pt1_lat) }]
351 set c_lat2 [expr { cos($pt2_lat) }]
352 set s_lat1 [expr { sin($pt1_lat) }]
353 set s_lat2 [expr { sin($pt2_lat) }]
354
355 set f1 [expr { $c_lat2 * sin($d_lng) }]
356 set f2 [expr { $c_lat2 * cos($d_lng) }]
357 set f3 [expr { $c_lat1 * $s_lat2 - $s_lat1 * $f2 }]
358 set c_ang [expr { atan2(sqrt($f1 * $f1 + $f3 * $f3), $s_lat1 * $s_lat2 + $c_lat1 * $f2) }]
359
360 # Calculate distance based on Earth radius approximation in km * 1000m
361 set dist [expr { 6371.439 * 1000.0 * $c_ang }]
344 lappend result [list $ukey $dist] 362 lappend result [list $ukey $dist]
345 } 363 }
346 } 364 }
347 365
348 # Sort the list by distance 366 # Sort the list by distance
354 lappend uresult [weather_get_str $weather_data([lindex $uval 0]) $weather_msg_list_nearest] 372 lappend uresult [weather_get_str $weather_data([lindex $uval 0]) $weather_msg_list_nearest]
355 } 373 }
356 374
357 # Print out the result 375 # Print out the result
358 set res [join $uresult " ; "] 376 set res [join $uresult " ; "]
359 weather_msg $upublic $unick $uchan $weather_msg_nearest_stations [list $d_lat $d_lng $res] 377 weather_msg $upublic $unick $uchan $weather_msg_nearest_stations [list $pt1_lat $pt1_lng $res]
360 return 0 378 return 0
361 } elseif {$rarg == "vakio" || $rarg == "default" || $rarg == "vakiot" || $rarg == "defaults"} { 379 } elseif {$rarg == "vakio" || $rarg == "default" || $rarg == "vakiot" || $rarg == "defaults"} {
362 # List or set the default weather station name patterns for this user 380 # List or set the default weather station name patterns for this user
363 381
364 # Access check 382 # Access check