#!/usr/bin/perl -w

#  Program to extract monthly precipitation from GHCN
#  and write it to a DBF table

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

my $prism = PRISM->new('PPT','in');
my $xfm = Transform->SPCS('USGSRR');

#  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==2) || die "Usage:  $0 <first year> <last year>\n";
my @yr = ($ARGV[0] .. $ARGV[1]);
(@yr>0) || die "Invalid year range\n";
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');

my $dbg = DBF->new('ID:A7','YEAR:4','ANNUAL',map{uc $mo{$_}} @mo);
my $dbp = DBF->new('ID:A7','YEAR:4','ANNUAL',map{uc $mo{$_}} @mo);

#  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 (`zgrep PRCP /pm/raid7/t/NCDC/GHCN/by_year/$yr.csv.gz`)
   {
      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;
   }
}

#
#  Summarize monthly data
#
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;
         }
      }
   }
}

#
#  Dump monthly and annual values
#
foreach my $stn (@stn)
{
   my ($lon,$lat) = $xfm->inverse($stn{$stn}{X},$stn{$stn}{Y});
   foreach my $yr (@yr)
   {
      my @hdr = ("C$stn",$yr);
      my ($SUM,$sum,@REC,@rec);
      foreach my $mo (@mo)
      {
         # PRISM
         my $IN = sprintf '%.3f' , $prism->eval($lon,$lat,$yr,$mo);
         $SUM += $IN;
         push @REC , $IN;
         #  Fill missing months from PRISM
         my $in = $IN;
         if ($ppt{$stn}{$yr}{$mo}{Mis} <= $MaxMiss && $ppt{$stn}{$yr}{$mo}{Gap} <= $MaxCons)
         {
            # Convert 1/10 mm to inches
            $in = sprintf '%.3f' , $ppt{$stn}{$yr}{$mo}{P}/254;
         }
         $sum += $in;
         push @rec,$in;
      }
      $dbg->add(@hdr,$sum,@rec);
      $dbp->add(@hdr,$SUM,@REC);
   }
}
$dbg->format();
$dbg->write('GHCN.dbf');
$dbp->format();
$dbp->write('PRISM.dbf');
