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

Last change on this file since 9836 was 9836, checked in by wiese, 19 months ago

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