#!/usr/bin/env perl # ======================================================================== # pppPlot.pl - plot BNC's PPP results using gnuplot # ======================================================================== # # Plot metrics: # - NEU displacements w.r.t. coordinates in Crd file # - A-priori + correction values of tropospheric zenith delay in [m], # - GNSS Receiver clock errors in [m], # - Elevations, given per satellite # - Ambiguities, given per satellite # - Ionosphere Delay [m], # - Receiver Bias [m] # - Code and phase residuals in [m], # # Author : Andrea Stuerze # Revision: $Header: trunk/BNC/scripts/pppPlot.pl 9947 2022-12-19 11:58:57Z stuerze $ # Changes : # ======================================================================== # Uses use strict; use warnings; BEGIN { use FindBin qw($Bin); use lib "$Bin"; } use diagnostics; use PDL; use FindBin qw($Bin); use Getopt::Long; use Chart::Gnuplot; use Data::Dumper qw(Dumper); use File::Basename; use Date::Manip; use Log::Log4perl qw(:easy); use PDF::API2; use Bnc; use Common; use constant { mm => 25.4 / 72, inch => 1 / 72, pt => 1, }; # There are 72 postscript points in an inch and there are 25.4 millimeters in an inch. # Logging Log::Log4perl->easy_init( { #file => "plotPPP.log", layout => '%d [%c l%L] %p: %m%n', level => $TRACE } ); # Options my ($prog) = fileparse($0); my $help = 0; my @plotTypes = (); my @logFiles = (); my $sampling = 1; GetOptions( 'help' => \$help, 'plotTypes=s' => \@plotTypes, 'logFiles=s' => \@logFiles, 'sampling=s' => \$sampling, ); HELP_MESSAGE() if $help; @plotTypes = qw(NEU) unless (@plotTypes); @plotTypes = map {uc} split ( /[, ]/, join ( ',', @plotTypes ) ); @logFiles = split ( /[, ]/, join ( ',', @logFiles ) ); unless (@logFiles) { ERROR "logfiles missing"; HELP_MESSAGE() } DEBUG("\n plotTpes: @plotTypes\n logfiles: @logFiles\n sampling: $sampling"); # ----------------------------------------------------------------------------- # Generate data sets for gnuplot # ----------------------------------------------------------------------------- # for pdf gerneration my ( $png, $page, $headline, $headline_text ); my $y0 = 180 / mm; my ( $x, $y, $width, $height ) = ( 40 / mm, $y0, 130 / mm, 80 / mm ); my $dy = $height + 10 / mm; # Loop over logfiles foreach my $file (@logFiles) { DEBUG "Parse logfile $file"; # ----------------------------------------------------------------------------- # Create pdf for plot results # ----------------------------------------------------------------------------- my ( $inputFilename, $inputDir, $inputSuffix ) = fileparse( $file, '\..*' ); my $pdf_name = sprintf ( "%s.pdf", $inputFilename ); my $pdf = PDF::API2->new( -file => "$inputDir$pdf_name" ); my $font1 = $pdf->corefont('Helvetica-Bold'); # ----------------------------------------------------------------------------- # Read logfile # ----------------------------------------------------------------------------- my ( $station, $file ) = Bnc::parsePPPLogfile( $file, $sampling ); my $EPOCHS = $file->{'EPOCHS'}; my $EPOCHS_G_CLK = $file->{'EPOCHS_G_CLK'}; my $EPOCHS_R_CLK = $file->{'EPOCHS_R_CLK'}; my $EPOCHS_E_CLK = $file->{'EPOCHS_E_CLK'}; my $EPOCHS_C_CLK = $file->{'EPOCHS_C_CLK'}; my %AMB = %{ $file->{'AMB'} }; my %RES = %{ $file->{'RES'} }; my %ELE = %{ $file->{'ELE'} }; my %ION = %{ $file->{'ION'} }; my %BIA = %{ $file->{'BIA'} }; # ----------------------------------------------------------------------------- # RMS computation # ----------------------------------------------------------------------------- my ( $mean, $prms, $median, $min, $max, $adev, $rms_n, $rms_e, $rms_u, $rms_trp ); my ( $n, $e, $u, $trp, $str_rms_n, $str_rms_e, $str_rms_u, $str_rms_trp ); $n = pdl( $file->{'N'} ); ( $mean, $prms, $median, $min, $max, $adev, $rms_n ) = stats($n); $e = pdl( $file->{'E'} ); ( $mean, $prms, $median, $min, $max, $adev, $rms_e ) = stats($e); $u = pdl( $file->{'U'} ); ( $mean, $prms, $median, $min, $max, $adev, $rms_u ) = stats($u); $trp = pdl( $file->{'TRPs'} ); ( $mean, $prms, $median, $min, $max, $adev, $rms_trp ) = stats($trp); $str_rms_n = sprintf ( " %.2f ", $rms_n ); $str_rms_e = sprintf ( " %.2f ", $rms_e ); $str_rms_u = sprintf ( " %.2f ", $rms_u ); $str_rms_trp = sprintf ( " %.2f ", $rms_trp ); DEBUG("RMS: North: $str_rms_n, East: $str_rms_e, Up: $str_rms_u, TRP: $str_rms_trp"); # ----------------------------------------------------------------------------- # Plot several data sets # ----------------------------------------------------------------------------- my $dataset; ######### NEU ##################### DEBUG "Plot NEU"; $page = $pdf->page(); $page->mediabox('A4'); $headline = sprintf ( "PPP results for station %s", $station ); $headline_text = $page->text; $headline_text->font( $font1, 11 / pt ); $headline_text->translate( 15 / mm, 280 / mm ); $headline_text->text($headline); $y = $y0; my $pngName = sprintf ( "%s_NEU.png", $station ); my $chartNEU = newChart($station); $chartNEU->set( output => $pngName ); $chartNEU->set( ylabel => "Displacements [m]", yrange => [ " -0.5 ", " 0.5 " ] ); #y2label => "Number of Satellites [-]", y2range => [" 0 ", " 20 "], y2tics => 'on', my $dataN = Chart::Gnuplot::DataSet->new( xdata => $EPOCHS, ydata => $file->{'N'}, title => "Displacements N, RMS + -$str_rms_n m", timefmt => '%s', style => "linespoints", ); my $dataE = Chart::Gnuplot::DataSet->new( xdata => $EPOCHS, ydata => $file->{'E'}, title => "Displacements E, RMS + -$str_rms_e m", timefmt => '%s', style => "linespoints", ); my $dataU = Chart::Gnuplot::DataSet->new( xdata => $EPOCHS, ydata => $file->{'U'}, title => "Displacements U, RMS + -$str_rms_u m", timefmt => '%s', style => "linespoints", ); my @datasets = ( $dataN, $dataE, $dataU ); $chartNEU->plot2d(@datasets); $png = $page->gfx(); LOGDIE("could not find image file: $!\n") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); #print "y= : $y \n" ######### TRP ##################### if ( grep ( $_ eq "ALL", @plotTypes ) ) { DEBUG "Plot TRPs"; my $pngName = sprintf ( "%s_TRP.png", $station ); my $chartTRP = newChart($station); $chartTRP->set( output => $pngName ); $chartTRP->set( ylabel => "Tropospheric Delay [m]", yrange => [ " 1.0 ", " 3.0 " ] ); my $dataTRP = Chart::Gnuplot::DataSet->new( xdata => $EPOCHS, ydata => $file->{'TRPs'}, title => "Tropospheric Delay, RMS + -$str_rms_trp m", timefmt => '%s', style => "linespoints", ); $chartTRP->plot2d($dataTRP); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); LOGDIE("could not find image file: $!\n") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); ######### GPS CLK ##################### if ( scalar @{ $file->{'G_CLKs'} } < 1 ) { DEBUG "No GPS CLKs found"; } else { DEBUG "Plot GPS CLKs"; $pngName = sprintf ( "%s_CLK_G.png", $station ); my $chartCLK = newChart($station); $chartCLK->set( output => $pngName ); $chartCLK->set( ylabel => "CLK [m]" ); $dataset = Chart::Gnuplot::DataSet->new( xdata => $EPOCHS_G_CLK, ydata => $file->{'G_CLKs'}, title => "GPS Receiver clock", timefmt => '%s', style => "linespoints", ); $chartCLK->plot2d($dataset); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); LOGDIE("could not find image file: $!\n") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); } ######### GLONASS CLK ##################### if ( scalar @{ $file->{'R_CLKs'} } < 1 ) { DEBUG "No GLONASS CLKs found"; } else { DEBUG "Plot GLONASS CLKs"; $pngName = sprintf ( "%s_CLK_R.png", $station ); my $chartCLK = newChart($station); $chartCLK->set( output => $pngName ); $chartCLK->set( ylabel => "CLK [m]" ); $dataset = Chart::Gnuplot::DataSet->new( xdata => $EPOCHS_R_CLK, ydata => $file->{'R_CLKs'}, title => "GLONASS Receiver clock", timefmt => '%s', style => "linespoints", ); $chartCLK->plot2d($dataset); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); LOGDIE("could not find image file: $!\n") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); } ######### Galileo CLK ##################### if ( scalar @{ $file->{'E_CLKs'} } < 1 ) { DEBUG "No Galileo CLKs found"; } else { DEBUG "Plot Galileo CLKs"; $pngName = sprintf ( "%s_CLK_E.png", $station ); my $chartCLK = newChart($station); $chartCLK->set( output => $pngName ); $chartCLK->set( ylabel => "CLK [m]" ); $dataset = Chart::Gnuplot::DataSet->new( xdata => $EPOCHS_E_CLK, ydata => $file->{'E_CLKs'}, title => "Galileo Receiver clock", timefmt => '%s', style => "linespoints", ); $chartCLK->plot2d($dataset); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); LOGDIE("could not find image file: $!\n") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); } ######### Beidou CLK ##################### if ( scalar @{ $file->{'C_CLKs'} } < 1 ) { DEBUG "No BDS CLKs found"; } else { DEBUG "Plot BDS CLKs"; $pngName = sprintf ( "%s_CLK_C.png", $station ); my $chartCLK = newChart($station); $chartCLK->set( output => $pngName ); $chartCLK->set( ylabel => "CLK [m]" ); $dataset = Chart::Gnuplot::DataSet->new( xdata => $EPOCHS_C_CLK, ydata => $file->{'C_CLKs'}, title => "BDS Receiver clock", timefmt => '%s', style => "linespoints", ); $chartCLK->plot2d($dataset); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); LOGDIE("could not find image file: $!\n") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); } ######### ELE ##################### DEBUG "Plot Elevations"; $page = $pdf->page(); $page->mediabox('A4'); $y = $y0 + $dy; $headline = sprintf ( "Satellite Elevations for station %s", $station ); $headline_text = $page->text; $headline_text->font( $font1, 11 / pt ); $headline_text->translate( 15 / mm, 280 / mm ); $headline_text->text($headline); my $chartELE = newChart($station); $chartELE->set( legend => { position => "outside right" } ); # SYSTEM #print Dumper \%ELE; foreach my $key_sys ( sort keys %ELE ) { # print "$key_sys \n"; #print Dumper $ELE{$key_sys}; my @datasets = (); # init datasets my $pngName = sprintf ( "%s_ELE_%s.png", $station, $key_sys ); $chartELE->set( output => $pngName ); $chartELE->set( ylabel => "Elevation [°]", yrange => [ " 0.0 ", " 90.0 " ] ); # SATELLITE foreach my $key_sat ( sort keys %{ $ELE{$key_sys} } ) { # print "$key_sat = $ELE{$key_sys}{$key_sat} \n"; $dataset = Chart::Gnuplot::DataSet->new( xdata => \@{ $ELE{$key_sys}{$key_sat}{EPOCH} }, # array of epochs ydata => \@{ $ELE{$key_sys}{$key_sat}{DATA} }, # array of elevations of one satellite title => "$key_sat", timefmt => '%s', #style => "linespoints", style => "linespoints", ); push ( @datasets, $dataset ); } $chartELE->plot2d(@datasets); # system ("display $pngName&"); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); die ("could not find image file: $!") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); } ######### AMB ##################### DEBUG "Plot Ambiguities"; $page = $pdf->page(); $page->mediabox('A4'); $y = $y0 + $dy; $headline = sprintf ( "Ambiguities for station %s", $station ); $headline_text = $page->text; $headline_text->font( $font1, 11 / pt ); $headline_text->translate( 15 / mm, 280 / mm ); $headline_text->text($headline); # AMBIGUITY_TYPE #print Dumper \%AMB; foreach my $key_ambType (%AMB) { #print "$key_ambType \n"; foreach my $key_sys ( sort keys %{ $AMB{$key_ambType} } ) { #print "$key_sys \n"; print Dumper $AMB{$key_ambType}; my ( @datasets_amb, @datasets_epo ); # init datasets my $pngName_amb = sprintf ( "%s_AMB_%s_%s.png", $station, $key_ambType, $key_sys ); my $pngName_epo = sprintf ( "%s_EPO_%s_%s.png", $station, $key_ambType, $key_sys ); my $chartAMB = Chart::Gnuplot->new( output => $pngName_amb, terminal => 'png', title => $station, ylabel => "Ambiguities $key_ambType [m]", # yrange => [" 0.0 ", " 90.0 "], xlabel => "Time [h]", timeaxis => 'x', xtics => { labelfmt => '%H:%M', rotate => '-270', }, legend => { position => "outside right", }, grid => 'on', ); my $chartEPO = Chart::Gnuplot->new( output => $pngName_epo, terminal => 'png', title => $station, ylabel => "Number of Epochs $key_ambType [-]", # yrange => [" 0.0 ", " 90.0 "], xlabel => "Time [h]", timeaxis => 'x', xtics => { labelfmt => '%H:%M', rotate => '-270', }, legend => { position => "outside right", }, grid => 'on', ); # SATELLITE foreach my $key_sat ( sort keys %{ $AMB{$key_ambType}{$key_sys} } ) { #print "$key_sat = $AMB{$key_ambType}{$key_sys}{$key_sat} \n"; # ambiguities my $dataset_amb = Chart::Gnuplot::DataSet->new( xdata => $AMB{$key_ambType}{$key_sys}{$key_sat}{EPOCH}, # array of epochs ydata => $AMB{$key_ambType}{$key_sys}{$key_sat}{DATA}, # array of ambiguities of one satellite title => "$key_sat", timefmt => '%s', style => "linespoints", #style => "linespoints", ); push ( @datasets_amb, $dataset_amb ); # number of epochs used for ambiguity my $dataset_epo = Chart::Gnuplot::DataSet->new( xdata => $AMB{$key_ambType}{$key_sys}{$key_sat}{EPOCH}, # array of epochs ydata => $AMB{$key_ambType}{$key_sys}{$key_sat}{NUMEPO}, # array of ambiguities of one satellite title => "$key_sat", timefmt => '%s', #style => "linespoints", style => "linespoints", ); push ( @datasets_epo, $dataset_epo ); } # ambiguities $chartAMB->plot2d(@datasets_amb); # system ("display $pngName_amb&"); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); LOGDIE("could not find image file: $!\n") unless -e $pngName_amb; $png->image( $pdf->image_png($pngName_amb), $x, $y, $width, $height ); # number of epochs used for ambiguity $chartEPO->plot2d(@datasets_epo); # system ("display $pngName_epo&"); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); LOGDIE("could not find image file $pngName_epo: $!\n") unless -e $pngName_epo; $png->image( $pdf->image_png($pngName_epo), $x, $y, $width, $height ); } } ######### ION ##################### if ( grep ( $_ eq "ALL", @plotTypes ) ) { DEBUG "Plot ION"; $page = $pdf->page(); $page->mediabox('A4'); $y = $y0 + $dy; $headline = sprintf ( "Ionosphere Delay for station %s", $station ); $headline_text = $page->text; $headline_text->font( $font1, 11 / pt ); $headline_text->translate( 15 / mm, 280 / mm ); $headline_text->text($headline); my $chartION = newChart($station); $chartION->set( ylabel => "Ionophere Delay [m]" ); $chartION->set( legend => { position => "outside right" } ); # SYSTEM foreach my $ksys ( sort keys %ION ) { #print "$key_sys \n"; #print Dumper $ION{$key_sys}; my @datasets; # init datasets my $pngName = sprintf ( "%s_ION_%s.png", $station, $ksys ); $chartION->set( output => $pngName ); # SATELLITE foreach my $sat ( sort keys %{ $ION{$ksys} } ) { #print "$key_sat = $ION{$key_sys}{$key_sat} \n"; my $dataset = Chart::Gnuplot::DataSet->new( xdata => $ION{$ksys}{$sat}{EPOCH}, # array of epochs ydata => $ION{$ksys}{$sat}{DATA}, # array of ionvations of one satellite title => "$sat", timefmt => '%s', #style => " linespoints ", style => "linespoints", ); push ( @datasets, $dataset ); } $chartION->plot2d(@datasets); #system ("display $pngName&"); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); die ("could not find image file: $!") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); } } ######### BIAS ##################### if ( grep ( $_ eq "ALL", @plotTypes ) ) { DEBUG "Plot BIAS"; $page = $pdf->page(); $page->mediabox('A4'); $y = $y0 + $dy; $headline = sprintf ( "Receiver Biases for station %s", $station ); $headline_text = $page->text; $headline_text->font( $font1, 11 / pt ); $headline_text->translate( 15 / mm, 280 / mm ); $headline_text->text($headline); my $chartBIAS = newChart($station); $chartBIAS->set( legend => { position => "outside right" } ); # BIAS_TYPE #print Dumper \%BIA; foreach my $key_biasType ( sort keys %BIA ) { #print "key_biasType: $key_biasType \n"; foreach my $key_sys ( sort keys %{ $BIA{$key_biasType} } ) { #print "key_sys: $key_sys \n"; #print Dumper $BIA{$key_biasType}; my @datasets; # init datasets my $pngName = sprintf ( "%s_BIAS_%s_%s.png", $station, $key_biasType, $key_sys ); $chartBIAS->set( output => $pngName ); $chartBIAS->set( ylabel => "Receiver Bias $key_biasType [m]" ); my $dataset = Chart::Gnuplot::DataSet->new( xdata => $BIA{$key_biasType}{$key_sys}{EPOCH}, ydata => $BIA{$key_biasType}{$key_sys}{DATA}, title => "$key_sys", timefmt => '%s', style => "linespoints", ); push ( @datasets, $dataset ); $chartBIAS->plot2d(@datasets); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); die ("could not find image file: $!") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); } } } ######### RES ##################### DEBUG "Plot Residuals"; $page = $pdf->page(); $page->mediabox('A4'); $y = $y0 + $dy; $headline = sprintf ( "Residuals for station %s", $station ); $headline_text = $page->text; $headline_text->font( $font1, 11 / pt ); $headline_text->translate( 15 / mm, 280 / mm ); $headline_text->text($headline); my $chartRES = newChart($station); $chartRES->set( legend => { position => "outside right" } ); # RESIDUAL_TYPE #print Dumper \%RES; foreach my $key_resType ( sort keys %RES ) { #print "key_resType: $key_resType \n"; #SYSTEM foreach my $key_sys ( sort keys %{ $RES{$key_resType} } ) { #print "key_sys: $key_sys \n"; #print Dumper $RES{$key_resType}; my @datasets; my $pngName = sprintf ( "%s_RES_%s_%s.png", $station, $key_resType, $key_sys ); $chartRES->set( output => $pngName ); $chartRES->set( ylabel => "Residuals $key_resType [m]" ); if ( $key_resType =~ /^c/ ) { $chartRES->set( yrange => [ " -10.0 ", " 10.0 " ] ); } elsif ( $key_resType =~ /^l/ ) { $chartRES->set( yrange => [ " -0.10 ", " 0.10 " ] ); } elsif ( $key_resType =~ /^GIM/ ) { $chartRES->set( yrange => [ " -6.0 ", " 6.0 " ] ); } # SATELLITE foreach my $key_sat ( sort keys %{ $RES{$key_resType}{$key_sys} } ) { #print "$key_sat = $RES{$key_resType}{$key_sys}{$key_sat} \n"; $dataset = Chart::Gnuplot::DataSet->new( xdata => $RES{$key_resType}{$key_sys}{$key_sat}{EPOCH}, # array of epochs ydata => $RES{$key_resType}{$key_sys}{$key_sat}{DATA}, # array of residuals of one satellite title => "$key_sat", timefmt => '%s', #style => "linespoints", style => "linespoints", ); push ( @datasets, $dataset ); } $chartRES->plot2d(@datasets); $y = $y - $dy; if ( $y < 30 / mm ) { $page = $pdf->page(); $page->mediabox('A4'); $y = $y0; } $png = $page->gfx(); LOGDIE("could not find image file: $!\n") unless -e $pngName; $png->image( $pdf->image_png($pngName), $x, $y, $width, $height ); } } } # end if ALL @plotTypes $pdf->save(); $pdf->end(); system ("rm *.png"); if (Common::amInteractiv ) { # system ("evince $inputDir/$pdf_name"); } } # ----- next logfile ----- # newChart returns a default chart. sub newChart { my $title = shift; return Chart::Gnuplot->new( terminal => 'png', title => $title, xlabel => "Time [h]", timeaxis => 'x', xtics => { labelfmt => '%H:%M', rotate => '-270' }, grid => 'on', ); } sub HELP_MESSAGE { print < seconds should be usefull -h, --help show help contents EXAMPLES: $prog --plotTypes ALL -s 30 -l /path/to/FFMJ186730.ppp 2021 andrea.stuerze\@bkg.bund.de EOI_HILFE exit; }