#!/usr/bin/perl 

####################################################################
#British Crown Copyright (c) 2009, the Met Office All rights reserved.
#
#Redistribution and use in source and binary forms, with or without 
#modification, are permitted provided that the following conditions 
#are met:
#
#Redistributions of source code must retain the above copyright notice, 
#this list of conditions and the following disclaimer. 
#
#Redistributions in binary form must reproduce the above copyright 
#notice, this list of conditions and the following disclaimer in the 
#documentation and/or other materials provided with the distribution. 
#
#Neither the name of the Met Office nor the names of its contributors 
#may be used to endorse or promote products derived from this software 
#without specific prior written permission. 
#
#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
#"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
#LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 
#A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 
#HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 
#SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
#LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
#DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
#THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
#(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
#OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#########################################################


# Gridder for Land Surface data

use strict;
use warnings;

# Path to directory containing the station files
# Change this to wherever you've unpacked the ZIP file.
my $Base_dir = "station_files";

# Accumulate the station anomalies in the structure of the grid
#  (monthly, 5x5 degrees)
# $GridAnom[189803][0][0] will be an array containing all the anomalies
# from March 1898 for longitude -180 to -175 and latitude -90 to -85
# $GridAnom[196909][27][15] Is September 1969, long -45 to -40, Lat -15 to -10
my $GridAnom = populate_grid($Base_dir);

# For each month, print out the gridded data - mean anomaly for the
#  grid cells with data, -1.000e+30 otherwise.
for ( my $Year = 1850 ; $Year <= 2009 ; $Year++ ) {
    for ( my $Month = 1 ; $Month <= 12 ; $Month++ ) {
        my $key = sprintf "%04d%02d", $Year, $Month;
        printf
          "  %04d    %2d     1    36 rows     72 columns. Missing=-1.000e+30\n",
          $Year, $Month;
        for ( my $Row = 35 ; $Row >= 0 ; $Row-- ) {
            for ( my $Column = 0 ; $Column < 72 ; $Column++ ) {
                if ( defined( $GridAnom->[$key][$Column][$Row] ) )
                {    # There are stations in this grid cell
                    printf "%10.3e ",
                      mean( @{ $GridAnom->[$key][$Column][$Row] } );
                }
                else {    # No stations in this grid cell
                    printf "%10.3e ", -1.0e+30;   # Arbitrary missing data value
                }
            }
            print "\n";
        }
    }
}

# Populate @GridA from the station data files
sub populate_grid {

    my $Data_source = shift;
    my @GridA;

    opendir( DIRH, "$Data_source" )
      or die "Can't open Base directory $Data_source";
    my @Directories = readdir(DIRH);
    closedir(DIRH);

    foreach my $Dir (@Directories) {

        # Include all directories except '.' and '..'
        if ( -d "$Data_source/$Dir" && $Dir !~ /^\./ ) {
            opendir( WMOREGION, "$Data_source/$Dir" )
              or die "Couldn't open data directory $Dir";
            my @Files = readdir(WMOREGION);
            close(WMOREGION);
            foreach (@Files) {

                #read only regular files, not directories
                if ( !-d "$Data_source/$Dir/$_" ) {

                    open( DIN, "$Data_source/$Dir/$_" )
                      or die "Couldn't open station file $_";
                    my $Station = read_station( \*DIN )
                      or die "Failed reading station file $_";
                    close(DIN);

                    # Ignore stations without the necessary metadata
                    if (   defined( $Station->{long} )
                        && defined( $Station->{lat} )
                        && defined( $Station->{start_year} )
                        && defined( $Station->{end_year} ) )
                    {

                        # Latitude and longitude indices for @GridA
                        my $LonI =
                          int( ( -$Station->{long} + 180 ) / 5 )
                          ;    # Station files have lon in degW.
                        my $LatI = int( ( $Station->{lat} + 90 ) / 5 );

                        # Push anomaly for each year and month to @GridA
                        for (
                            my $i = $Station->{start_year} ;
                            $i <= $Station->{end_year} ;
                            $i++
                          )
                        {
                            for ( my $j = 1 ; $j <= 12 ; $j++ ) {
                                my $key = sprintf "%04d%02d", $i, $j;
                                if ( defined( $Station->{anomalies}->{$key} ) )
                                {
                                    push @{ $GridA[$key][$LonI][$LatI] },
                                      $Station->{anomalies}->{$key};
                                }
                            }
                        }
                    }

                }

            }
        }
    }
    return \@GridA;
}

# Read in a station record from a station file
sub read_station {

    # File handle to read from is passed as argument
    my $fh = shift;

    # Define a data structure to hold the station data
    my %data = (
        header          => undef,
        number          => undef,
        name            => undef,
        country         => undef,
        start_year      => undef,    # first year data available
        end_year        => undef,    # last year data available
        first_good_year => undef,
        lat             => undef,
        long            => undef,
        normals         => undef,
        sds             => undef,
        temperatures    => undef,    # keys are dates (YYYY/MM)
        anomalies       => undef,    # keys are dates (YYYY/MM)
    );

    # Populate the data structure from the file
    while (<$fh>) {
        if ( $_ =~ /^Obs:/ ) { last; }    # move on to reading obs
        if ( $_ =~ /^Number= \s*(\d+)/ ) { $data{number} = $1; next; }
        if ( $_ =~ /^Name= (.+)$/ )          { $data{name}       = $1; next; }
        if ( $_ =~ /^Country= (.+)$/ )       { $data{country}    = $1; next; }
        if ( $_ =~ /^Lat= \s*(\S+)/ )        { $data{lat}        = $1; next; }
        if ( $_ =~ /^Long= \s*(\S+)/ )       { $data{long}       = $1; next; }
        if ( $_ =~ /^Start year= \s*(\d+)/ ) { $data{start_year} = $1; next; }
        if ( $_ =~ /^End year= \s*(\d+)/ )   { $data{end_year}   = $1; next; }
        if ( $_ =~ /^First Good year= \s*(\d+)/ ) {
            $data{first_good_year} = $1;
            next;
        }
        if ( $_ =~ /^Normals=/ ) {

            #read in 12 monthly normals
            my @fields = split;
            for ( my $i = 1 ; $i <= 12 ; $i++ ) {
                $data{normals}[$i] = $fields[$i];
                if ( $data{normals}[$i] < -90 ) {
                    undef( $data{normals}[$i] );
                }
            }
            next;
        }
        if ( $_ =~ /^Standard deviations=/ ) {

            #read in 12 monthly standard deviations
            my @fields = split;
            for ( my $i = 1 ; $i <= 12 ; $i++ ) {
                $data{sds}[$i] = $fields[ $i + 1 ];
                if ( $data{sds}[$i] < -90 ) {
                    undef( $data{sds}[$i] );
                }
            }
            next;
        }

    }

    # Read in all the obs and make the anomalies
    while (<$fh>) {

        #read in 12 monthly observations for each year
        my @fields = split;
        for ( my $i = 1 ; $i <= 12 ; $i++ ) {
            if ( $fields[$i] < -90 ) { next; }
            my $key = sprintf "%04d%02d", $fields[0], $i;
            $data{temperatures}{$key} = $fields[$i];

        # Round anomalies to nearest 0.1C - but skip them if too far from normal
            if (   defined( $data{normals}[$i] )
                && $data{normals}[$i] > -90
                && defined( $data{sds}[$i] )
                && $data{sds}[$i] > -90
                && abs( $data{temperatures}{$key} - $data{normals}[$i] ) <=
                ( $data{sds}[$i] * 5 ) )
            {
                $data{anomalies}{$key} = sprintf "%5.1f",
                  $data{temperatures}{$key} - $data{normals}[$i];
            }
        }
    }
    return \%data;    # Return reference to structure
}

sub mean {
    if ( scalar(@_) == 0 ) { return; }
    my $Total = 0;
    foreach (@_) { $Total += $_; }
    $Total /= scalar(@_);
    return $Total;
}