source: ntrip/trunk/BNC/scripts/cmbPlot.pl@ 10270

Last change on this file since 10270 was 10119, checked in by stuerze, 17 months ago

minor changes

  • Property svn:executable set to *
File size: 17.8 KB
Line 
1#!/usr/bin/env perl
2
3# ========================================================================
4# cmbPlot.pl
5# ========================================================================
6#
7# Purpose: Script using gnuplot to plot BNC's combination results
8#
9# Comment: stored values:
10# in case of "single epoch" and "filter" combination
11# - residuals per satellite and AC
12# - resultant combined clock correction in [seconds]
13# - full clock in [m]
14# additional for "filter" combination
15# - estimated clock offset common for all satellites of a system and one AC
16# - estimated clock offset per satellite and AC
17#
18# Author: Andrea Stuerze
19# ========================================================================
20
21# Uses
22use strict;
23use warnings;
24use File::Spec;
25use Getopt::Long;
26use Chart::Gnuplot;
27use PDL;
28use Data::Dumper qw(Dumper);
29use diagnostics;
30use File::Basename;
31use Date::Manip;
32
33
34my ($home) = glob "~";
35# -----------------------------------------------------------------------------
36# required to create a pdf
37# -----------------------------------------------------------------------------
38use PDF::API2;
39use constant {
40 mm => 25.4 / 72,
41 inch => 1 / 72,
42 pt => 1,
43}; # There are 72 postscript points in an inch and there are 25.4 millimeters in an inch.
44
45# -----------------------------------------------------------------------------
46# Options
47# -----------------------------------------------------------------------------
48# init options
49my $help = 0;
50my $plotType = "";
51my @logFiles = ();
52
53# read options from command line
54GetOptions(
55 'help' => \$help,
56 'plotType=s' => \$plotType,
57 'logFiles=s' => \@logFiles,
58);
59HELP_MESSAGE() if $help;
60@logFiles = split ( /[, ]+/, join ( ',', @logFiles ) );
61unless (@logFiles) {
62 print "ERROR: logfiles missing\n";
63 exit;
64}
65print " plotTpe: $plotType\n logfiles: @logFiles\n\n";
66
67# -----------------------------------------------------------------------------
68# generate data sets for gnuplot
69# -----------------------------------------------------------------------------
70#my $old_epochSec = 0;
71my $epochSec = 0;
72my $epochDiff = 0;
73
74# for pdf generation
75my ( $png, $page, $headline, $headline_text );
76my $y0 = 180 / mm;
77my ( $x, $y, $width, $height ) = ( 40 / mm, $y0, 130 / mm, 80 / mm );
78my $dy = $height + 10 / mm;
79
80# -----------------------------------------------------------------------------
81# create pdf for plot results
82# -----------------------------------------------------------------------------
83my ( $inputFilename, $inputDir ) = fileparse( $logFiles[0] ); # , '\..*'
84my $pdf_name = sprintf ( "%s.pdf", $inputFilename );
85my $pdf = PDF::API2->new( -file => "$inputDir$pdf_name" );
86my $font1 = $pdf->corefont('Helvetica-Bold');
87
88# -----------------------------------------------------------------------------
89# read logfile(s) and save data
90# -----------------------------------------------------------------------------
91#chdir ( "$home/tmp" ) || die "Could not change to dir $home/tmp\n";
92
93my ( %CLK_CORR, %CLK_FULL, %RES, %OFFSET_AC, %OFFSET_SAT ); # %CLK
94foreach (@logFiles) {
95 #my $file = File::Spec->catfile( $inputDir, $_ );
96 my $file = $_;
97 print "Parse logfile $file ...\n";
98 my ( @array, $sys, $sat, $epo, $val, $ac );
99 open ( INPUT, "<", $file ) || die "Could not open file '$file': $!";
100 while (<INPUT>) {
101 if ( $_ =~ /Clk Corr / ) { # 2015-04-06 23:59:30.000 Clk Corr G06 0.2371 +- 0.0291
102 @array = split ( /\s+/, $_ );
103 $epo = $array[0] . "T" . $array[1];
104 $sat = $array[4];
105 $sys = substr ( $sat, 0, 1 );
106 $val = $array[5];
107 push ( @{ $CLK_CORR{$sys}{$sat}{EPOCH} }, $epo );
108 push ( @{ $CLK_CORR{$sys}{$sat}{DATA} }, $val );
109 }
110 elsif ( $_ =~ /Full Clock / ) { # 2015-04-06 23:59:30.000 Full Clock G04 81 -3205.8546
111 @array = split ( /\s+/, $_ );
112 $epo = $array[0] . "T" . $array[1];
113 $sat = $array[4];
114 $sys = substr ( $sat, 0, 1 );
115 $val = $array[6];
116 push ( @{ $CLK_FULL{$sys}{$sat}{EPOCH} }, $epo );
117 push ( @{ $CLK_FULL{$sys}{$sat}{DATA} }, $val );
118 }
119 elsif ( $_ =~ /res = / ) { # 2015-04-05 23:59:40.000 NRCAN G01 res = -0.0003
120 @array = split ( /\s+/, $_ );
121 $epo = $array[0] . "T" . $array[1];
122 $ac = $array[2];
123 $sat = $array[3];
124 $sys = substr ( $sat, 0, 1 );
125 $val = $array[6];
126 push ( @{ $RES{$ac}{$sys}{$sat}{EPOCH} }, $epo );
127 push ( @{ $RES{$ac}{$sys}{$sat}{DATA} }, $val );
128 }
129 elsif ( $_ =~ /AC Offset / ) { # AC Offset DLR C 1.9002 +- 35.3112
130 @array = split ( /\s+/, $_ );
131 $epo = $array[0] . "T" . $array[1];
132 $sys = $array[5];
133 $ac = $array[4];
134 $val = $array[6];
135 push ( @{ $OFFSET_AC{$sys}{$ac}{EPOCH} }, $epo );
136 push ( @{ $OFFSET_AC{$sys}{$ac}{DATA} }, $val );
137 }
138 elsif ( $_ =~ /Sat Offset / ) { # 2015-04-05 23:59:40.000 Sat Offset BKG G01 0.5483 +- 100.0000
139 @array = split ( /\s+/, $_ );
140 $epo = $array[0] . "T" . $array[1];
141 $ac = $array[4];
142 $sat = $array[5];
143 $sys = substr ( $sat, 0, 1 );
144 $val = $array[6];
145 push ( @{ $OFFSET_SAT{$ac}{$sys}{$sat}{EPOCH} }, $epo );
146 push ( @{ $OFFSET_SAT{$ac}{$sys}{$sat}{DATA} }, $val );
147 }
148 }
149 close INPUT;
150} # ----- end foreach logfile -----
151
152# -----------------------------------------------------------------------------
153# plot several data sets
154# -----------------------------------------------------------------------------
155######### CLOCK CORR #####################
156print "Plot CLOCK CORR ...\n";
157$page = $pdf->page();
158$page->mediabox('A4');
159$headline = sprintf ( "Clock Corrections (%s)", $inputFilename );
160$headline_text = $page->text;
161$headline_text->font( $font1, 11 / pt );
162$headline_text->translate( 15 / mm, 280 / mm );
163$headline_text->text($headline);
164$y = $y0 + $dy;
165my $yrange = 20.0;
166if ( $plotType eq "KF" ) {
167 $yrange = 10.0;
168}
169
170# SYSTEM
171foreach my $key_sys ( sort keys %CLK_CORR ) {
172 print " Plot CLOCK CORR for $key_sys\n";
173 my @datasets; # init datasets
174 my $pngName = sprintf ( "%s_CLK_CORR_%s.png", $inputFilename, $key_sys );
175 my $chart = Chart::Gnuplot->new(
176 output => $pngName,
177 terminal => 'png',
178 title => "$key_sys", # G,R
179 ylabel => "Clock corrections [s]",
180 yrange => [ -$yrange, $yrange ],
181 xlabel => "Time [h]",
182 timeaxis => 'x',
183 xtics => { labelfmt => '%H:%M', rotate => '-270' },
184 legend => { position => "outside below", },
185 grid => 'on',
186 );
187
188 # SATELLITE
189 foreach my $key_sat ( sort keys %{ $CLK_CORR{$key_sys} } ) {
190 my $colour = sprintf ( "%d", ( substr $key_sat, -2 ) + 1 );
191 my $dataset = Chart::Gnuplot::DataSet->new(
192 xdata => $CLK_CORR{$key_sys}{$key_sat}{EPOCH},
193 ydata => $CLK_CORR{$key_sys}{$key_sat}{DATA},
194 title => "$key_sat",
195 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
196 linetype => $colour,
197 style => " linespoints ",
198 );
199 push ( @datasets, $dataset );
200 }
201 $chart->plot2d(@datasets);
202 $y = $y - $dy;
203 if ( $y < 30 / mm ) {
204 $page = $pdf->page();
205 $page->mediabox('A4');
206 $y = $y0;
207 }
208 $png = $page->gfx();
209 die ("Unable to find image file: $!") unless -e $pngName;
210 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
211}
212
213######### CLOCK FULL #####################
214print "Plot FULL CLOCK ...\n";
215$page = $pdf->page();
216$page->mediabox('A4');
217$headline = sprintf ( "Full Clock (%s)", $inputFilename );
218$headline_text = $page->text;
219$headline_text->font( $font1, 11 / pt );
220$headline_text->translate( 15 / mm, 280 / mm );
221$headline_text->text($headline);
222$y = $y0 + $dy;
223
224# SYSTEM
225foreach my $key_sys ( sort keys %CLK_FULL ) {
226 my @datasets; # init datasets
227 my $pngName = sprintf ( "%s_CLK_FULL_%s.png", $inputFilename, $key_sys );
228 my $chart = Chart::Gnuplot->new(
229 output => $pngName,
230 terminal => 'png',
231 title => "$key_sys",
232 ylabel => "Full Clock [m]", #yrange => [" -3 ", " 3 "],
233 xlabel => "Time [h]",
234 timeaxis => 'x',
235 xtics => { labelfmt => '%H:%M', rotate => '-270', },
236 legend => { position => "outside below", },
237 grid => 'on',
238 );
239
240 # SATELLITE
241 foreach my $key_sat ( sort keys %{ $CLK_FULL{$key_sys} } ) {
242 my $colour = sprintf ( "%d", ( substr $key_sat, -2 ) + 1 );
243 my $dataset = Chart::Gnuplot::DataSet->new(
244 xdata => $CLK_FULL{$key_sys}{$key_sat}{EPOCH},
245 ydata => $CLK_FULL{$key_sys}{$key_sat}{DATA},
246 title => "$key_sat",
247 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
248 linetype => $colour,
249 style => " linespoints ",
250 );
251 push ( @datasets, $dataset );
252 }
253 $chart->plot2d(@datasets);
254 $y = $y - $dy;
255 if ( $y < 30 / mm ) {
256 $page = $pdf->page();
257 $page->mediabox('A4');
258 $y = $y0;
259 }
260 $png = $page->gfx();
261 die ("Unable to find image file: $!") unless -e $pngName;
262 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
263}
264
265######### OFFSET SAT #####################
266if ( $plotType eq "KF" ) {
267 print "Plot OFFSET SAT ...\n";
268 $page = $pdf->page();
269 $page->mediabox('A4');
270 $headline = sprintf ( "AC and satellite specific offset (%s)", $inputFilename );
271 $headline_text = $page->text;
272 $headline_text->font( $font1, 11 / pt );
273 $headline_text->translate( 15 / mm, 280 / mm );
274 $headline_text->text($headline);
275 $y = $y0 + $dy;
276
277 # AC
278 foreach my $key_ac ( sort keys %OFFSET_SAT ) { #print "$key_ac \n";
279 #SYSTEM
280 foreach my $key_sys ( sort keys %{ $OFFSET_SAT{$key_ac} } ) {
281 my @datasets; # init datasets
282 my $pngName = sprintf ( "%s_OFFSETSAT_%s_%s.png", $inputFilename, $key_ac, $key_sys );
283 my $chartOFFSETSAT = Chart::Gnuplot->new(
284 output => $pngName,
285 terminal => 'png',
286 title => "$key_ac - $key_sys",
287 ylabel => "AC and satellite specific offset [m]", # yrange => [" -0.1 ", " 0.1 "],
288 xlabel => "Time [h]",
289 timeaxis => 'x',
290 xtics => { labelfmt => '%H:%M', rotate => '-270' },
291 legend => { position => "outside below", },
292 grid => 'on',
293 );
294
295 # SATELLITE
296 foreach my $key_sat ( sort keys %{ $OFFSET_SAT{$key_ac}{$key_sys} } ) {
297 my $colour = sprintf ( "%d", ( substr $key_sat, -2 ) + 1 );
298 my $dataset = Chart::Gnuplot::DataSet->new(
299 xdata => $OFFSET_SAT{$key_ac}{$key_sys}{$key_sat}{EPOCH},
300 ydata => $OFFSET_SAT{$key_ac}{$key_sys}{$key_sat}{DATA},
301 title => "$key_sat",
302 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
303 linetype => $colour,
304 style => " linespoints ",
305 );
306 push ( @datasets, $dataset );
307 }
308 $chartOFFSETSAT->plot2d(@datasets);
309 $y = $y - $dy;
310 if ( $y < 30 / mm ) {
311 $page = $pdf->page();
312 $page->mediabox('A4');
313 $y = $y0;
314 }
315 $png = $page->gfx();
316 die ("Unable to find image file: $!") unless -e $pngName;
317 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
318 }
319 }
320}
321######### OFFSET AC #####################
322print "Plot OFFSET AC ...\n";
323$page = $pdf->page();
324$page->mediabox('A4');
325$headline = sprintf ( "Common AC offsets (%s)", $inputFilename );
326$headline_text = $page->text;
327$headline_text->font( $font1, 11 / pt );
328$headline_text->translate( 15 / mm, 280 / mm );
329$headline_text->text($headline);
330$y = $y0 + $dy;
331
332# SYSTEM
333foreach my $key_sys ( sort keys %OFFSET_AC ) {
334 my @datasets; # init datasets
335 my $pngName = sprintf ( "%s_OFFSETAC_%s.png", $inputFilename, $key_sys );
336 my $chart = Chart::Gnuplot->new(
337 output => $pngName,
338 terminal => 'png',
339 title => "$key_sys",
340 ylabel => "Common Offset [m]", #yrange => ["-20.0", "20.0"],
341 xlabel => "Time [h]",
342 timeaxis => 'x',
343 xtics => { labelfmt => '%H:%M', rotate => '-270' },
344 legend => { position => "outside below", },
345 grid => 'on',
346 );
347 # SATELLITE
348 foreach my $key_ac ( sort keys %{ $OFFSET_AC{$key_sys} } ) {
349 my $dataset = Chart::Gnuplot::DataSet->new(
350 xdata => $OFFSET_AC{$key_sys}{$key_ac}{EPOCH},
351 ydata => $OFFSET_AC{$key_sys}{$key_ac}{DATA},
352 title => "$key_ac",
353 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
354 style => " linespoints ",
355 );
356 push ( @datasets, $dataset );
357 }
358 $chart->plot2d(@datasets);
359 $y = $y - $dy;
360 if ( $y < 30 / mm ) {
361 $page = $pdf->page();
362 $page->mediabox('A4');
363 $y = $y0;
364 }
365 $png = $page->gfx();
366 die ("Unable to find image file: $!") unless -e $pngName;
367 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
368}
369######### RES #####################
370print "Plot RES ...\n";
371$page = $pdf->page();
372$page->mediabox('A4');
373$headline = sprintf ( "Residuals (%s)", $inputFilename );
374$headline_text = $page->text;
375$headline_text->font( $font1, 11 / pt );
376$headline_text->translate( 15 / mm, 280 / mm );
377$headline_text->text($headline);
378$y = $y0 + $dy;
379$yrange = 5.0;
380if ( $plotType eq "KF" ) {
381 $yrange = 0.1;
382}
383# AC
384foreach my $key_ac ( sort keys %RES ) { #print "$key_ac \n";
385 #SYSTEM
386 foreach my $key_sys ( sort keys %{ $RES{$key_ac} } ) {
387 my @datasets; # init datasets
388 my $pngName = sprintf ( "%s_RES_%s_%s.png", $inputFilename, $key_ac, $key_sys );
389 my $chart = Chart::Gnuplot->new(
390 output => $pngName,
391 terminal => 'png',
392 title => "$key_ac - $key_sys",
393 ylabel => "Residuals [m]",
394 yrange => [ " -$yrange ", "$yrange " ],
395 xlabel => "Time [h]",
396 timeaxis => 'x',
397 xtics => { labelfmt => '%H:%M', rotate => '-270', },
398 legend => { position => "outside below", },
399 grid => 'on',
400 );
401 # SATELLITE
402 foreach my $key_sat ( sort keys %{ $RES{$key_ac}{$key_sys} } ) {
403 my $colour = sprintf ( "%d", ( substr $key_sat, -2 ) + 1 );
404 my $dataset = Chart::Gnuplot::DataSet->new(
405 xdata => $RES{$key_ac}{$key_sys}{$key_sat}{EPOCH},
406 ydata => $RES{$key_ac}{$key_sys}{$key_sat}{DATA},
407 title => "$key_sat",
408 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
409 style => " linespoints ",
410 );
411 push ( @datasets, $dataset );
412 }
413 $chart->plot2d(@datasets);
414 $y = $y - $dy;
415 if ( $y < 30 / mm ) {
416 $page = $pdf->page();
417 $page->mediabox('A4');
418 $y = $y0;
419 }
420 $png = $page->gfx();
421 die ("Unable to find image file: $!") unless -e $pngName;
422 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
423 }
424}
425
426$pdf->save();
427$pdf->end();
428system ("rm $inputDir/*png");
429system("evince $inputDir/$pdf_name&");
430
431#foreach my $t(@array) {print"$t \n ";}
432#print Dumper \%AMB;
433#########################################
434sub HELP_MESSAGE {
435 print <<EOI_HILFE;
436cmbPlot.pl - Script to pplot BNC's PPP logfiles
437
438CALL:
439 cmbPlot.pl [OPTIONS]
440
441OPTIONS:
442 --plotType KF (Kalman Filter combination) or SE (Single Epoch Combination)
443 --logFiles comma separated list of BNC's combination logfiles
444 --help show help contents
445
446DESCRIPTION:
447
448EOI_HILFE
449 exit;
450}
Note: See TracBrowser for help on using the repository browser.