#!/usr/bin/perl -w

#  Program to extract annual precipitation from
#  GHCN and PRISM data sets

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

my $prism = PRISM->new('PPT','in');
my $xfm = Transform->SPCS('USGSRR');
my $dbf = DBF->open('../data/dri.dbf');

#  Maximum number of missing days permitted
my $MaxMiss = 10;
my $MaxCons = 5;

my %stn = (
   '050109' =>  {ID=>'USC00050109',DRI=>'050109',X=> 480549,Y=>14607776,N=>'Akron 4 E'},
   '051121' =>  {ID=>'USC00051121',DRI=>'051121',X=> 710588,Y=>14263754,N=>'Burlington'},
   '051564' =>  {ID=>'USC00051564',DRI=>'051564',X=> 686112,Y=>14112695,N=>'Cheyenne Wells'},
   '054082' =>  {ID=>'USC00054082',DRI=>'054082',X=> 724056,Y=>14755644,N=>'Holyoke'},
   '054413' =>  {ID=>'USC00054413',DRI=>'054413',X=> 738747,Y=>14901009,N=>'Julesburg'},
   '059243' =>  {ID=>'USC00059243',DRI=>'059243',X=> 749903,Y=>14572326,N=>'Wray'},
   '141179' =>  {ID=>'USC00141179',DRI=>'141179',X=>1831189,Y=>14483074,N=>'Burr Oak 1 N'},
   '141699' =>  {ID=>'USC00141699',DRI=>'141699',X=>1055629,Y=>14314282,N=>'Colby 1SW'},
   '143527' =>  {ID=>'USC00143527',DRI=>'143527',X=>1545538,Y=>14113573,N=>'Hays 1 S'},
   '143837' =>  {ID=>'USC00143837',DRI=>'143837',X=>1230500,Y=>14292659,N=>'Hoxie'},
   '145363' =>  {ID=>'USC00145363',DRI=>'145363',X=>2009059,Y=>14213125,N=>'Minneapolis'},
   '145856' =>  {ID=>'USC00145856',DRI=>'145856',X=>1406128,Y=>14430033,N=>'Norton 9 SSE'},
   '145906' =>  {ID=>'USC00145906',DRI=>'145906',X=>1209838,Y=>14462978,N=>'Oberlin1 E'},
   '146378' =>  {ID=>'USC00146378',DRI=>'146374',X=>1546746,Y=>14441255,N=>'Phillipsburg 1 SSE'},
   '147093' =>  {ID=>'USC00147093',DRI=>'147093',X=> 848848,Y=>14453532,N=>'Saint Francis'},
   '148495' =>  {ID=>'USC00148495',DRI=>'148495',X=>1390477,Y=>14170432,N=>'Wakeeny'},
   '250640' =>  {ID=>'USC00250640',DRI=>'250640',X=>1407487,Y=>14575688,N=>'Beaver City'},
   '250810' =>  {ID=>'USC00250810',DRI=>'250810',X=>1464389,Y=>14714821,N=>'Bertrand'},
   '252065' =>  {ID=>'USC00252065',DRI=>'252065',X=>1128713,Y=>14616300,N=>'Culbertson'},
   '252690' =>  {ID=>'USC00252690',DRI=>'252690',X=>1394783,Y=>14703279,N=>'Elwood 8 S'},
   '253365' =>  {ID=>'USC00253365',DRI=>'253365',X=>1322774,Y=>14868020,N=>'Gothenburg'},
   '253735' =>  {ID=>'USC00253735',DRI=>'253735',X=>2036109,Y=>14595960,N=>'Hebron'},
   '253910' =>  {ID=>'USC00253910',DRI=>'253910',X=>1538380,Y=>14684054,N=>'Holdredge'},
   '254110' =>  {ID=>'USC00254110',DRI=>'254110',X=> 903844,Y=>14725259,N=>'Imperial'},
   '255090' =>  {ID=>'USC00255090',DRI=>'255090',X=> 935167,Y=>14845850,N=>'Madrid'},
   '255310' =>  {ID=>'USW00094040',DRI=>'255310',X=>1188038,Y=>14603001,N=>'McCook'},
   '255565' =>  {ID=>'USC00255565',DRI=>'255565',X=>1654313,Y=>14714193,N=>'Minden'},
   '256480' =>  {ID=>'USC00256480',DRI=>'256480',X=>1050642,Y=>14660550,N=>'Palisade'},
   '256585' =>  {ID=>'USC00256585',DRI=>'256585',X=> 993099,Y=>14941433,N=>'Paxton'},
   '257070' =>  {ID=>'USC00257070',DRI=>'257070',X=>1775580,Y=>14562825,N=>'Red Cloud'},
   '258255' =>  {ID=>'USC00258255',DRI=>'258255',X=>1016296,Y=>14588511,N=>'Stratton'},
   '258320' =>  {ID=>'USC00258320',DRI=>'258320',X=>1901742,Y=>14533481,N=>'Superior'},
   '258735' =>  {ID=>'USC00258735',DRI=>'258735',X=>1677566,Y=>14653524,N=>'Upland'},
   '259020' =>  {ID=>'USC00259020',DRI=>'259020',X=> 968206,Y=>14705184,N=>'Wauneta 3 NW'},
   );
my @stn = sort keys %stn;
my %ren;
foreach my $stn (@stn)
{
   $ren{$stn{$stn}{ID}} = $stn;
}

#  Process command line
(@ARGV>0) || die "Usage:  $0 <year>\n";
my $yrs = shift @ARGV;
my @yr = ($yrs =~/^(\d+)-(\d+)$/) ? ($1 .. $2) : ($yrs);
my %mo = ('01'=>'Jan','02'=>'Feb','03'=>'Mar','04'=>'Apr','05'=>'May','06'=>'Jun','07'=>'Jul','08'=>'Aug','09'=>'Sep','10'=>'Oct','11'=>'Nov','12'=>'Dec');
my @mo = ('01','02','03','04','05','06','07','08','09','10','11','12');


#  Initialize all stations, years and months
my %ppt;
foreach my $stn (@stn)
{
   foreach my $yr (@yr)
   {
      foreach my $mo (@mo)
      {
         
         $ppt{$stn}{$yr}{$mo} = {N=>0,P=>0};
      }
   }
}

#
#  Read GHCN data
#
foreach my $yr (@yr)
{
   #  Get only precipitation entries from GHCN data
   foreach my $line (`grep PRCP GHCN`)
   {
      my ($id,$date,$elem,$val) = split(/,/,$line);
      #  Retain only selected stations
      exists($ren{$id}) || next;
      #  Map GHCN id to station number
      my $stn = $ren{$id};
      #  Pause for paranoia
      ($elem eq 'PRCP') || die "Invalid element $elem\n";
      #  Decode date
      my $yr = substr($date,0,4);
      my $mo = substr($date,4,2);
      my $dy = substr($date,6,2);
      #  Sum precip in 1/10 mm
      $ppt{$stn}{$yr}{$mo}{N}++;
      $ppt{$stn}{$yr}{$mo}{D}{$dy}++;
      $ppt{$stn}{$yr}{$mo}{P} += $val;
   }
}

#
#  Check for missing days
#
foreach my $stn (@stn)
{
   foreach my $yr (@yr)
   {
      foreach my $mo (@mo)
      {
         $ppt{$stn}{$yr}{$mo}{Mis} = DaysInMonth($yr,$mo) - $ppt{$stn}{$yr}{$mo}{N};
         if ($ppt{$stn}{$yr}{$mo}{N}>0)
         {
            my @dy = sort keys %{$ppt{$stn}{$yr}{$mo}{D}};
            my $gap =  $dy[0]-1;
            while (@dy>1)
            {
               ($dy[1]-$dy[0]-1>$gap) && ($gap = $dy[1]-$dy[0]-1);
               shift @dy;
            }
            (DaysInMonth($yr,$mo)-$dy[0]-1 > $gap) && ($gap = DaysInMonth($yr,$mo)-$dy[0]-1);
            $ppt{$stn}{$yr}{$mo}{Gap} = $gap;
         }
         else
         {
            $ppt{$stn}{$yr}{$mo}{Gap} = 0;
         }
      }
   }
}

#
#  Fill incomplete or missing months from PRISM
#
my %fill;
foreach my $stn (@stn)
{
   my ($lon,$lat) = $xfm->inverse($stn{$stn}{X},$stn{$stn}{Y});
   foreach my $yr (@yr)
   {
      my $sum;
      foreach my $mo (@mo)
      {
         #  Fill missing months from PRISM
         if ($ppt{$stn}{$yr}{$mo}{Mis} > $MaxMiss || $ppt{$stn}{$yr}{$mo}{Gap} > $MaxCons)
         {
            my $in = $prism->eval($lon,$lat,$yr,$mo);
            my $mean = $dbf->lookup('ID',$stn{$stn}{DRI},'N',$mo{$mo}) || die "Missing data $stn $mo{$mo}\n";
            $fill{$stn}{$yr}{$mo} = {N=>$ppt{$stn}{$yr}{$mo}{N},M=>$ppt{$stn}{$yr}{$mo}{Mis},C=>$ppt{$stn}{$yr}{$mo}{Gap},P=>$ppt{$stn}{$yr}{$mo}{P}/254,R=>$in,D=>$mean};
            $sum += defined($in) ? $in : $mean;
         }
         else
         # Convert 1/10 mm to inches
         {
            $sum += $ppt{$stn}{$yr}{$mo}{P}/254;
         }
      }
      $ppt{$stn}{$yr} = $sum;
   }
}

open(HTM , ">fill.htm") || die "Cannot open file fill.htm\n";
print HTM "<table border>\n";
print HTM "<tr><th colspan=2>Station</th><th rowspan=2>Month</th><th colspan=2>Missing Days</th><th colspan=3>Precipitation (in)</th>\n";
print HTM "<tr><th>ID</th><th>Name</th><th>Total</th><th>Consecutive</th><th>GHCN</th><th>PRISM</th><th>MEAN</th>\n";
my $k=0;
foreach my $stn (sort keys %fill)
{
   $k++ && printf HTM "<tr><td colspan=7>&nbsp;</td></tr>\n";
   foreach my $yr (sort keys %{$fill{$stn}})
   {
      my @mo = sort keys %{$fill{$stn}{$yr}};
      my $n = scalar(@mo);
      foreach my $mo (@mo)
      {
         printf HTM "<tr align=right>";
         ($mo == $mo[0]) && print HTM "<td rowspan=$n align=left>$stn</td><td rowspan=$n align=left>$stn{$stn}{N}</td>";
         print HTM "<td align=center>$mo{$mo}</td><td>$fill{$stn}{$yr}{$mo}{M}</td>";
         if ($fill{$stn}{$yr}{$mo}{N})
         {
            printf HTM "<td>$fill{$stn}{$yr}{$mo}{C}</td><td>%.3f</td>" , $fill{$stn}{$yr}{$mo}{P};
         }
         else
         {
            print HTM "<td>---</td><td>---</td>";
         }
         if (defined($fill{$stn}{$yr}{$mo}{R}))
         {
            printf HTM "<td bgcolor=yellow>%.3f</td><td>%.3f</td></tr>\n" , $fill{$stn}{$yr}{$mo}{R} , $fill{$stn}{$yr}{$mo}{D};
         }
         else
         {
            printf HTM "<td>---</td><td bgcolor=yellow>%.3f</td></tr>\n" , $fill{$stn}{$yr}{$mo}{D};
         }
      }
   }
}
print HTM "</table>\n";
close(HTM);

#  Output data
foreach my $yr (@yr)
{
   print $yr;
   foreach my $stn (@stn)
   {
      printf "%8.2f" , $ppt{$stn}{$yr};
   }
   print "\n";
}
