1#
2# This file is part of:
3#
4#  gpsman --- GPS Manager: a manager for GPS receiver data
5#
6# Copyright (c) 1998-2013 Miguel Filgueiras migfilg@t-online.de
7#
8#    This program is free software; you can redistribute it and/or modify
9#      it under the terms of the GNU General Public License as published by
10#      the Free Software Foundation; either version 3 of the License, or
11#      (at your option) any later version.
12#
13#      This program is distributed in the hope that it will be useful,
14#      but WITHOUT ANY WARRANTY; without even the implied warranty of
15#      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16#      GNU General Public License for more details.
17#
18#      You should have received a copy of the GNU General Public License
19#      along with this program.
20#
21#  File: gdata.tcl
22#  Last change:  6 October 2013
23#
24
25
26## manipulation of geo-referenced data structures
27#   consisting of a list of pairs $coord1 $col, in increasing order of $coord1
28#   where $col is a list of pairs $coord2 $data, in increasing order of $coord2
29#   and $data a list of information for each individual at $coord1 $coord2
30
31
32proc AddSeqToGData {gdata coord1 coord2 info} {
33    # add information to geo-referenced data structure
34    # sequential search
35
36    set ne [list $coord1 [list [list $coord2 [list $info]]]]
37    set i 0
38    foreach e $gdata {
39	if { [set le [lindex $e 0]] > $coord1 } {
40	    break
41	} elseif { $le == $coord1 } {
42	    set nee [list $coord2 [list $info]] ; set col [lindex $e 1]
43	    set j 0
44	    foreach ee $col {
45		if { [set lee [lindex $ee 0]] > $coord2 } {
46		    break
47		} elseif { $lee == $coord2 } {
48		    set nee [list $coord2 [linsert [lindex $ee 1] 0 $info]]
49		    set col [lreplace $col $j $j $nee]
50		    return [lreplace $gdata $i $i [list $coord1 $col]]
51		}
52		incr j
53	    }
54	    set col [linsert $col $j $nee]
55	    return [lreplace $gdata $i $i [list $coord1 $col]]
56	}
57	incr i
58    }
59    return [linsert $gdata $i $ne]
60}
61
62proc AddToGData {gdata coord1 coord2 info} {
63    # add information to geo-referenced data structure
64    # binary search on top level list
65
66    set ne [list $coord1 [list [list $coord2 [list $info]]]]
67    if { [set b [llength $gdata]] == 0 } {
68	return [list $ne]
69    }
70    set a 0
71    while 1 {
72	set i [expr ($a+$b)/2]
73	set e [lindex $gdata $i]
74	if { [set le [lindex $e 0]] > $coord1 } {
75	    if { $a == $i } {
76		break
77	    }
78	    set b $i
79	} elseif { $le == $coord1 } {
80	    set nee [list $coord2 [list $info]] ; set col [lindex $e 1]
81	    set j 0
82	    foreach ee $col {
83		if { [set lee [lindex $ee 0]] > $coord2 } {
84		    break
85		} elseif { $lee == $coord2 } {
86		    set nee [list $coord2 [linsert [lindex $ee 1] 0 $info]]
87		    set col [lreplace $col $j $j $nee]
88		    return [lreplace $gdata $i $i [list $coord1 $col]]
89		}
90		incr j
91	    }
92	    set col [linsert $col $j $nee]
93	    return [lreplace $gdata $i $i [list $coord1 $col]]
94	} elseif { $a == $i } {
95	    incr a
96	    break
97	} else { set a $i }
98    }
99    return [linsert $gdata $a $ne]
100}
101
102proc LookupGData {gdata coord1 coord2} {
103    # find information in geo-referenced data structure
104    # return "" if nothing found
105    # binary search on top level list
106
107    if { [set b [llength $gdata]] == 0 } {
108	return ""
109    }
110    set a 0
111    while 1 {
112	set i [expr ($a+$b)/2]
113	set e [lindex $gdata $i]
114	if { [set le [lindex $e 0]] > $coord1 } {
115	    if { $a == $i } { break }
116	    set b $i
117	} elseif { $le == $coord1 } {
118	    foreach ee [lindex $e 1] {
119		if { [set lee [lindex $ee 0]] > $coord2 } {
120		    return ""
121		} elseif { $lee == $coord2 } {
122		    return [lindex $ee 1]
123		}
124	    }
125	    return ""
126	} elseif { $a == $i } { break } else { set a $i }
127    }
128    return ""
129}
130
131proc LookupQuadrGData {gdata coord1 coord2 dc1 dc2} {
132    # find information in geo-referenced data structure
133    #  for each point with coordinates ($coord10,$coord20) if ($coord1,$coord2)
134    #  belongs to the quadrangle centred on it and having as corners
135    #  ($coord10-$dc1/2,$coord20-$dc2/2) ($coord10+$dc1/2,$coord20+$dc2/2)
136    #  with $dc1 and $dc2 > 0
137    # return list of information for each point
138    # binary search on top level list
139
140    if { $dc1 <= 0 || $dc2 <= 0 || [set b [llength $gdata]] == 0 } {
141	return ""
142    }
143    set hdc1 [expr 0.5*$dc1] ; set hdc2 [expr 0.5*$dc2]
144    set a 0 ; set res ""
145    while 1 {
146	set i [expr ($a+$b)/2]
147	set e [lindex $gdata $i]
148	if { [set le [lindex $e 0]]-$hdc1 > $coord1 } {
149	    if { $a == $i } { return $res }
150	    set b $i
151	} elseif { $le+$hdc1 >= $coord1 } {
152	    foreach ee [lindex $e 1] {
153		if { [set lee [lindex $ee 0]]-$hdc2 > $coord2 } {
154		    break
155		} elseif { $lee+$hdc2 >= $coord2 } {
156		    set res [concat $res [lindex $ee 1]]
157		}
158	    }
159	    break
160	} elseif { $a == $i } { break } else { set a $i }
161    }
162    for { set j [expr $i+1] } { $j != $b } { incr j } {
163	set e [lindex $gdata $j]
164	if { [set le [lindex $e 0]]-$hdc1 > $coord1 } { break }
165	if {  $le+$hdc1 >= $coord1 } {
166	    foreach ee [lindex $e 1] {
167		if { [set lee [lindex $ee 0]]-$hdc2 > $coord2 } {
168		    break
169		} elseif { $lee+$hdc2 >= $coord2 } {
170		    set res [concat $res [lindex $ee 1]]
171		}
172	    }
173	}
174    }
175    incr a -1
176    for { set j [expr $i-1] } { $j != $a } { incr j -1 } {
177	set e [lindex $gdata $j]
178	if { [set le [lindex $e 0]]+$hdc1 < $coord1 } { break }
179	if {  $le-$hdc1 <= $coord1 } {
180	    foreach ee [lindex $e 1] {
181		if { [set lee [lindex $ee 0]]-$hdc2 > $coord2 } {
182		    break
183		} elseif { $lee+$hdc2 >= $coord2 } {
184		    set res [concat $res [lindex $ee 1]]
185		}
186	    }
187	}
188    }
189    return $res
190}
191
192#### geo-referencing time-stamped information
193
194proc GeoAdapt {type ix datum lts} {
195    # geo-reference time-stamped information in list $lts adapting it to
196    #  item of given $type and index $ix
197    # if $type==TR, interpolate/extrapolate positions from time-stamps of TPs
198    # if $type==GR, take each WP in sequence from GR after sorting time-stamps
199    #  $lts is a list of tuples whose first element is the time-stamp in
200    #       seconds from $YEAR0, sorted by time-stamps
201    #  $datum is the datum for the computed positions
202    # return "" on error, or list with tuples obtained from the tuples in $lst
203    #  by inserting the latitude and longitude (decimal degrees, for $datum)
204    #  and altitude (metre) at the beginning
205    global TRTPoints TRDatum DataIndex GRConts WPPosn WPDatum WPAlt
206
207    set glts ""
208    if { $type == "TR" } {
209	foreach v "la lo s alt" n "TPlatd TPlongd TPsecs TPalt" {
210	    set ${v}ix $DataIndex($n)
211	    set prev$v undef
212	}
213	if { $TRDatum($ix) != $datum } {
214	    set tps [ChangeTPsDatum $TRTPoints($ix) $TRDatum($ix) $datum]
215	} else { set tps $TRTPoints($ix) }
216	set tp [lindex $tps 0] ; set tps [lreplace $tps 0 0]
217	foreach x "la lo s alt" {
218	    set tpt$x [lindex $tp [set ${x}ix]]
219	}
220	set tptalt [lindex $tptalt 0]
221	foreach tuple $lts {
222	    set secs [lindex $tuple 0]
223	    while { $secs > $tpts && $tps != "" } {
224		set tp [lindex $tps 0] ; set tps [lreplace $tps 0 0]
225		foreach x "la lo s alt" {
226		    set prev$x [set tpt$x]
227		    set tpt$x [lindex $tp [set ${x}ix]]
228		}
229		set tptalt [lindex $tptalt 0]
230	    }
231	    # default to postion of current TP if in trouble
232	    #  this may imply several tuples to have the same position
233	    set la $tptla ; set lo $tptlo ; set alt $tptalt
234	    if { $secs != $tpts } {
235		if { $prevla != "undef" } {
236		    # liner interpolation or extrapolation
237		    if { ! [catch \
238			   {set ratio [expr 1.0*($secs-$prevs)/($tpts-$prevs)]}]
239		     } {
240			set la [expr $prevla+$ratio*($tptla-$prevla)]
241			set lo [expr $prevlo+$ratio*($tptlo-$prevlo)]
242			if { $tptalt == "" || $prevalt == "" } {
243			    set alt ""
244			} else {
245			    set alt [expr $prevalt+$ratio*($tptalt-$prevalt)]
246			}
247		    }
248		} elseif { $secs < $tpts && $tps != "" } {
249		    set tp2 [lindex $tps 0]
250		    foreach x "la lo s alt" {
251			set tp2$x [lindex $tp2 [set ${x}ix]]
252		    }
253		    set tp2alt [lindex $tp2alt 0]
254		    if { ! [catch \
255			    {set ratio [expr 1.0*($secs-$tp2s)/($tpts-$tp2s)]}]
256		     } {
257			set la [expr $tp2la+$ratio*($tptla-$tp2la)]
258			set lo [expr $tp2lo+$ratio*($tptlo-$tp2lo)]
259			if { $tptalt == "" || $tp2alt == "" } {
260			    set alt ""
261			} else {
262			    set alt [expr $tp2alt+$ratio*($tptalt-$tp2alt)]
263			}
264		    }
265		}
266	    }
267	    lappend glts [linsert $tuple 0 $la $lo $alt]
268	}
269    } else {
270	# $type == "GR"
271	set wpns ""
272	foreach p $GRConts($ix) {
273	    if { [lindex $p 0] == "WP" } {
274		set wpns [lindex $p 1]
275		break
276	    }
277	}
278	if { $wpns == "" } {
279	    GMMessage $MESS(allundef)
280	    return ""
281	}
282	foreach tuple $lts {
283	    set secs [lindex $tuple 0]
284	    while { $wpns != "" && \
285			[set wpix [IndexNamed WP [lindex $wpns 0]]] == -1 } {
286		GMMessage [format $MESS(undefinedWP) [lindex $wpns 0]]
287		set wpns [lreplace $wpns 0 0]
288	    }
289	    if { $wpns == "" } {
290		# $la, ...  are defined because $wpns was initially non-empty
291		# not enough WPs: use last one for the remaining tuples
292		lappend glts [linsert $tuple 0 $la $lo $alt]
293		continue
294	    }
295	    set la [lindex $WPPosn($wpix) 0] ; set lo [lindex $WPPosn($wpix) 1]
296	    if { $WPDatum($wpix) != $datum } {
297		set p [ToDatum $la $lo $WPDatum($wpix) $datum]
298		set la [lindex $p 0] ; set lo [lindex $p 1]
299	    }
300	    set alt [lindex $WPAlt($wpix) 0]
301	    lappend glts [linsert $tuple 0 $la $lo $alt]
302	    set wpns [lreplace $wpns 0 0]
303	}
304    }
305    return $glts
306}
307
308