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

Last change on this file since 10491 was 10491, checked in by stuerze, 7 weeks ago

changes regarding CMB plot

  • Property svn:executable set to *
File size: 18.0 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 if ($val != 0.0) {
108 push ( @{ $CLK_CORR{$sys}{$sat}{EPOCH} }, $epo );
109 push ( @{ $CLK_CORR{$sys}{$sat}{DATA} }, $val );
110 }
111 }
112 elsif ( $_ =~ /Full Clock / ) { # 2015-04-06 23:59:30.000 Full Clock G04 81 -3205.8546
113 @array = split ( /\s+/, $_ );
114 $epo = $array[0] . "T" . $array[1];
115 $sat = $array[4];
116 $sys = substr ( $sat, 0, 1 );
117 $val = $array[6];
118 push ( @{ $CLK_FULL{$sys}{$sat}{EPOCH} }, $epo );
119 push ( @{ $CLK_FULL{$sys}{$sat}{DATA} }, $val );
120 }
121 elsif ( $_ =~ /res = / ) { # 2015-04-05 23:59:40.000 NRCAN G01 res = -0.0003
122 @array = split ( /\s+/, $_ );
123 $epo = $array[0] . "T" . $array[1];
124 $ac = $array[2];
125 $sat = $array[3];
126 $sys = substr ( $sat, 0, 1 );
127 $val = $array[6];
128 push ( @{ $RES{$ac}{$sys}{$sat}{EPOCH} }, $epo );
129 push ( @{ $RES{$ac}{$sys}{$sat}{DATA} }, $val );
130 }
131 elsif ( $_ =~ /AC Offset / ) { # AC Offset DLR C 1.9002 +- 35.3112
132 @array = split ( /\s+/, $_ );
133 $epo = $array[0] . "T" . $array[1];
134 $sys = $array[5];
135 $ac = $array[4];
136 $val = $array[6];
137 if ($val != 0.0) {
138 push ( @{ $OFFSET_AC{$sys}{$ac}{EPOCH} }, $epo );
139 push ( @{ $OFFSET_AC{$sys}{$ac}{DATA} }, $val );
140 }
141 }
142 elsif ( $_ =~ /Sat Offset / ) { # 2015-04-05 23:59:40.000 Sat Offset BKG G01 0.5483 +- 100.0000
143 @array = split ( /\s+/, $_ );
144 $epo = $array[0] . "T" . $array[1];
145 $ac = $array[4];
146 $sat = $array[5];
147 $sys = substr ( $sat, 0, 1 );
148 $val = $array[6];
149 if ($val != 0.0) {
150 push ( @{ $OFFSET_SAT{$ac}{$sys}{$sat}{EPOCH} }, $epo );
151 push ( @{ $OFFSET_SAT{$ac}{$sys}{$sat}{DATA} }, $val );
152 }
153 }
154 }
155 close INPUT;
156} # ----- end foreach logfile -----
157
158# -----------------------------------------------------------------------------
159# plot several data sets
160# -----------------------------------------------------------------------------
161######### CLOCK CORR #####################
162print "Plot CLOCK CORR ...\n";
163$page = $pdf->page();
164$page->mediabox('A4');
165$headline = sprintf ( "Clock Corrections (%s)", $inputFilename );
166$headline_text = $page->text;
167$headline_text->font( $font1, 11 / pt );
168$headline_text->translate( 15 / mm, 280 / mm );
169$headline_text->text($headline);
170$y = $y0 + $dy;
171my $yrange = 20.0;
172if ( $plotType eq "KF" ) {
173 $yrange = 10.0;
174}
175
176# SYSTEM
177foreach my $key_sys ( sort keys %CLK_CORR ) {
178 print " Plot CLOCK CORR for $key_sys\n";
179 my @datasets; # init datasets
180 my $pngName = sprintf ( "%s_CLK_CORR_%s.png", $inputFilename, $key_sys );
181 my $chart = Chart::Gnuplot->new(
182 output => $pngName,
183 terminal => 'png',
184 title => "$key_sys", # G,R
185 ylabel => "Clock corrections [s]",
186 yrange => [ -$yrange, $yrange ],
187 xlabel => "Time [h]",
188 timeaxis => 'x',
189 xtics => { labelfmt => '%H:%M', rotate => '-270' },
190 legend => { position => "outside below", },
191 grid => 'on',
192 );
193
194 # SATELLITE
195 foreach my $key_sat ( sort keys %{ $CLK_CORR{$key_sys} } ) {
196 my $colour = sprintf ( "%d", ( substr $key_sat, -2 ) + 1 );
197 my $dataset = Chart::Gnuplot::DataSet->new(
198 xdata => $CLK_CORR{$key_sys}{$key_sat}{EPOCH},
199 ydata => $CLK_CORR{$key_sys}{$key_sat}{DATA},
200 title => "$key_sat",
201 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
202 linetype => $colour,
203 style => " linespoints ",
204 );
205 push ( @datasets, $dataset );
206 }
207 $chart->plot2d(@datasets);
208 $y = $y - $dy;
209 if ( $y < 30 / mm ) {
210 $page = $pdf->page();
211 $page->mediabox('A4');
212 $y = $y0;
213 }
214 $png = $page->gfx();
215 die ("Unable to find image file: $!") unless -e $pngName;
216 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
217}
218
219######### CLOCK FULL #####################
220print "Plot FULL CLOCK ...\n";
221$page = $pdf->page();
222$page->mediabox('A4');
223$headline = sprintf ( "Full Clock (%s)", $inputFilename );
224$headline_text = $page->text;
225$headline_text->font( $font1, 11 / pt );
226$headline_text->translate( 15 / mm, 280 / mm );
227$headline_text->text($headline);
228$y = $y0 + $dy;
229
230# SYSTEM
231foreach my $key_sys ( sort keys %CLK_FULL ) {
232 my @datasets; # init datasets
233 my $pngName = sprintf ( "%s_CLK_FULL_%s.png", $inputFilename, $key_sys );
234 my $chart = Chart::Gnuplot->new(
235 output => $pngName,
236 terminal => 'png',
237 title => "$key_sys",
238 ylabel => "Full Clock [m]", #yrange => [" -3 ", " 3 "],
239 xlabel => "Time [h]",
240 timeaxis => 'x',
241 xtics => { labelfmt => '%H:%M', rotate => '-270', },
242 legend => { position => "outside below", },
243 grid => 'on',
244 );
245
246 # SATELLITE
247 foreach my $key_sat ( sort keys %{ $CLK_FULL{$key_sys} } ) {
248 my $colour = sprintf ( "%d", ( substr $key_sat, -2 ) + 1 );
249 my $dataset = Chart::Gnuplot::DataSet->new(
250 xdata => $CLK_FULL{$key_sys}{$key_sat}{EPOCH},
251 ydata => $CLK_FULL{$key_sys}{$key_sat}{DATA},
252 title => "$key_sat",
253 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
254 linetype => $colour,
255 style => " linespoints ",
256 );
257 push ( @datasets, $dataset );
258 }
259 $chart->plot2d(@datasets);
260 $y = $y - $dy;
261 if ( $y < 30 / mm ) {
262 $page = $pdf->page();
263 $page->mediabox('A4');
264 $y = $y0;
265 }
266 $png = $page->gfx();
267 die ("Unable to find image file: $!") unless -e $pngName;
268 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
269}
270
271######### OFFSET SAT #####################
272if ( $plotType eq "KF" ) {
273 print "Plot OFFSET SAT ...\n";
274 $page = $pdf->page();
275 $page->mediabox('A4');
276 $headline = sprintf ( "AC and satellite specific offset (%s)", $inputFilename );
277 $headline_text = $page->text;
278 $headline_text->font( $font1, 11 / pt );
279 $headline_text->translate( 15 / mm, 280 / mm );
280 $headline_text->text($headline);
281 $y = $y0 + $dy;
282
283 # AC
284 foreach my $key_ac ( sort keys %OFFSET_SAT ) { #print "$key_ac \n";
285 #SYSTEM
286 foreach my $key_sys ( sort keys %{ $OFFSET_SAT{$key_ac} } ) {
287 my @datasets; # init datasets
288 my $pngName = sprintf ( "%s_OFFSETSAT_%s_%s.png", $inputFilename, $key_ac, $key_sys );
289 my $chartOFFSETSAT = Chart::Gnuplot->new(
290 output => $pngName,
291 terminal => 'png',
292 title => "$key_ac - $key_sys",
293 ylabel => "AC and satellite specific offset [m]", # yrange => [" -0.1 ", " 0.1 "],
294 xlabel => "Time [h]",
295 timeaxis => 'x',
296 xtics => { labelfmt => '%H:%M', rotate => '-270' },
297 legend => { position => "outside below", },
298 grid => 'on',
299 );
300
301 # SATELLITE
302 foreach my $key_sat ( sort keys %{ $OFFSET_SAT{$key_ac}{$key_sys} } ) {
303 my $colour = sprintf ( "%d", ( substr $key_sat, -2 ) + 1 );
304 my $dataset = Chart::Gnuplot::DataSet->new(
305 xdata => $OFFSET_SAT{$key_ac}{$key_sys}{$key_sat}{EPOCH},
306 ydata => $OFFSET_SAT{$key_ac}{$key_sys}{$key_sat}{DATA},
307 title => "$key_sat",
308 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
309 linetype => $colour,
310 style => " linespoints ",
311 );
312 push ( @datasets, $dataset );
313 }
314 $chartOFFSETSAT->plot2d(@datasets);
315 $y = $y - $dy;
316 if ( $y < 30 / mm ) {
317 $page = $pdf->page();
318 $page->mediabox('A4');
319 $y = $y0;
320 }
321 $png = $page->gfx();
322 die ("Unable to find image file: $!") unless -e $pngName;
323 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
324 }
325 }
326}
327######### OFFSET AC #####################
328print "Plot OFFSET AC ...\n";
329$page = $pdf->page();
330$page->mediabox('A4');
331$headline = sprintf ( "Common AC offsets (%s)", $inputFilename );
332$headline_text = $page->text;
333$headline_text->font( $font1, 11 / pt );
334$headline_text->translate( 15 / mm, 280 / mm );
335$headline_text->text($headline);
336$y = $y0 + $dy;
337
338# SYSTEM
339foreach my $key_sys ( sort keys %OFFSET_AC ) {
340 my @datasets; # init datasets
341 my $pngName = sprintf ( "%s_OFFSETAC_%s.png", $inputFilename, $key_sys );
342 my $chart = Chart::Gnuplot->new(
343 output => $pngName,
344 terminal => 'png',
345 title => "$key_sys",
346 ylabel => "Common Offset [m]", #yrange => ["-20.0", "20.0"],
347 xlabel => "Time [h]",
348 timeaxis => 'x',
349 xtics => { labelfmt => '%H:%M', rotate => '-270' },
350 legend => { position => "outside below", },
351 grid => 'on',
352 );
353 # SATELLITE
354 foreach my $key_ac ( sort keys %{ $OFFSET_AC{$key_sys} } ) {
355 my $dataset = Chart::Gnuplot::DataSet->new(
356 xdata => $OFFSET_AC{$key_sys}{$key_ac}{EPOCH},
357 ydata => $OFFSET_AC{$key_sys}{$key_ac}{DATA},
358 title => "$key_ac",
359 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
360 style => " linespoints ",
361 );
362 push ( @datasets, $dataset );
363 }
364 $chart->plot2d(@datasets);
365 $y = $y - $dy;
366 if ( $y < 30 / mm ) {
367 $page = $pdf->page();
368 $page->mediabox('A4');
369 $y = $y0;
370 }
371 $png = $page->gfx();
372 die ("Unable to find image file: $!") unless -e $pngName;
373 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
374}
375######### RES #####################
376print "Plot RES ...\n";
377$page = $pdf->page();
378$page->mediabox('A4');
379$headline = sprintf ( "Residuals (%s)", $inputFilename );
380$headline_text = $page->text;
381$headline_text->font( $font1, 11 / pt );
382$headline_text->translate( 15 / mm, 280 / mm );
383$headline_text->text($headline);
384$y = $y0 + $dy;
385$yrange = 5.0;
386if ( $plotType eq "KF" ) {
387 $yrange = 0.1;
388}
389# AC
390foreach my $key_ac ( sort keys %RES ) { #print "$key_ac \n";
391 #SYSTEM
392 foreach my $key_sys ( sort keys %{ $RES{$key_ac} } ) {
393 my @datasets; # init datasets
394 my $pngName = sprintf ( "%s_RES_%s_%s.png", $inputFilename, $key_ac, $key_sys );
395 my $chart = Chart::Gnuplot->new(
396 output => $pngName,
397 terminal => 'png',
398 title => "$key_ac - $key_sys",
399 ylabel => "Residuals [m]",
400 yrange => [ " -$yrange ", "$yrange " ],
401 xlabel => "Time [h]",
402 timeaxis => 'x',
403 xtics => { labelfmt => '%H:%M', rotate => '-270', },
404 legend => { position => "outside below", },
405 grid => 'on',
406 );
407 # SATELLITE
408 foreach my $key_sat ( sort keys %{ $RES{$key_ac}{$key_sys} } ) {
409 my $colour = sprintf ( "%d", ( substr $key_sat, -2 ) + 1 );
410 my $dataset = Chart::Gnuplot::DataSet->new(
411 xdata => $RES{$key_ac}{$key_sys}{$key_sat}{EPOCH},
412 ydata => $RES{$key_ac}{$key_sys}{$key_sat}{DATA},
413 title => "$key_sat",
414 timefmt => '%Y-%m-%dT%H:%M:%S', # '%H:%M:%S',
415 style => " linespoints ",
416 );
417 push ( @datasets, $dataset );
418 }
419 $chart->plot2d(@datasets);
420 $y = $y - $dy;
421 if ( $y < 30 / mm ) {
422 $page = $pdf->page();
423 $page->mediabox('A4');
424 $y = $y0;
425 }
426 $png = $page->gfx();
427 die ("Unable to find image file: $!") unless -e $pngName;
428 $png->image( $pdf->image_png($pngName), $x, $y, $width, $height );
429 }
430}
431
432$pdf->save();
433$pdf->end();
434system ("rm $inputDir/*png");
435system("evince $inputDir/$pdf_name&");
436
437#foreach my $t(@array) {print"$t \n ";}
438#print Dumper \%AMB;
439#########################################
440sub HELP_MESSAGE {
441 print <<EOI_HILFE;
442cmbPlot.pl - Script to pplot BNC's PPP logfiles
443
444CALL:
445 cmbPlot.pl [OPTIONS]
446
447OPTIONS:
448 --plotType KF (Kalman Filter combination) or SE (Single Epoch Combination)
449 --logFiles comma separated list of BNC's combination logfiles
450 --help show help contents
451
452DESCRIPTION:
453
454EOI_HILFE
455 exit;
456}
Note: See TracBrowser for help on using the repository browser.