Differences

This shows you the differences between two versions of the page.

--- bufr.pm:mergebufrtemp.pl [2015-02-05 07:58:08]
pals
+++ bufr.pm:mergebufrtemp.pl [2022-05-31 09:29:31]
@@ Line 1: / Line 1: @@
-<code perl>
-#!/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') {
-    splice(@$desc, $ii, 10, @pressure_level_desc);
-    if ($desc1->[$i1] eq '031001') {
-        # No more pressure levels in message 1, so fetch from message 2
-        splice(@$data, $ii, 10, @$data2[$i2 .. $i2+9]);
-        $i2 += 10;
-    } elsif ($desc2->[$i2] eq '031001') {
-        # No more pressure levels in message 2, so fetch from message 1
-        splice(@$data, $ii, 10, @$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
-        splice(@$data, $ii, 10, @$data1[$i1 .. $i1+9]);
-        $i1 += 10;
-    } else {
-        # pick the highest pressure, i.e. data2
-        splice(@$data, $ii, 10, @$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) {
-        splice(@$desc, $ii, 7, @windshear_level_desc);
-        if ($nws1 == 0) {
-            # Fetch from message 2
-            splice(@$data, $ii, 7, @$data2[$i2 .. $i2+6]);
-            $i2 += 7;
-            $nws2--;
-        } elsif ($nws2 == 0) {
-            # Fetch from message 1
-            splice(@$data, $ii, 7, @$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
-            splice(@$data, $ii, 7, @$data1[$i1 .. $i1+6]);
-            $i1 += 7;
-            $i2 += 7;
-            $nws2--;
-        } elsif ($data1->[$i1+2] > $data2->[$i2+2]) {
-            # pick the highest pressure, i.e. data1
-            splice(@$data, $ii, 7, @$data1[$i1 .. $i1+6]);
-            $i1 += 7;
-            $nws1--;
-        } else {
-            # pick the highest pressure, i.e. data2
-            splice(@$data, $ii, 7, @$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) = @_;
-    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;
-}
-=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<--he> or C<-he> 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
-</code>