#!/usr/bin/perl -w # (C) Copyright 2015, met.no # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA # 02110-1301, USA. # pod included at end of file use strict; use Getopt::Long; use Pod::Usage qw(pod2usage); use File::Slurp qw(write_file); use Geo::BUFR; # Will be used if neither --tablepath nor $ENV{BUFR_TABLES} is set use constant DEFAULT_TABLE_PATH => '/usr/local/lib/bufrtables'; # Parse command line options our %option = (); GetOptions( \%option, 'help', 'outfile=s', 'tablepath=s', ) or pod2usage(-verbose => 0); # User asked for help pod2usage(-verbose => 1) if $option{help}; # Make sure there are 2 input files pod2usage(-verbose => 0) unless @ARGV == 2; my $infile1 = $ARGV[0]; my $infile2 = $ARGV[1]; # Set BUFR table path if ($option{tablepath}) { # Command line option --tablepath overrides all Geo::BUFR->set_tablepath($option{tablepath}); } elsif ($ENV{BUFR_TABLES}) { # If no --tablepath option, use the BUFR_TABLES environment variable Geo::BUFR->set_tablepath($ENV{BUFR_TABLES}); } else { # If all else fails, use the libbufr bufrtables Geo::BUFR->set_tablepath(DEFAULT_TABLE_PATH); } my $bufr1 = Geo::BUFR->new(); $bufr1->fopen($infile1); my $bufr2 = Geo::BUFR->new(); $bufr2->fopen($infile2); my ($data1, $desc1) = $bufr1->next_observation(); my ($data2, $desc2) = $bufr2->next_observation(); die "Unexpanded descriptors in section 3 for $infile1 must be 309052, is " . $bufr1->get_descriptors_unexpanded() if $bufr1->get_descriptors_unexpanded() ne '309052'; die "Unexpanded descriptors in section 3 for $infile2 must be 309052, is " . $bufr2->get_descriptors_unexpanded() if $bufr2->get_descriptors_unexpanded() ne '309052'; die "Can only handle single subset BUFR messages, but $infile1 contains " . $bufr1->get_number_of_subsets() . " subsets" if $bufr1->get_number_of_subsets() != 1; die "Can only handle single subset BUFR messages, but $infile2 contains " . $bufr2->get_number_of_subsets() . " subsets" if $bufr2->get_number_of_subsets() != 1; my $bufr3 = Geo::BUFR->new(); $bufr3->copy_from($bufr1,'metadata'); my @pressure_level_desc = ( '004086', # Long time period or displacement [s] '008042', # Extended vertical sounding significance [flag table] '007004', # Pressure [PA] '010009', # Geopotential height [GPM] '005015', # Latitude displacement (high accuracy) [degree] '006015', # Longitude displacement (high accuracy) [degree] '012101', # Temperature/dry-bulb temperature [K] '012103', # Dew-point temperature [K] '011001', # Wind direction [degree true] '011002', # Wind speed [m/s] ); my @windshear_level_desc = ( '004086', # Long time period or displacement [s] '008042', # Extended vertical sounding significance [flag table] '007004', # Pressure [PA] '005015', # Latitude displacement (high accuracy) [degree] '006015', # Longitude displacement (high accuracy) [degree] '011061', # Absolute wind shear in 1 km layer below [m/s] '011062', # Absolute wind shear in 1 km layer above [m/s] ); my ($data, $desc); # First handle section 4 up to where pressure level data starts for my $ii (0 .. 27) { $desc->[$ii] = $desc1->[$ii]; if (defined($data1->[$ii])) { $data->[$ii] = $data1->[$ii]; } else { $data->[$ii] = $data2->[$ii]; } } # Then handle pressure levels (and eddrf) $desc->[28] = '031002'; $data->[28] = undef; # Will be filled with $eddrf later my $eddrf = 0; # extended delayed descriptor replication factor my ($ii,$i1,$i2) = (29,29,29); while ($desc1->[$i1] ne '031001' || $desc2->[$i2] ne '031001') { push(@$desc, @pressure_level_desc); if ($desc1->[$i1] eq '031001') { # No more pressure levels in message 1, so fetch from message 2 push(@$data, @$data2[$i2 .. $i2+9]); $i2 += 10; } elsif ($desc2->[$i2] eq '031001') { # No more pressure levels in message 2, so fetch from message 1 push (@$data, @$data1[$i1 .. $i1+9]); $i1 += 10; } elsif ($data1->[$i1+2] eq $data2->[$i2+2]) { # Same pressure level, should be merged into one # evss = extended vertical sounding significance my $evss1 = $data1->[$i1+1]; my $evss2 = $data2->[$i2+1]; my $evss = combine_evss($evss1,$evss2); push(@$data, defined($data1->[$i1]) ? $data1->[$i1] : $data2->[$i2]); push(@$data, $evss); for (2 .. 9) { push(@$data, defined($data1->[$i1 + $_]) ? $data1->[$i1 + $_] : $data2->[$i2 + $_]); } $i1 += 10; $i2 += 10; } elsif ($data1->[$i1+2] > $data2->[$i2+2]) { # pick the highest pressure, i.e. data1 push(@$data, @$data1[$i1 .. $i1+9]); $i1 += 10; } else { # pick the highest pressure, i.e. data2 push(@$data, @$data2[$i2 .. $i2+9]); $i2 += 10; } $eddrf++; $ii += 10; } $data->[28] = $eddrf; # Handle windshear levels my $iws = $ii; $desc->[$iws] = '031001'; $data->[$iws] = undef; # Will be filled later if ($data1->[$i1] == 0 && $data2->[$i2] == 0) { $data->[$iws] = 0; # No wind shear levels } else { my $nws = 0; # Number of wind shear levels, to be calculated my $nws1 = $data1->[$i1++]; my $nws2 = $data2->[$i2++]; $data->[$ii++] = 0; while ($nws1 > 0 || $nws2 > 0) { push(@$desc, @windshear_level_desc); if ($nws1 == 0) { # Fetch from message 2 push(@$data, @$data2[$i2 .. $i2+6]); $i2 += 7; $nws2--; } elsif ($nws2 == 0) { # Fetch from message 1 push(@$data, @$data1[$i1 .. $i1+6]); $i1 += 7; $nws1--; } elsif ($data1->[$i1+2] eq $data2->[$i2+2]) { # see no reason why data1 should differ from data2 here, # so we use data1 push(@$data, @$data1[$i1 .. $i1+6]); $i1 += 7; $i2 += 7; $nws2--; } elsif ($data1->[$i1+2] > $data2->[$i2+2]) { # pick the highest pressure, i.e. data1 push(@$data, @$data1[$i1 .. $i1+6]); $i1 += 7; $nws1--; } else { # pick the highest pressure, i.e. data2 push(@$data, @$data2[$i2 .. $i2+6]); $i2 += 7; $nws2--; } $nws++; $ii += 7; } $data->[$iws] = $nws; } my ($data_refs, $desc_refs); # One subset only $data_refs->[1] = $data; $desc_refs->[1] = $desc; # Print the encoded BUFR message my $msg = $bufr3->encode_message($data_refs, $desc_refs); if ($option{outfile}) { write_file($option{outfile}, {binmode => ':raw'}, $msg); } else { print $msg; } ## Combine two extended vertical sounding significance (evss) into one sub combine_evss { my ($evss1,$evss2) = @_; if (defined($evss1) && defined($evss2)) { my $binary1 = pack "N", $evss1; # Packed as 32 bits in big-endian order my $binary2 = pack "N", $evss2; my $binary = $binary1 | $binary2; # bitwise or return unpack "N", $binary; } elsif (defined($evss1)) { return $evss1; } else { return $evss2; } } =pod =head1 SYNOPSIS mergebufrtemp.pl <bufr file 1> <bufr file 2> [--outfile <file>] [--tablepath <path to BUFR tables>] [--help] =head1 DESCRIPTION Will merge (the first) BUFR message in <bufr file 1> with (the first) BUFR message in <bufr file 2> into one single BUFR message, printing to STDOUT unless C<--outfile> is set. The 2 input BUFR messages are assumed to be single subset BUFR TEMPs utilizing TM305092 (this is checked) containing data for the same station and termin (this is not checked). Metadata are fetched from input file 1, the same applies for conflicting data in section 4 (e.g. if BUFR launch time differs) except that input file 2 is used for data missing in input file 1. Execute without arguments for Usage, with option C<--help> for some additional info. =head1 OPTIONS --outfile <filename> Will print to <filename> instead of STDOUT --tablepath <path to BUFR tables> Set path to BUFR tables (overrides $ENV{BUFR_TABLES}) --help Display Usage and explain the options used. Almost the same as consulting perldoc mergebufrtemp.pl Options may be abbreviated, e.g. C<--h> or C<-h> for C<--help>. To avoid having to use the C<--tablepath> option, you are adviced to set the environment variable BUFR_TABLES to the directory where your BUFR tables are located (unless the default path provided by mergebufrtemp.pl works for you). =head1 AUTHOR Pål Sannes E<lt>pal.sannes@met.noE<gt> =head1 COPYRIGHT Copyright (C) 2015 met.no =cut