#!/bin/perl -w

package GRID;

# GRID.pm - Module for GRID (WMO code 47) utility routines
#
# Copyright (C) 2002 W. M. Connolley
# wmc@bas.ac.uk / http://www.antarctica.ac.uk/met/wmc
#
# $Author: wmc $:$Date: 2002/07/22 21:01:11 $:$Revision: 1.2 $           
#
# 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.
# http://www.gnu.org/copyleft/gpl.html
#
# Defensive note: this is the first perl I have written using
# modules, objects or multidim arrays. Some of the code was rerewritten
# until I got the syntax right. So beware.

# To-do: something to dump the output in format suitable for reading elsewhere
#      : need to be able to *read* GRIB data and *write* GRID...

require Exporter;
@ISA = qw(Exporter);
@EXPORT = qw();
@EXPORT_OK = ();

use strict;

#
# Define various hashes corresponding to tables (see section D of
# the WMO manual on codes.
# Warning! Most of these tables are incomplete with only enough
# entries for the common codes
#

my %table291 = qw(00 not_present
                  01 pressure
                  02 geopotential_heigth
                  04 temperature
                  23 wind_u
                  24 wind_v
                 );
 
my %table4252 = qw(1 hour    2 day    3 month);

# Hash for looking up F1F2
my %centres = qw(98 ECMWF);

# ----------------------------------------------

sub new {

  my $self = {};
  bless $self;

  $self->{Error_message}="";
  $self->{Error}=0;

  $self->{Warn_message}="";
  $self->{Warn}=0;

  $self->{Skipped}=0;

  $self->{F1F2}="";
  $self->{NNN}="";

  $self->{DEBUG}=0;
 
  $self->{data2}=[[]];
# Ahem...
  $self->{data_max}=-1e30;
  $self->{data_min}=1e30;

  return $self;

};

# ----------------------------------------------

# Simply declare that we are now a GRIB object!
#
# Need to add some translation

sub grid2grib {

  return bless shift,"GRIB"

}

# ----------------------------------------------

sub warn {

  my $self = shift;
  my $message = shift;

  $self->{Warn_message}=$message;
  $self->{Warn}=1;

  return(0)

};

# ----------------------------------------------

sub error {

  my $self = shift;
  my $message = shift;

  $self->{Error_message}=$message;
  $self->{Error}=1;

  return(0)

};

# ----------------------------------------------

sub set_message {

  my $self = shift;
  
  $self->{message} = shift;

};

# ----------------------------------------------

sub check5 {

  my $self = shift;
  my $message = shift;
  my $group;

  while($group = shift) {
    if (length($group) != 5) { 
      return $self->error("[check5] $message group <$group>")
    }
  };

  return 1;

};

# ----------------------------------------------

# At the moment, this just calls encode
#
# Presumably, one day it will do some setup too...

sub as_grid {

  my $self = shift;

  my $GRID = $self->encode();

  return $GRID

};
 
# ----------------------------------------------

# Return the GRID message of a GRID object

sub encode {

  my $self = shift;

  my $rmessage = $self->{message};
  my $message = $$rmessage; $message=~s/\s+\n\s+/\n/g;
  my @message = split('\n',$message);

  my (@groups);				# Groups in a line
  my $extra;				# Extra stuff (which we don't expect)
  my $junk;				# Guess
  my $line="";                  	# One line of message
  my $ident;				# 111, etc

  if ($self->{DEBUG}) { print "starting encode\n" };

  my $GRID = encode_grid($self,$line) or return(0); 
  my $e111 = encode_111($self,shift @message) or return(0);
  my $e333 = encode_333($self,shift @message) or return(0);
# Encode the data lines
  my $rrarray=$self->{data2};
  my $data;
  for (my $ind=0; $ind < $self->{nana}; $ind++) {
    $data .= encode_data($self,$ind,$rrarray->[$ind]) or return(0);
  };
  (my $e555 = $GRID) =~ s/GRID/555/;

#
# All is OK (we think) - return success
#
  if ($self->{DEBUG}) { print "about to return success\n" };

  return "$GRID\n$e111\n$e333\n$data$e555\n777=";

};

sub encode_grid {

my $self = shift;

my $GRID;

#
# Encode GRID line which is
#
# GRID F1F2NNN 1nnntnt
#
# F1F2 = ident of centre (98 for ECMWF)
# NNN  = catalog number of grid (???)
#
# nn   = serial number of this segment
# ntnt = total number of segments (we expect this to be 1)
#
  if ($self->{DEBUG}) { print "starting encode of GRID line\n" };

  $GRID = "GRID $self->{F1F2}$self->{NNN} 1$self->{nn}$self->{ntnt}";

  return $GRID;

};

#
# End of encode_grid
#

sub encode_111 {

  my $self = shift;
  my $e111;

#
# Encode 111 line which is
#
# 111 1a1a1a2a2 (2p1p1p2p2) 3H1H1H1H1 (4....) (5....) 6JJMM 7YYGcGc (8utttt) (9....) 0mmgrgr
#
# a1a1     tens and digits of first var type
# a2a2     ditto, second
#
# p1p1     first ref level in 10's of hPa
# p2p2     second... 99 == "not used"
#
# H1H1H1H1 level (m)
#
# JJ       last 2 digits of year
# MM       month number
# 
# YY       day number (yes, really!)
# GcGc     hour
#
# ut       time unit (table 4252: 1=hour, 2=day, 3=month)
# ttt      time from GcGc to forecast in ut's
#
# mm       model used (table 2677: 10-19=num. anal., 99=not mentioned)
# grgr     grid geometry (table 1487: 99=defined in wmo vol b)
#

  if ($self->{DEBUG}) { print "starting encode of 111 line\n" };

  $e111 = "111 1$self->{a1a1}$self->{a2a2}";

  if (defined $self->{p1p1}) { $e111 .= "2$self->{p1p1}$self->{p2p2}" };

  $e111 .= " 3$self->{h1h1h1h1}";

# Could add 4.... and 5.... optionally sometime

  if (length($self->{yy}) ne 2) { return $self->error("Error: year too long: $self->{yy}") };
  $e111 .= " 6$self->{yy}$self->{mm} 7$self->{dd}$self->{hh}";
 
  if (defined $self->{ut}) { $e111 .= "  8$self->{ut}$self->{ttt}" };

# Also the 9....

  $e111 .= " 0$self->{mm1}$self->{grgr}";

  if ($self->{DEBUG}) { print "  returning line: $e111\n" };

  return $e111;

};

sub encode_333 {

  my $self = shift;

  my $e333;

#
# Encode 333 line which is
#
# 333 1nananpis 2n1n2q1q2 (3usnrr rrrrr) and probably other opts but ignore
#
# nana     number of data lines
# np       grid points per data group
# is       sign, table 1851: 0=sx inc; 1=all +ve, 2=all -ve
#
# n1,n2    number of digits for a parameter
# q1       message contraction, table 3462: 0=spaces between data groups, 2=no
# q2       data contraction, table 3463: 0=all data location groups included
#
# u        scale, table 4200: 0=1, ... 4=10^4, 5=0.1, 6=10^-2, ... 9=10^-5
# sn       sign of ref values (?), table 3845: 0=+ve, 1=-ve
# nrrrrrr  reference value
#
  
  if ($self->{DEBUG}) { print "starting encode of 333 line\n" };

  $e333 = "333 1$self->{nana}$self->{np}$self->{is} 2$self->{n1}$self->{n2}$self->{q1}$self->{q2}";

  if (defined $self->{u}) { $e333 .= " 3$self->{u}$self->{sn}$self->{rr} $self->{rrrrr}" };

  if ($self->{DEBUG}) { print "  returning line: $e333\n" };

  return $e333;

};

sub encode_data {

  my $self = shift;
  my $ind  = shift;		# index number of data line
  my $rdata= shift;		# Ref to data array 

  my ($data,$n);

  if ($self->{DEBUG}) { print "encode_data: $ind... " };

  if ($ind == $self->{pole_row_index}) {
    $data = "98";
    $n = 1;
  } else {
    $data = sprintf("%.2i",$ind+1);
    $n = $self->{ngng}
  };
# $ind*10 is a kludge until I understand how "y" is coded better
  $data .= sprintf("%.2i",$n)." $self->{iaiaia}".sprintf("%.3i",$ind*10)." ";

# Decode data
  my $val;
  for (my $i=0; $i < $n; $i++) {
    $val = ($rdata->[$i])-$self->{rrrrrrr};
    $data .= sprintf("%$self->{n1}.$self->{n1}i",$val)
  }; 

  return $data."\n";

};

# ----------------------------------------------

sub decode {

  my $self = shift;

  my $rmessage = $self->{message};
  my $message = $$rmessage; $message=~s/\s+\n\s+/\n/g;
  my @message = split('\n',$message);

  my (@groups);				# Groups in a line
  my $extra;				# Extra stuff (which we don't expect)
  my $junk;				# Guess
  my $line="";                  	# One line of message
  my $ident;				# 111, etc

  if ($self->{DEBUG}) { print "starting decode\n" };

  while (defined ($line=shift @message)) {
    if ($line !~ /^GRID/) { $self->{Skipped}++ } else { last };
    if ($line =~ /GRID/) { $self->warn("found GRID in line: <$line>") }
  };

  if (!$line) {
    return $self->error("No GRID found in $self->{Skipped} lines")
  };

# Decode the GRID line
  decode_grid($self,$line) or return(0); 
#        the 111
  decode_111($self,shift @message) or return(0);
#        the 333
  decode_333($self,shift @message) or return(0);
#        the data lines
  my $rrarray=$self->{data2};
  for (my $ind=0; $ind < $self->{nana}; $ind++) {
    $rrarray->[$ind]=[];
    decode_data($self,shift @message,$ind,$rrarray->[$ind]) or return(0);
  };
  if (defined $self->{pole_row_index}) { decode_fill_pole($self,$rrarray->[$self->{pole_row_index}]) };

# We could go on and do the redundant 555 group and the 666/777 terminator
# but for now I see no reason too

#
# All is OK (we think) - return success
#
  if ($self->{DEBUG}) { print "about to return success\n" };
  return 1;

};

sub decode_grid {

my $self = shift;
my $line = shift;

#
# Decode GRID line which is
#
# GRID F1F1NNN 1nnntnt
#
# F1F2 = ident of centre (98 for ECMWF)
# NNN  = catalog number of grid (???)
#
# nn   = serial number of this segment
# ntnt = total number of segments (we expect this to be 1)
#
  if ($self->{DEBUG}) { print "starting decode of GRID line <$line>\n" };
  my ($junk,$g1,$g2,$extra) = split(/\s+/,$line);

  if ($extra) { 
    return $self->error("on GRID line: Wasn't expecting extra: $extra")
  };

  if (length($g1) != 5) { return $self->error("F1F2NNN group ($g1) length wrong") };

  ($self->{F1F2},$self->{NNN}) = ($g1=~/(..)(...)/);
  if (!($self->{centre}=$centres{$self->{F1F2}})) { $self->{centre}="Unknown: $self->{F1F2}" };

  if (length($g2) != 5) { return $self->error("nnntnt group length wrong") };

  ($self->{nn},$self->{ntnt}) = ($g2=~/.(..)(..)/);
  if ($self->{ntnt} > 1 or $self->{nn} > 1) {
    return $self->error("This message is split into $self->{ntnt} segments (this is $self->{nn}) which I don't understand")
  };

  return 1;

};

#
# End of decode_grid
#

sub decode_111 {

  my $self = shift;
  my $line = shift;

#
# Decode 111 line which is
#
# 111 1a1a1a2a2 (2p1p1p2p2) 3H1H1H1H1 (4....) (5....) 6JJMM 7YYGcGc (8utttt) (9....) 0mmgrgr
#
# a1a1     tens and digits of first var type
# a2a2     ditto, second
#
# p1p1     first ref level in 10's of hPa
# p2p2     second... 99 == "not used"
#
# H1H1H1H1 level (m)
#
# JJ       last 2 digits of year
# MM       month number
# 
# YY       day number (yes, really!)
# GcGc     hour
#
# ut       time unit (table 4252: 1=hour, 2=day, 3=month)
# ttt      time from GcGc to forecast in ut's
#
# mm       model used (table 2677: 10-19=num. anal., 99=not mentioned)
# grgr     grid geometry (table 1487: 99=defined in wmo vol b)
#

  if ($self->{DEBUG}) { print "starting decode of 111 line <$line>\n" };
  my ($ident,@groups) = split(/\s+/,$line);
  my $junk;

  $self->check5("111 line",@groups) or return (0);

  my $g1 = shift @groups;
  ($self->{a1a1},$self->{a2a2}) = ($g1=~/.(..)(..)/);

  $g1 = shift @groups;
  ($junk,$self->{p1p1},$self->{p2p2}) = ($g1=~/(.)(..)(..)/);
  if ($junk != "2") {
    undef $self->{p1p1};
    unshift @groups, $g1
  };

  $g1 = shift @groups;
  ($junk,$self->{h1h1h1h1}) = ($g1=~/(.)(....)/);
  if ($junk != "3") {
    undef $self->{h1h1h1h1};
    unshift @groups, $g1
  };

  $g1 = shift @groups;
  ($junk,$self->{yy},$self->{mm}) = ($g1=~/(.)(..)(..)/);
  if ($junk != "6") {
    return $self->error("111 line wanted 6.... group but got <$g1>")
  };

  $g1 = shift @groups;
  ($junk,$self->{dd},$self->{hh}) = ($g1=~/(.)(..)(..)/);
  if ($junk != "7") {
   return $self->error("111 line wanted 7.... group but got <$g1>")
  };

  $g1 = shift @groups;
  ($junk,$self->{ut},$self->{ttt}) = ($g1=~/(.)(.)(...)/);
  if ($junk != "8") {
    undef $self->{ut};
    unshift @groups, $g1
  };

  $g1 = shift @groups;
  ($junk,$self->{mm1},$self->{grgr}) = ($g1=~/(.)(..)(..)/);
  if ($junk != "0") {
    return $self->error("111 line wanted 0.... group but got <$g1>")
  };

  return 1;
};

sub decode_333 {

  my $self = shift;
  my $line = shift;

#
# Decode 333 line which is
#
# 333 1nananpis 2n1n2q1q2 (3usnrr rrrrr) and probably other opts but ignore
#
# nana     number of data lines
# np       grid points per data group
# is       sign, table 1851: 0=sx inc; 1=all +ve, 2=all -ve
#
# n1,n2    number of digits for a parameter
# q1       message contraction, table 3462: 0=spaces between data groups, 2=no
# q2       data contraction, table 3463: 0=all data location groups included
#
# u        scale, table 4200: 0=1, ... 4=10^4, 5=0.1, 6=10^-2, ... 9=10^-5
# sn       sign of ref values (?), table 3845: 0=+ve, 1=-ve
# nrrrrrr  reference value
#
  
  my ($ident,@groups) = split(/\s+/,$line);
  my ($junk,$g1,$g2);
 
  if ($self->{DEBUG}) { print "starting decode of 333 line <$line>\n" };

  $self->check5("333 line",@groups) or return (0);

  $g1 = shift @groups;
  ($junk,$self->{nana},$self->{np},$self->{is}) = ($g1=~/(.)(..)(.)(.)/);

  $g1 = shift @groups;
  ($junk,$self->{n1},$self->{n2},$self->{q1},$self->{q2}) = ($g1=~/(.)(.)(.)(.)(.)/);
  if ($junk != "2") {
    return $self->error("333 line wanted 2.... group but got <$g1>")
  };

# nb - these 2 groups (3.... .....) may not be present
  $g1 = shift @groups; 
  $g2 = shift @groups;
  if ($g1) {
    ($junk,$self->{u},$self->{sn},$self->{rr}) = ($g1=~/(.)(.)(.)(..)/)
  };
  if ($g2) {
    ($self->{rrrrr}) = ($g2=~/(.....)/)
  };
  if (!$g1 or $junk != "3") {
# In this case, no ref value, so set to zero
    $self->{rrrrrrr}=0;
  } else {
# Set from the 3.... group
    $self->{rrrrrrr}=$self->{rr}.$self->{rrrrr};
    if ($self->{sn} == 1) { $self->{rrrrrrr}=-$self->{rrrrrrr} };
  };

  return 1;
};

sub decode_data {

  my $self = shift;
  my $line = shift;		# data line
  my $ind  = shift;		# index number of data line
  my $rdata= shift;		# Ref to data array to fill

  my ($k1ng, $iaj2, $k1k1, $ngng, $iaiaia, $j2j2j2, $data);

# Split line into k1ng (line index and npts line)
#                 iaj2 (lon and lat coords)
#                 data (data)
  ($k1ng,$iaj2,$data) = split /\s+/, $line;

# Split k1ng into k1k1 (index number)
#                 ngng (npts in line)
  ($k1k1,$ngng) = ($k1ng=~/(..)(..)/);
  if ($k1k1 eq "98") {
    $self->{pole_row_index} = $ind
  };
  if ($k1k1 != $ind+1 and $k1k1 ne "98") { 
    return $self->error("data line $ind thinks its $k1k1 <$line>") 
  };
# Define ngng (number of pts value) for use at pole later
  if (!defined $self->{ngng} and $k1k1 ne "98") { $self->{ngng}=$ngng };

# Split i1j2 into iaiaia (longitude of first point in line)
#                 j2j2j2 (latitude, in something!)
  ($self->{iaiaia},$self->{j2j2j2}) = ($iaj2=~/(...)(...)/);

# Decode data
  for (my $i=0; $i<$ngng; $i++) {
    my $datapt=substr($data,$i*$self->{n1},$self->{n1})+$self->{rrrrrrr};
    $rdata->[$i]=$datapt;
    if ($datapt>$self->{data_max}) { $self->{data_max}=$datapt };
    if ($datapt<$self->{data_min}) { $self->{data_min}=$datapt };
  }; 
  return 1;

};

# ----------------------------------------------

sub decode_fill_pole {

  my $self = shift;
  my $rdata= shift;		# Ref to data array to fill

  if ($self->{DEBUG}) { print "  filling pole row with $rdata->[0]\n" };

  for (my $i=0; $i < $self->{ngng}; $i++) {
    $rdata->[$i] = $rdata->[0]
  }

}

# ----------------------------------------------

# Convert data from the GRID-native type 2d array format to a long 1-d
# array compatible with GRIB
#
# Needs decode_fill_pole called first

sub data_2d_to_1d {

  my $self = shift;

  $self->{data} = [];

  my $ind = 0;
  for (my $i=0; $i < $self->{nana}; $i++) {
    for (my $j=0; $j < $self->{ngng}; $j++) {
      $self->{data}->[$ind++] = $self->{data2}->[$i]->[$j]
    }
  }

};

# ----------------------------------------------

#
# sub to return a string reprenting the state of the GRID object.
#
# use: print $grid->as_text("opt2"=>"value1", "opt2"=>"value2", ...)
#
# options:
#   "data" => "range"		just print the range of the data
#   "header" => "one-line"	one-line summary of the record

sub as_text {

  my $self = shift;
  my %option = (
                 header => "one-line"
               , data => "range"
               , @_
               );

  my ($forecast, $v1, $v2, $method);

#
# First off, print any warnings that the decode generated
#

  if ($self->{Warn}) { print $self->{Warn_message} };

#
# Parameter type - eg pressure
#
  if (!($v1=$table291{$self->{a1a1}})) { $v1="unknown ($self->{a1a1})" };
  if (!($v2=$table291{$self->{a2a2}})) { $v2="unknown ($self->{a2a2})" };

#
# Maybe a forecast
#
  if ($self->{ut}) {
    $forecast="  Forecast valid at +$self->{ttt} $table4252{$self->{ut}}"
  } else {
    $forecast=""
  };

#
# Model type
#
  if ($self->{mm1} == 99) {
    $method="NWP"
  } elsif (int($self->{mm1}/10) == 1) {
    $method="numerical analysis"
  } else {
    $method="unknown"
  };

#
# Reference level(s)
#

  my $reflevs="";
  if ($self->{p1p1}) {
    $reflevs=$reflevs.($self->{p1p1}*10);
    if ($self->{p2p2} != 99) {$reflevs.=" and ".($self->{p2p2}*10) }
  } elsif ($self->{h1h1h1h1}) {
    my $rl=$self->{h1h1h1h1};
    if ($rl == "0000") { $rl="MSL" };
    $reflevs=$reflevs.$rl
  } else {
    $reflevs=""
  };

#
# Now construct return string
#

  my $rets = 
  "GRID begins after $self->{Skipped} lines\n"
  .
  "  Centre: $self->{centre} and grid type $self->{NNN};  Segment number: $self->{nn} of $self->{ntnt}\n"
  .
  "  Date  : $self->{yy}/$self->{mm}/$self->{dd}:$self->{hh} $forecast\n"
  .
  "  Var1,2: $v1, $v2;    Reflevs: $reflevs\n"
  .
  "  Method: $method;       Grid geometry: $self->{grgr}\n" 
  .
  "  Refval: $self->{rrrrrrr}\n"
  .
  "  Grid size: $self->{nana} by $self->{ngng}\n"
  ;

  if ($option{header} eq "one-line") {
    $rets="$self->{centre} ".
          "$self->{yy}/$self->{mm}/$self->{dd}:$self->{hh}".
          $forecast . 
          " $v1 $reflevs" .
          "\n"
  };

  if ($option{data} eq "range") {
    $rets .= "Data range: [$self->{data_min},$self->{data_max}]\n"
  } else {
    $rets .= "\nData\n----\n\n";
    for (my $i=0; $i<$self->{nana}; $i++) {
      my $ra = $self->{data2}->[$i];
      my $prets = "line $i: ";
      for (my $j=0; $j<0+@$ra; $j++) {
        $prets .= $ra->[$j]." ";
      };
      $rets .= $prets."\n";
    };
  };

  return $rets;

};

# ----------------------------------------------

1;