source: ntrip/trunk/BNC/scripts/compareSP3Plot.pl@ 10076

Last change on this file since 10076 was 9864, checked in by stuerze, 2 years ago

minor changes

  • Property svn:executable set to *
File size: 14.3 KB
RevLine 
[9830]1#!/usr/bin/env perl
2# ========================================================================
[9864]3# compareSP3Plot.pl
[9830]4# ========================================================================
5#
6# Purpose: Script using gnuplot to plot BNC's SP3 Comparison Results
7#
8# Comment: stored values:
9# To compare satellite clocks provided by the two files, BNC first converts
10# coordinate differences dX,dY,dZ into along track, out-of-plane, and radial
11# components.
12# It then corrects the clock differences for the radial components of
13# coordinate differences. RMS values of clock differences are
14# finally calculated after introducing
15# - at first one offset 'per epoch for all satellites' and
16# - secondly one offset 'per satellite for all epochs'.
17# ========================================================================
18# Uses
19use strict;
20use warnings;
21use Getopt::Long;
22use Chart::Gnuplot;
23use Data::Dumper qw(Dumper);
24use diagnostics;
25use File::Basename;
26use Date::Manip;
[9864]27use Time::Piece 1.30;
[9830]28use List::MoreUtils qw( minmax );
29use experimental 'smartmatch';
30
31# -----------------------------------------------------------------------------
32# required to create a pdf
33# -----------------------------------------------------------------------------
34use PDF::API2;
35use constant {mm => 25.4 / 72, inch => 1 / 72,pt => 1,
36};# There are 72 postscript points in an inch and there are 25.4 millimeters in an inch.
37# -----------------------------------------------------------------------------
38# Options
39# -----------------------------------------------------------------------------
40# init options
41my $help = 0;
42my @plotTypes = ();
43my @logFiles = ();
44# read options from command line
45GetOptions(
46 'help' => \$help,
47 'plotTypes=s' => \@plotTypes,
48 'logFiles=s' => \@logFiles,
49);
50help() if $help;
51@plotTypes = qw(NEU) unless (@plotTypes);# Default values
52@plotTypes = map {uc} split ( /[, ]/, join ( ',', @plotTypes ) );# Reformat options
53@logFiles = split ( /[, ]/, join ( ',', @logFiles ) );
54unless (@logFiles) {
55 print "ERROR: logfiles missing\n";
56 exit;
57}
58print " plotTpes: @plotTypes\n logfiles: @logFiles\n\n";
59# -----------------------------------------------------------------------------
60# generate data sets for gnuplot
61# -----------------------------------------------------------------------------
[9864]62my $epoSec = 0;
[9830]63# for pdf generation
64my ($png, $page, $headline, $headline_text);
65my $y0 = 190/mm;
66my ($x, $y, $width, $height) = (40/mm, $y0, 130/mm, 70/mm);
67my $dy = $height+7/mm;
68# file by file ..
69foreach my $file (@logFiles) {
70 my (%COMPARISON);
71 my (@array);
[9864]72 my ($sys, $sat, $epo, $radial, $along, $out, $drange, $clk, $clkRed);
[9830]73 # -----------------------------------------------------------------------------
74 # create pdf for plot results
75 # -----------------------------------------------------------------------------
76 my($inputFilename, $inputDir, $inputSuffix) = fileparse($file, '\..*');
77 my $pdf_name = sprintf("%s.pdf", $inputFilename);
78 my $pdf = PDF::API2->new(-file => "$inputDir$pdf_name");
79 my $font1 = $pdf->corefont('Helvetica-Bold');
80 # -----------------------------------------------------------------------------
81 # read file and save data
82 # -----------------------------------------------------------------------------
83 open ( INPUT, "<$file" ) || die $!;
84 print "Parse logfile $file ...\n";
85 while (<INPUT>) {#56985.000000 G01 0.0160 -0.0242 0.0961 0.0082 -0.0077 1
86 if ( $_ !~ /^!/) { #print "$_\n";
87 @array = split(/\s+/, $_);
[9864]88 $epo = $array[0]; #print "$epo\n";
89 $epo =~s/(\.\d*)//; #print "$epo\n";
90 $sys = substr $array[1], 0, 1;
[9830]91 $sat = $array[1];#print "$sat\n";
92 $radial = $array[2];#print "$radial\n";
93 $along = $array[3];#print "$along\n";
94 $out = $array[4];#print "$out\n";
95 $clk = $array[5];#print "$clk\n";
96 $clkRed = $array[6];#print "$clkRed\n";
97 if ($clk !~ /^\.$/) {
98 push @{$COMPARISON{$sys}{$sat}{EPOCH }}, $epo;
99 push @{$COMPARISON{$sys}{$sat}{RADIAL}}, $radial;
100 push @{$COMPARISON{$sys}{$sat}{ALONG }}, $along;
101 push @{$COMPARISON{$sys}{$sat}{OUT }}, $out;
102 push @{$COMPARISON{$sys}{$sat}{CLK }}, $clk;
103 push @{$COMPARISON{$sys}{$sat}{CLKRED}}, $clkRed;
104 }
105 }
106 }
107 close INPUT;
108
109 # -----------------------------------------------------------------------------
110 # plot several data sets
111 # -----------------------------------------------------------------------------
112 ######### RAO #####################
113 if ("ALL" ~~ @plotTypes ) {
114 print "Plot RADIAL, ALONG TRACK and OUT OF PLANE COMPONENTS ...\n";
115 my $yrange = .6;
116 $page = $pdf->page(); $page->mediabox('A4');
117 $headline = sprintf("Orbit Differences in radial, along-track and out-of-plane component (%s)", $inputFilename);
118 $headline_text = $page->text;
119 $headline_text->font($font1, 11/pt);
120 $headline_text->translate(15/mm, 280/mm);
121 $headline_text->text($headline);
122 $y = $y0+$dy;
123 #SYSTEM
124 foreach my $key_sys (sort keys %COMPARISON) {
125 print " for $key_sys\n";
126 my (@datasets_r, @datasets_a, @datasets_o); # init datasets
127 my $pngName_r = sprintf ("%s_RADIAL_%s.png", $inputFilename, $key_sys);
128 my $pngName_a = sprintf ("%s_ALONG_%s.png", $inputFilename, $key_sys);
129 my $pngName_o = sprintf ("%s_OUT_%s.png", $inputFilename, $key_sys);
130 my $chart_r = Chart::Gnuplot->new(
131 output => $pngName_r, terminal => 'png', title => "Radial component",
132 ylabel => "Orbit Difference [m]", yrange => ["-$yrange", "$yrange"],
133 xlabel => "Time [h]", timeaxis => 'x', xtics => {labelfmt => '%H:%M', rotate => '-270',},
134 legend => { position => "outside right",},
135 grid => 'on',);
136 my $chart_a = Chart::Gnuplot->new(
137 output => $pngName_a, terminal => 'png', title => "Along track component",
138 ylabel => "Orbit Difference [m]", yrange => ["-$yrange", "$yrange"],
139 xlabel => "Time [h]", timeaxis => 'x', xtics => {labelfmt => '%H:%M', rotate => '-270',},
140 legend => { position => "outside right",},
141 grid => 'on',);
142 my $chart_o = Chart::Gnuplot->new(
143 output => $pngName_o, terminal => 'png', title => "Out of plane component",
144 ylabel => "Orbit Difference [m]", yrange => ["-$yrange", "$yrange"],
145 xlabel => "Time [h]", timeaxis => 'x', xtics => {labelfmt => '%H:%M', rotate => '-270',},
146 legend => { position => "outside right",},
147 grid => 'on',);
148 #SATELLITE
149 foreach my $key_sat (sort keys %{$COMPARISON{$key_sys}}) {
150 my $colour = sprintf("%d", (substr $key_sat, -2)+1);
151 #RADIAL
152 my $dataset_r = Chart::Gnuplot::DataSet->new(
153 xdata => \@{$COMPARISON{$key_sys}{$key_sat}{EPOCH}},
154 ydata => \@{$COMPARISON{$key_sys}{$key_sat}{RADIAL}},
155 title => "$key_sat",
156 timefmt => '%Y-%m-%d_%H:%M:%S',
157 linetype => $colour,
158 #style => " linespoints ",);
159 style => " lines ",);
160 push (@datasets_r, $dataset_r);
161 # ALONG TRACK
162 my $dataset_a = Chart::Gnuplot::DataSet->new(
163 xdata => \@{$COMPARISON{$key_sys}{$key_sat}{EPOCH}},
164 ydata => \@{$COMPARISON{$key_sys}{$key_sat}{ALONG}},
165 title => "$key_sat",
166 timefmt => '%Y-%m-%d_%H:%M:%S',
167 linetype => $colour,
168 style => " lines ",);
169 push (@datasets_a, $dataset_a);
170 # OUT OF PLANE
171 my $dataset_o = Chart::Gnuplot::DataSet->new(
172 xdata => \@{$COMPARISON{$key_sys}{$key_sat}{EPOCH}},
173 ydata => \@{$COMPARISON{$key_sys}{$key_sat}{OUT}},
174 title => "$key_sat",
175 timefmt => '%Y-%m-%d_%H:%M:%S',
176 linetype => $colour,
177 style => " lines",);
178 push (@datasets_o, $dataset_o);
179 }
180 $chart_r->plot2d(@datasets_r);
181 $chart_a->plot2d(@datasets_a);
182 $chart_o->plot2d(@datasets_o);
183 # RADIAL
184 $y = $y - $dy;
185 if($y < 30/mm) {
186 $page = $pdf->page(); $page->mediabox('A4');
187 $y = $y0;
188 }
189 $png = $page->gfx(); die("Unable to find image file: $!") unless -e $pngName_r;
190 $png->image($pdf->image_png($pngName_r), $x, $y, $width, $height);
191 # ALONG
192 $y = $y - $dy;
193 if($y < 30/mm) {
194 $page = $pdf->page(); $page->mediabox('A4');
195 $y = $y0;
196 }
197 $png = $page->gfx(); die("Unable to find image file: $!") unless -e $pngName_a;
198 $png->image($pdf->image_png($pngName_a), $x, $y, $width, $height);
199 # OUT
200 $y = $y - $dy;
201 if($y < 30/mm) {
202 $page = $pdf->page(); $page->mediabox('A4');
203 $y = $y0;
204 }
205 $png = $page->gfx(); die("Unable to find image file: $!") unless -e $pngName_o;
206 $png->image($pdf->image_png($pngName_o), $x, $y, $width, $height);
207 }
208 }
209 ######## CLK / CLKRED #####################
210 if (("ALL" ~~ @plotTypes ) || ("CLK" ~~ @plotTypes)) {
211 print "Plot CLOCK and CLOCK REDUCED BY RADIAL COMPONENT ...\n";
[9864]212 my $yrange = 0.6;
[9830]213 $page = $pdf->page(); $page->mediabox('A4');
214 $headline = sprintf("Clock residuals and reduced clock residuals(%s)", $inputFilename);
215 $headline_text = $page->text;
216 $headline_text->font($font1, 11/pt);
217 $headline_text->translate(15/mm, 280/mm);
218 $headline_text->text($headline);
219 $y = $y0+$dy;
220 #SYSTEM
221 foreach my $key_sys (sort keys %COMPARISON) {
222 print " for $key_sys\n";
223 my (@datasets_c, @datasets_r); # init datasets
224 my $pngName_c = sprintf ("%s_CLK_%s.png", $inputFilename, $key_sys);
225 my $pngName_r = sprintf ("%s_CLKRED_%s.png", $inputFilename, $key_sys);
226 my $chart_c = Chart::Gnuplot->new(
227 output => $pngName_c, terminal => 'png',title => "Clock only",
[9864]228 ylabel => "Clock Residuals [m]", yrange => ["-$yrange", "$yrange"],
[9830]229 xlabel => "Time [h]", timeaxis => 'x', xtics => {labelfmt => '%H:%M', rotate => '-270',},
230 legend => { position => "outside right",},
231 grid => 'on',);
232 my $chart_r = Chart::Gnuplot->new(
233 output => $pngName_r, terminal => 'png',title => "Clock minus radial component",
[9864]234 ylabel => "Clock Residuals [m]", yrange => ["-$yrange", "$yrange"],
[9830]235 xlabel => "Time [h]", timeaxis => 'x', xtics => {labelfmt => '%H:%M', rotate => '-270',},
236 legend => { position => "outside right",},
237 grid => 'on',);
238 #SATELLITE
239 foreach my $key_sat (sort keys %{$COMPARISON{$key_sys}}) {
240 my $colour = sprintf("%d", (substr $key_sat, -2)+1);
241 # CLK
242 my $dataset_c = Chart::Gnuplot::DataSet->new(
243 xdata => \@{$COMPARISON{$key_sys}{$key_sat}{EPOCH}},
244 ydata => \@{$COMPARISON{$key_sys}{$key_sat}{CLK}},
245 title => "$key_sat",
246 timefmt => '%Y-%m-%d_%H:%M:%S',
247 linetype => $colour,
248 style => " lines ",);
249 push (@datasets_c, $dataset_c);
250 my ( $min, $max ) = minmax @{$COMPARISON{$key_sys}{$key_sat}{CLK}};
251 if ($max > 1.0) {print "$key_sat: max clk = $max\n"};
252 # CLKRED
253 my $dataset_r = Chart::Gnuplot::DataSet->new(
254 xdata => \@{$COMPARISON{$key_sys}{$key_sat}{EPOCH}},
255 ydata => \@{$COMPARISON{$key_sys}{$key_sat}{CLKRED}},
256 title => "$key_sat",
257 timefmt => '%Y-%m-%d_%H:%M:%S',
258 linetype => $colour,
259 style => " lines ",);
260 push (@datasets_r, $dataset_r);
261 }
262 $chart_c->plot2d(@datasets_c);
263 $chart_r->plot2d(@datasets_r);
264 # CLK
265 $y = $y - $dy;
266 if($y < 30/mm) {
267 $page = $pdf->page(); $page->mediabox('A4');
268 $y = $y0;
269 }
270 $png = $page->gfx(); die("Unable to find image file: $!") unless -e $pngName_c;
271 $png->image($pdf->image_png($pngName_c), $x, $y, $width, $height);
272 # CLKRED
273 $y = $y - $dy;
274 if($y < 30/mm) {
275 $page = $pdf->page(); $page->mediabox('A4');
276 $y = $y0;
277 }
278 $png = $page->gfx(); die("Unable to find image file: $!") unless -e $pngName_r;
279 $png->image($pdf->image_png($pngName_r), $x, $y, $width, $height);
280 }
281 }
282 $pdf->save();
283 $pdf->end();
284 system("rm *png");
[9864]285 #system("evince $inputDir/$pdf_name&");
[9830]286}
287
288
289#########################################
290sub help {
291 print <<EOI_HILFE;
292CALL:
293 compareSp3Plot.pl [paramter]
294
295PARAMETER:
296
297 --plotTypes possible: ALL => RAO and CLOCK; CLK => CLOCK and CLOCKRED
298 --logFiles comma separated list of BNC's combination logfiles
299 --help show help contents
300
301DESCRIPTION:
302
[9864]303 Script to plot BNC's SP3 Comparison Results
[9830]304
305EOI_HILFE
306exit;
307}
Note: See TracBrowser for help on using the repository browser.