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

Last change on this file since 9830 was 9830, checked in by stuerze, 19 months ago

minor changes within plot scripts
  • Property svn:executable set to *
File size: 17.9 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 = 5.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 right", },
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 right", },
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 right", },
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 right", },
345 grid => 'on',
346 );
347
348 # SATELLITE
349 foreach my $key_ac ( sort keys %{ $OFFSET_AC{$key_sys} } ) {
350 my $dataset = Chart::Gnuplot::DataSet->new(
351 xdata => $OFFSET_AC{$key_sys}{$key_ac}{EPOCH},
352 ydata => $OFFSET_AC{$key_sys}{$key_ac}{DATA},
353 title => "$key_ac",
354 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
355 style => " linespoints ",
356 );
357
358 push ( @datasets, $dataset );
359 }
360 $chart->plot2d(@datasets);
361 $y = $y - $dy;
362 if ( $y < 30 / mm ) {
363 $page = $pdf->page();
364 $page->mediabox('A4');
365 $y = $y0;
366 }
367 $png = $page->gfx();
368 die ("Unable to find image file: $!") unless -e $pngName;
369 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
370}
371
372######### RES #####################
373print "Plot RES ...\n";
374$page = $pdf->page();
375$page->mediabox('A4');
376$headline = sprintf ( "Residuals (%s)", $inputFilename );
377$headline_text = $page->text;
378$headline_text->font( $font1, 11 / pt );
379$headline_text->translate( 15 / mm, 280 / mm );
380$headline_text->text($headline);
381$y = $y0 + $dy;
382$yrange = 4.0;
383if ( $plotType eq "KF" ) {
384 $yrange = 0.1;
385}
386
387# AC
388foreach my $key_ac ( sort keys %RES ) { #print "$key_ac \n";
389 #SYSTEM
390 foreach my $key_sys ( sort keys %{ $RES{$key_ac} } ) {
391 my @datasets; # init datasets
392 my $pngName = sprintf ( "%s_RES_%s_%s.png", $inputFilename, $key_ac, $key_sys );
393 my $chart = Chart::Gnuplot->new(
394 output => $pngName,
395 terminal => 'png',
396 title => "$key_ac - $key_sys",
397 ylabel => "Residuals [m]",
398 yrange => [ " -$yrange ", "$yrange " ],
399 xlabel => "Time [h]",
400 timeaxis => 'x',
401 xtics => { labelfmt => '%H:%M', rotate => '-270', },
402 legend => { position => "outside right", },
403 grid => 'on',
404 );
405
406 # SATELLITE
407 foreach my $key_sat ( sort keys %{ $RES{$key_ac}{$key_sys} } ) {
408 my $colour = sprintf ( "%d", ( substr $key_sat, -2 ) + 1 );
409 my $dataset = Chart::Gnuplot::DataSet->new(
410 xdata => $RES{$key_ac}{$key_sys}{$key_sat}{EPOCH},
411 ydata => $RES{$key_ac}{$key_sys}{$key_sat}{DATA},
412 title => "$key_sat",
413 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
414 style => " linespoints ",
415 );
416 push ( @datasets, $dataset );
417 }
418 $chart->plot2d(@datasets);
419 $y = $y - $dy;
420 if ( $y < 30 / mm ) {
421 $page = $pdf->page();
422 $page->mediabox('A4');
423 $y = $y0;
424 }
425 $png = $page->gfx();
426 die ("Unable to find image file: $!") unless -e $pngName;
427 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
428 }
429}
430
431$pdf->save();
432$pdf->end();
433system ("rm $inputDir/*png");
434system("evince $inputDir/$pdf_name&");
435
436#foreach my $t(@array) {print"$t \n ";}
437#print Dumper \%AMB;
438#########################################
439sub HELP_MESSAGE {
440 print <<EOI_HILFE;
441cmbPlot.pl - Script to pplot BNC's PPP logfiles
442
443CALL:
444 cmbPlot.pl [OPTIONS]
445
446OPTIONS:
447 --plotType KF (Kalman Filter combination) or SE (Single Epoch Combination)
448 --logFiles comma separated list of BNC's combination logfiles
449 --help show help contents
450
451DESCRIPTION:
452
453EOI_HILFE
454 exit;
455}
Note: See TracBrowser for help on using the repository browser.