#!/usr/bin/perl -w

#  Program to extract daily Min/Max Temperature and Precipitation
#  from NCDC TD3200 (Cooperative Summary of the Day) data sets
#  and fill missing values
#
#  Willem A. Schreuder
#  August 17, 2004

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

#  Where to end
my $jd1 = (@ARGV>0) ? JulianDay($ARGV[0],12,31) : 0;
 
#  TD3200 parameters to extract
my %par  = (TMIN=>1, TMAX=>2 , PRCP=>3);
my %unit = (PRCP=>'HI' , TMIN=>'F' , TMAX=>'F');
my %conv = ('HI'=>0.01 , 'F'=>1);

#  Alternate stations to fill data
my %alt  = (
   '050109'=>['059295','057950'],   #  Akron 4E > Yuma > Sterling
   '255310'=>['252065','251415'],   #  McCook > Culbertson > Cambridge
   '257070'=>['253037']);           #  RedCloud > Franklin #2

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

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

#
#  Read and process TD3200 input data file
#
my %ttp;
open(TD , "<TD3200") || die "Cannot open data file TD3200\n";
while (my $line = <TD>)
{
   #  Get station and data type
   my $stn  = substr($line,3,6);
   exists($load{$stn}) || next;
   my $type = substr($line,11,4);
   my $unit = substr($line,15,2);
   $unit =~ s/ //;
   exists($par{$type}) || next;
   ($unit{$type} eq $unit) || die "Odd units $type $unit\n$line\n";
   #  Get year, month and number of days
   my $yr = substr($line,17,4);
   my $mo = substr($line,21,2);
   my $n  = substr($line,27,3);
   #  Extract daily data
   for (my $k=0;$k<$n;$k++)
   {
      my $str = substr($line,30+12*$k,12);
      my $dy  = substr($str, 0,2);
      my $hr  = substr($str, 2,2);
      my $val = substr($str, 4,6);
      my $fl1 = substr($str,10,1);
      my $fl2 = substr($str,11,1);
      #  Validate data
      ($val !~ /[ -]\d\d\d\d\d/ || $val eq ' 99999' || $val eq '-99999' || $fl1 eq 'S' || $fl2 eq '2' || $fl2 eq '3') && next;
      #  Store by Julian date
      my $jd = JulianDay($yr,$mo,$dy);
      my $data = "$val$fl1$fl2";
      ($jd>$jd1) && ($jd1 = $jd);
      (exists($ttp{$stn}{$jd}{$type}) && (substr($ttp{$stn}{$jd}{$type},0,6) ne substr($data,0,6))) &&
         warn "Duplicate $stn $mo/$dy/$yr $jd $type: $ttp{$stn}{$jd}{$type} $data\n";
      $ttp{$stn}{$jd}{$type} = $data;
   }
}
close(TD);

#
#  Convert data to clean format
#
foreach my $stn (sort keys %ttp)
{
   foreach my $jd (sort {$a <=> $b} keys %{$ttp{$stn}})
   {
      exists($ttp{$stn}{$jd}{TMIN}) && ($ttp{$stn}{$jd}{TMIN} = sprintf '%d'   ,      substr($ttp{$stn}{$jd}{TMIN},0,6));
      exists($ttp{$stn}{$jd}{TMAX}) && ($ttp{$stn}{$jd}{TMAX} = sprintf '%d'   ,      substr($ttp{$stn}{$jd}{TMAX},0,6));
      exists($ttp{$stn}{$jd}{PRCP}) && ($ttp{$stn}{$jd}{PRCP} = sprintf '%.2f' , 0.01*substr($ttp{$stn}{$jd}{PRCP},0,6));
   }
}

#
#  Fill in gaps in all parameters
#
#  ET Stations  (twice)
foreach my $stn ((sort keys %alt),(sort keys %alt))
{
   #  Tmin, Tmax and PRCP
   foreach my $item (sort keys %fill)
   {
      my ($jd0,@jd) = sort {$a <=> $b} keys %{$ttp{$stn}};
      defined($jd0) || die "No $stn $item\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 (sort keys %alt)
{
   foreach my $item (sort keys %fill)
   {
      my ($jd0,@jd) = sort {$a <=> $b} keys %{$ttp{$stn}};
      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 0
#
foreach my $stn (sort keys %alt)
{
   foreach my $jd (sort {$a <=> $b} keys %{$ttp{$stn}})
   {
      foreach my $item (sort keys %fill)
      {
         exists($ttp{$stn}{$jd}{$item}) || ($ttp{$stn}{$jd}{$item} = 0);
      }
   }
}

#
#  Trap TMIN>TMAX
#
foreach my $stn (sort keys %alt)
{
   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 (keys %alt)
{
   open(TXT , ">$stn.txt") || die "Cannot open file $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 %4d %4d %6.2f\n" , $m , $d , $y , @{$ttp{$stn}{$jd}}{'TMIN','TMAX','PRCP'};
   }
   close(TXT);
}
