#!/usr/bin/perl -w

#  Program to extract daily Min/Max Temperature and Precipitation
#  from GHCN data set and fill missing values
#
#  Willem A. Schreuder
#  August 17, 2004
#  Adapted from TD3200 to GHCN Oct 3, 2012

use strict;
use lib '/pm/grflib/perl';
use Calendar;
use DBF;

#  Where to end
my $jd1 = (@ARGV>0) ? JulianDay($ARGV[0],12,31) : 0;
 
#  GHCN parameters to extract
#   PRCP = Precipitation (tenths of mm)
#   TMAX = Maximum temperature (tenths of degrees C)
#   TMIN = Minimum temperature (tenths of degrees C)
my %par  = (TMIN=>'TC', TMAX=>'TC' , PRCP=>'TM');

#  Alternate stations to fill data
my %alt  = (
   USC00050109=>['USC00059295','USC00057950'],   #  Akron 4E > Yuma > Sterling
   USW00094040=>['USC00252065','USC00251415'],   #  McCook > Culbertson > Cambridge
   USC00257070=>['USC00253037']);                #  RedCloud > Franklin #2
my @stn = keys %alt;
my %name = (
   USC00050109=>'050109',
   USW00094040=>'255310',
   USC00257070=>'257070',
   );

#  Days to fill by interpolation
my %fill = (PRCP=>2 , TMIN=>15 , TMAX=>15);

#  Stations to extract from database
my %load;
foreach my $stn (@stn)
{
   $load{$stn}++;
   foreach my $alt (@{$alt{$stn}})
   {
      $load{$alt}++;
   }
}

#
#  Read and process GHCN input data file
#
my %ttp;
open(TD , "<GHCN") || die "Cannot open data file GHCN\n";
while (my $line = <TD>)
{
   chomp $line;
   my ($stn,$date,$elem,$val) = split(/,/,$line);
   #  Do only relevant stations
   exists($load{$stn}) || next;
   #  Limit fields
   exists($par{$elem}) || next;
   #  Get date
   my $yr = substr($date,0,4);
   my $mo = substr($date,4,2);
   my $dy = substr($date,6,2);
   my $jd = JulianDay($yr,$mo,$dy);
   ($jd>$jd1) && ($jd1 = $jd);
   #  Extract daily data and convert units
   #  1/10 mm to in
   if ($par{$elem} eq 'TM')
   {
      $val = sprintf '%.2f' , $val/254;
   }
   #  1/10 C to F
   else
   {
      $val = sprintf '%.1f' , 9/50*$val+32;
   }
   #  Check for duplicates
   (exists($ttp{$stn}{$jd}{$elem}) && ($ttp{$stn}{$jd}{$elem} != $val)) &&
         warn "Duplicate $stn $mo/$dy/$yr $jd $elem: $ttp{$stn}{$jd}{$elem} $val\n";
   #  Store
   $ttp{$stn}{$jd}{$elem} = $val;
}
close(TD);

#
#  Fill in gaps in all parameters
#
#  ET Stations  (twice)
foreach my $stn (@stn,@stn)
{
   #  Fill with alterantives if station is completely missing
   my ($jd0,@jd) = sort {$a <=> $b} keys %{$ttp{$stn}};
   $jd0 || next;
   #  Tmin, Tmax and PRCP
   foreach my $item (sort keys %fill)
   {
      my ($jd0,@jd) = sort {$a <=> $b} keys %{$ttp{$stn}};
      defined($jd0) || die "No data $stn\n";
      ($jd[-1]<$jd1) && (push @jd, $jd1);
      #  Initial gap
      while (!exists($ttp{$stn}{$jd0}{$item}))
      {
         foreach my $alt (@{$alt{$stn}})
         {
            exists($ttp{$alt}{$jd0}{$item}) && ($ttp{$stn}{$jd0}{$item} = $ttp{$alt}{$jd0}{$item});
         }
         $jd0 = shift @jd;
      }
      #  Internal gap
      foreach my $jd (@jd)
      {
         exists($ttp{$stn}{$jd}{$item}) || next;
         my $n = $jd-$jd0-1;
         #  Gap too large - fill from alternates
         if ($n>$fill{$item})
         {
            for (my $j=$jd0+1;$j<$jd;$j++)
            {
               foreach my $alt (@{$alt{$stn}})
               {
                  exists($ttp{$alt}{$j}{$item}) && ($ttp{$stn}{$j}{$item} = $ttp{$alt}{$j}{$item});
               }
            }
         }
         #  Fillable gap
         elsif ($n>0)
         {
            my $f0 = ($item eq 'PRCP') ? 0 : $ttp{$stn}{$jd0}{$item};
            my $f1 = ($item eq 'PRCP') ? 0 : ($ttp{$stn}{$jd}{$item}-$ttp{$stn}{$jd0}{$item})/($jd-$jd0);
            for (my $j=$jd0+1;$j<$jd;$j++)
            {
               $ttp{$stn}{$j}{$item} = $f1*($j-$jd0) + $f0;
            }
         }
         $jd0 = $jd;
      }
      #  Terminal gap
      for (my $j=$jd0+1;$j<$jd[-1];$j++)
      {
         foreach my $alt (@{$alt{$stn}})
         {
            exists($ttp{$alt}{$j}{$item}) && ($ttp{$stn}{$j}{$item} = $ttp{$alt}{$j}{$item});
         }
      }
   }
}

#
#  Print range of missing data
#
sub missing
{
   my ($stn,$item,$jd0,$jd1) = @_;
   my $n = $jd1-$jd0+1;
   if ($n==1)
   {
      my ($y0,$m0,$d0) = CalendarDay($jd0);
      warn "Missing $stn $item $m0/$d0/$y0\n";
   }
   elsif ($n>1)
   {
      my ($y0,$m0,$d0) = CalendarDay($jd0);
      my ($y1,$m1,$d1) = CalendarDay($jd1);
      warn "Missing $stn $item $n $m0/$d0/$y0 - $m1/$d1/$y1\n";
   }
}

#
#  See what is still missing
#
foreach my $stn (@stn)
{
   foreach my $item (sort keys %fill)
   {
      my ($jd0,@jd) = sort {$a <=> $b} keys %{$ttp{$stn}};
      $jd0 || next;
      foreach my $jd (@jd)
      {
         exists($ttp{$stn}{$jd}{$item}) || next;
         missing($stn,$item,exists($ttp{$stn}{$jd0}{$item})?$jd0+1:$jd0 , $jd-1);
         $jd0 = $jd;
      }
      missing($stn,$item,exists($ttp{$stn}{$jd0}{$item})?$jd0+1:$jd0 , $jd[-1]);
   }
}

#
#  Set missing values to long term average
#
foreach my $stn (@stn)
{
   my $dbf = DBF->open("../data/$name{$stn}.dbf") || die "Cannot open ../data/$name{$stn}.dbf\n";
   my ($jd0) = sort {$a <=> $b} keys %{$ttp{$stn}};
   (@ARGV>0 && !$jd0) && ($jd0 = JulianDay($ARGV[0],1,1));
   foreach my $jd ($jd0 .. $jd1)
   {
      foreach my $item (sort keys %fill)
      {
         if (!exists($ttp{$stn}{$jd}{$item}))
         {
            my ($y,$m,$d) = CalendarDay($jd);
            $ttp{$stn}{$jd}{$item} = $dbf->lookup('DATE',100*$m+$d,'N',$item);
         }
      }
   }
}

#
#  Trap TMIN>TMAX
#
foreach my $stn (@stn)
{
   foreach my $jd (sort {$a <=> $b} keys %{$ttp{$stn}})
   {
      ($ttp{$stn}{$jd}{TMIN} > $ttp{$stn}{$jd}{TMAX}) &&
         ($ttp{$stn}{$jd}{TMIN} = $ttp{$stn}{$jd}{TMAX} = ($ttp{$stn}{$jd}{TMIN} + $ttp{$stn}{$jd}{TMAX})/2);
   }
}

#
#  Save output to files by station
#
foreach my $stn (@stn)
{
   open(TXT , ">$name{$stn}.txt") || die "Cannot open file $name{$stn}.txt\n";
   print TXT "Date       Tmin Tmax Precip\n";
   foreach my $jd (sort {$a <=> $b} keys %{$ttp{$stn}})
   {
      my ($y,$m,$d) = CalendarDay($jd);
      printf TXT "%.2d/%.2d/%.4d %6.1f %6.1f %6.2f\n" , $m , $d , $y , @{$ttp{$stn}{$jd}}{'TMIN','TMAX','PRCP'};
   }
   close(TXT);
}
