#!/usr/local/bin/perl
#
# program to generate snotel graphs
#
#

use GD; 
$numargs = @ARGV;
# no parameters.. show instructions.
if ($numargs == 0) {
   print ("  \n");
   print ("To invoke interactively :\n");
   print ("program usage : g1snot.pl ra\n");
   print ("  \n");
   print ("Command line usage : not available yet\n");
   print ("program usage : g1snot.pl station sdate edate XXX\n");
   print ("  station - snotel shef id\n");
   print ("  sdate - start date (yyyymmdd)\n");
   print ("  edate - end date (yyyymmdd)\n");
   print ("  XXX - PCPN, SWE, TEMP. 1=yes. 0=no\n");
   print ("  \n");
   die ("Try again.\n");
}

# invoked interatively?
if ($numargs == 1) {
   $param = $ARGV[0];
   if ($param eq "ra") {
      $param = 1;
# request user input
# set defaults if no input
      print ("Enter the Station Shef id (Def = SEPI1) :");
      $station = <STDIN>;
      $station = "\U$station\E";
      chop ($station);
      if (length($station) < 1) { $station = "SEPI1"; }
      print ("Enter the Start Date (YYYYMMDD Def = 19961001) :");
      $sdate = <STDIN>;
      chop ($sdate);
      if (length($sdate) < 1) { $sdate = "19961001"; }
      print ("Enter the End Date (YYYYMMDD Def = 19970201) :");
      $edate = <STDIN>;
      chop ($edate);
      if (length($edate) < 1) { $edate = "19970201"; }
   } else {
      print ("$numargs \n");
      print ("Improperly started.  Instructions follow.\n");
      print ("  \n");
      print ("To invoke interactively :\n");
      print ("program usage : g1snot.pl ra\n");
      print ("  \n");
      print ("Command line usage : not available yet\n");
#      print ("program usage : g1snot.pl \n");
      print ("  \n");
      die ("Try again.\n");
   }
}

# invoked with command line args
if ($numargs == 11) {
   $station = $ARGV[0];
   $sdate = $ARGV[1];
   $edate = $ARGV[2];
   $dtype = $ARGV[3];
   $ymn = $ARGV[4];
   $ymx = $ARGV[5];
   $y2offset = $ARGV[6];
   $ymx2 = $ymx - $y2offset;
   $ymn2 = $ymn - $y2offset;
   $ymn3 = $ARGV[7];
   $ymx3 = $ARGV[8];
   $Plot = $ARGV[9];
   $Pplot = substr($Plot, 0 ,1);
   $Splot = substr($Plot, 1 ,1);
   $Tplot = substr($Plot, 2 ,1);
   $Sauto = $ARGV[10];
#print ("$station $sdate $edate $dtype $ymn $ymx $y2offset $ymn3 $ymx3 $Plot $Sauto\n");
}

# Number of days between start and end
$days = &ndays($sdate, $edate);
#print ("$days\n");

# hourly or daily if < 30 days
if (($days < 30) && ($param == 1)) {
   print ("(H)ourly or (D)aily values (Def = D) :");
   $dtype = <STDIN>;
   chop ($dtype);
   $dtype =~ tr/a-z/A-Z/;
   if (($dtype ne "H") && ($dtype ne "h")) { $dtype = "D"; }
} else {
   $dtype =~ tr/a-z/A-Z/;
   if ($days > 25) { $dtype = "D"; }
   if (($dtype ne "H") && ($dtype ne "h")) { $dtype = "D"; }
}   

# Initialize the graph data arrays to missing values
if ($dtype eq "H") {
   $npoints = $days * 24;
} else {
   $npoints = $days;
}
for ($i=0; $i<=$npoints; $i++) {
   $pcpn[$i] = "-999";
   $swec[$i] = "-999";
   $obat[$i] = "-999";
   $mxat[$i] = "-999";
   $mnat[$i] = "-999";
   $avat[$i] = "-999";
}

# extract Station Metadata info from inventory file
$filename = "/wrcc4/krwrcc/inventories/snotinv.shl";
unless (open(SNOTINV, $filename )) {
   if (!(-e $filename)) {
      die ("File $filename does not exist.\n");
   } elsif (!(-r $filename)) {
      die ("You are not allowed to read $filename.\n");
   } else {
      die ("File $filename can not be opened.\n");
   }
}
while (! eof(SNOTINV)) {
   $line = <SNOTINV>;
   if (substr($line, 0, 5) eq $station) {
      $snotelmeta = $line;
   }
}
close (SNOTINV);

# Set state info
if (substr($snotelmeta, 13, 2) eq "02") { $State = "AZ";}
if (substr($snotelmeta, 13, 2) eq "04") { $State = "CA";}
if (substr($snotelmeta, 13, 2) eq "05") { $State = "CO";}
if (substr($snotelmeta, 13, 2) eq "10") { $State = "ID";}
if (substr($snotelmeta, 13, 2) eq "24") { $State = "MT";}
if (substr($snotelmeta, 13, 2) eq "26") { $State = "NV";}
if (substr($snotelmeta, 13, 2) eq "29") { $State = "NM";}
if (substr($snotelmeta, 13, 2) eq "35") { $State = "OR";}
if (substr($snotelmeta, 13, 2) eq "39") { $State = "SD";}
if (substr($snotelmeta, 13, 2) eq "42") { $State = "UT";}
if (substr($snotelmeta, 13, 2) eq "45") { $State = "WA";}
if (substr($snotelmeta, 13, 2) eq "48") { $State = "WY";}

# Open snotel station data file
$stwrcc = substr($snotelmeta, 13, 2) . "0" . substr($snotelmeta, 15, 3);
#$filename = "/wrcc10b/krwrcc/snotel/srfls/tx/$stwrcc";
$filename = "/clihome2/gmwrcc/programs/snotexbin.e $stwrcc $sdate";
open(SNOTEL, "$filename|") || die ("Can't retrieve data\n");
#unless (open(SNOTEL, $filename )) {
#   if (!(-e $filename)) {
#      die ("File $filename does not exist.\n");
#   } elsif (!(-r $filename)) {
#      die ("You are not allowed to read $filename.\n");
#   } else {
#      die ("File $filename can not be opened.\n");
#   }
#}

# initialize extremes
$pcpnmx = -999;
$pcpnmn = 999;
$swecmx = -999;
$swecmn = 999;
$mxatmx = -999;
$mnatmn = 999;

# load data
@snotdata = <SNOTEL>;
close (SNOTEL);
# scan data arrays and load data from desired interval
# into working graph arrays
# BE SURE to cross load hour 00 and hour 24.
foreach $snotdata(@snotdata) {
   $curdate = substr($snotdata, 6, 8);
   # is this hour 00 or hour 24?
   $curhour = substr($snotdata, 14, 2);
   if (($curhour eq "24") || ($curhour eq "00") || ($dtype eq "H")) {
      # is date of current line between the start and end date?
      if (($curdate >= $sdate) && ($curdate <= $edate)) {
         # Yes, now find the array index that would 
         # correspond to this date and time.
         if ($dtype eq "H") {
            # displaying hourly data
            $index = &ndays($sdate, $curdate);
            $index = (($index-1) * 24) + $curhour;
            if ($index > (($days + 1) * 24)) { 
               die ("date conversion error in program (Hourly index)\n"); 
            }
         } else {
            # displaying daily data
            $index = &ndays($sdate, $curdate);
            if ($index > $days) {
               die ("date conversion error in program (Daily index)\n");
            }
         }
         # Load the graph arrays (don't load missing values from current line)
         # PCPN
         if (substr($snotdata, 23, 1) ne "M") { 
            $pcpn[$index] = substr($snotdata, 18, 5) /100;
            if ($pcpn[$index] < $pcpnmn) { $pcpnmn = $pcpn[$index]; }
            if ($pcpn[$index] > $pcpnmx) { $pcpnmx = $pcpn[$index]; }
         }
         # SWE 
         if (substr($snotdata, 37, 1) ne "M") { 
            $swec[$index] = substr($snotdata, 32, 5) /10;
            if ($swec[$index] < $swecmn) { $swecmn = $swec[$index]; }
            if ($swec[$index] > $swecmx) { $swecmx = $swec[$index]; }
         }
         # MXAT
         if (substr($snotdata, 44, 1) ne "M") { 
            $mxat[$index] = substr($snotdata, 39, 5);
            if ($mxat[$index] > $mxatmx) { $mxatmx = $mxat[$index]; }
         }
         # MNAT
         if (substr($snotdata, 51, 1) ne "M") { 
            $mnat[$index] = substr($snotdata, 46, 5);
            if ($mnat[$index] < $mnatmn) { $mnatmn = $mnat[$index]; }
         }
         # AVAT
         if (substr($snotdata, 58, 1) ne "M") { 
            $avat[$index] = substr($snotdata, 53, 5);
         }
         # OBAT
         if (substr($snotdata, 65, 1) ne "M") { 
            $obat[$index] = substr($snotdata, 60, 5);
         }
#    print ("$pcpn[$index]\n");
#    print ("$index\n");
      }
   }
}   

# determine y-axis scale limits if auto scale option is selected.
if ($Sauto eq "1") {
   # if Pcpn will be plotted, use as controlling limits.
   if ($Pplot eq "1") {
      $mx = $pcpnmx;
      $mn = $pcpnmn;
      if ($Splot eq "1") {
         if ($swecmx > $mx) { $mx = $swecmx; }
         if ($swecmn < $mn) { $mn = $swecmn; }
      }
      if (($mx - $mn) < 10) {
         $ystep = 1;
      } elsif (($mx - $mn) < 20) {
         $ystep = 2;
      } elsif (($mx - $mn) < 50) {
         $ystep = 5;
      } elsif (($mx - $mn) < 200) {
         $ystep = 10;
      } else {
         $ystep = 20;
      }
      for ($i=0; $i<$ystep; $i++) {
         if ((int($mn-$i)/$ystep) == int((int($mn-$i)/$ystep))) {
            $ymn = int($mn) - $i;
         }
      }
      for ($i=1; $i<$ystep+1; $i++) {
         if ((int($mx+$i)/$ystep) == int((int($mx+$i)/$ystep))) {
            $ymx = int($mx) + $i;
         }
      }
     print ("$mx $mn $ystep $ymn $ymx\n");
      if ($ymn < 1) { $ymn = 0; }
      $y2offset = 0;
      $ymx2 = $ymx - $y2offset;
      $ymn2 = $ymn - $y2offset;
      $ymx3 = int(($mxatmx + 10)/10) * 10;
      $i = int(($ymx3 - $mnatmn)/($ymx - $ymn)) + 1;
      $ymn3 = $ymx3 - ($i * ($ymx - $ymn));
   }
   # if Pcpn not plotted, but SWE is use as controlling limits.
   if (($Splot eq "1") && ($Pplot eq "0")) {
      if (($swecmx - $swecmn) < 10) {
         $ystep = 1;
      } elsif (($swecmx - $swecmn) < 20) {
         $ystep = 2;
      } elsif (($swecmx - $swecmn) < 50) {
         $ystep = 5;
      } elsif (($swecmx - $swecmn) < 200) {
         $ystep = 10;
      } else {
         $ystep = 20;
      }
      for ($i=0; $i<$ystep; $i++) {
         if ((int($swecmn-$i)/$ystep) == int((int($swecmn-$i)/$ystep))) {
            $ymn = int($swecmn) - $i;
         }
      }
      for ($i=0; $i<$ystep; $i++) {
         if ((int($swecmx+$i)/$ystep) == int((int($swecmx+$i)/$ystep))) {
            $ymx = int($swecmx) + $i;
         }
      }
      if ($ymn < 1) { $ymn = 0; }
      $y2offset = 0;
      $ymx2 = $ymx - $y2offset;
      $ymn2 = $ymn - $y2offset;
      $ymx3 = int(($mxatmx + 10)/10) * 10;
      $i = int(($ymx3 - $mnatmn)/($ymx - $ymn)) + 1;
      $ymn3 = $ymx3 - ($i * ($ymx - $ymn));
   }
   # if Pcpn and SWE not plotted, but Temp is use as controlling limits.
   if (($Tplot eq "1") && ($Splot eq "0") && ($Pplot eq "0")) {
      if (($mxatmx - $mnatmn) < 10) {
         $ystep = 1;
      } elsif (($mxatmx - $mnatmn) < 20) {
         $ystep = 2;
      } elsif (($mxatmx - $mnatmn) < 50) {
         $ystep = 5;
      } elsif (($mxatmx - $mnatmn) < 200) {
         $ystep = 10;
      } else {
         $ystep = 20;
      }
      for ($i=0; $i<$ystep; $i++) {
         if ((int($mnatmn-$i)/$ystep) == int((int($mnatmn-$i)/$ystep))) {
            $ymn3 = int($mnatmn) - $i;
         }
      }
      for ($i=0; $i<$ystep; $i++) {
         if ((int($mxatmx+$i)/$ystep) == int((int($mxatmx+$i)/$ystep))) {
            $ymx3 = int($mxatmx) + $i;
         }
      }
   }
}   
print ("$station $sdate $edate $dtype $ymn $ymx $y2offset $ymn3 $ymx3 $Plot $Sauto\n");
print ("$pcpnmn $pcpnmx $swecmn $swecmx $mnatmn $mxatmx\n");

# Set Graph Boundaries
if ($param == 1) {
      print ("Extremes for selected period.\n");
      print ("Pcpn Min = $pcpnmn\n");
      print ("Pcpn Max = $pcpnmx\n");
      print ("SWE  Min = $swecmn\n");
      print ("SWE  Max = $swecmx\n");
      print ("AirT Min = $mnatmn\n");
      print ("AirT Max = $mxatmx\n");
      # Set Graph Boundaries
      # PCPN scale
      print ("Enter the Precipitation Scale Min Val. (Def = 0) :");
      $ymn = <STDIN>;
      chop ($ymn);
      if (length($ymn) < 1) { $ymn = "0"; }
      $ymn = int($ymn);
      print ("Enter the Precipitation Scale Max Val. (Def = 50) :");
      $ymx = <STDIN>;
      chop ($ymx);
      if (length($ymx) < 1) { $ymx = "50"; }
      $ymx = int($ymx);
      # SWE scale
      print ("Enter an offset for the SWE scale Val. (Def = 0)  \n");
      print ("ie. if the PCPN min is 20 and a 5 SWE min        \n");
      print ("    is desired then enter   15                  :");
      $y2offset = <STDIN>;
      chop ($y2offset);
      if (length($y2offset) < 1) { $y2offset = "0"; }
      $ymx2 = $ymx - $y2offset;
      $ymn2 = $ymn - $y2offset;
      # TEMP scale
      print ("Enter the Temperature Scale Min Val. (Def = -20) :");
      $ymn3 = <STDIN>;
      if (substr($ymn3, 0, 1) eq "\n")  { 
         $ymn3 = "-20";
      } else {
         chop ($ymn3);
         $ymn3 = int($ymn3);
      }
      print ("Enter the Temperature Scale Max Val.  (Def = 80) :");
      $ymx3 = <STDIN>;
      if (substr($ymx3, 0, 1) eq "\n")  { 
         $ymx3 = "80";
      } else {
         chop ($ymx3);
         $ymx3 = int($ymx3);
      }
   print ("$station $sdate $edate $ymx $ymn $ymx2 $ymn2 $ymx3 $ymn3\n");
}

if ($dtype eq "H") {
   $xmx = $days * 24;
} else {
   $xmx = $days;
}
$xmn = 0;

# compute y-axis ratio and x-axis ratio
$yratio = 190 / ($ymx - $ymn);
$yratio2 = 190 / ($ymx2 - $ymn2);
$yratio3 = 190 / ($ymx3 - $ymn3);
$xratio = 367 / ($xmx - $xmn);
#   print ("$yratio $ymx $ymn\n");

# Open GIF Image
$im = new GD::Image(510,290);

# allocate colors
# allocate white (also becomes background color)
$white = $im->colorAllocate(255, 255, 255);
# allocate black
$black = $im->colorAllocate(0, 0, 0);
# allocate red
$red = $im->colorAllocate(255, 0, 0);
# allocate green
$green = $im->colorAllocate(0,255,0);
# allocate yellow
$yellow = $im->colorAllocate(255,255,0);
# allocate blue
$blue = $im->colorAllocate(0,0,255);
# allocate gray
$gray = $im->colorAllocate(191,191,191);
# allocate dark gray
$dgray = $im->colorAllocate(207,207,207);

# draw frame and graph area and shade graph area
$im->interlaced('true');
$im->rectangle(0,0,509,289,$black);
$im->rectangle(70,38,437,228,$black);
#$im->fill(100,200,$gray);

# compute Y-axis minor intervals
if (($Pplot eq "1") || ($Splot eq "1")) {
   $intleft = $ymx - $ymn;
   $intminor = 1;
   if ($intleft > 10) { $intminor = 2; }
   if ($intleft > 20) { $intminor = 5; }
   if ($intleft > 50) { $intminor = 10; }
} elsif ($Tplot eq "1") {
   $intleft = $ymx3 - $ymn3;
   $intminor = 1;
   if ($intleft > 10) { $intminor = 2; }
   if ($intleft > 20) { $intminor = 5; }
   if ($intleft > 50) { $intminor = 10; }
}

if (($Pplot eq "1") || ($Splot eq "1")) {
   # Y-axis left legend
   $ileft = 0;
   $intleft = $ymx - $ymn;
   $label = $ymn;
   while ($ileft <= $intleft) {
      while (length($label) < 3) { $label = " " .$label; }
      $im->string(gdTinyFont,50,223 - int($ileft*$yratio),"$label",$black);
      $ileft = $ileft + $intminor;
      $label = $label + $intminor;
   }
}

# Y-axis left title
if (($Pplot eq "1") && ($Splot eq "1")) {
   $im->stringUp(gdMediumBoldFont,20,210,"Precip. ",$green);
   $im->stringUp(gdMediumBoldFont,20,154,"/",$black);
   $im->stringUp(gdMediumBoldFont,20,147," SWC ",$blue);
   if ($y2offset == 0) {
      $im->stringUp(gdMediumBoldFont,20,112,"(in.)",$black);
   } else {
      if (length($y2offset) > 3) { $y2offset = substr($y2offset, 0, 3); }
      $label = "$y2offset";
      if ($y2offset > 0) { $label = "+" . $label; }
      $im->stringUp(gdMediumBoldFont,20,112,$label,$blue); 
      $im->stringUp(gdMediumBoldFont,20,70,"(in.)",$black);
   }
}
if (($Pplot eq "1") && ($Splot eq "0")) {
   $im->stringUp(gdMediumBoldFont,20,210,"Precipitation  ",$green);
   $im->stringUp(gdMediumBoldFont,20,98,"(in.)",$black);
}
if (($Pplot eq "0") && ($Splot eq "1")) {
   $im->stringUp(gdMediumBoldFont,20,210,"Snow Water Content ",$blue);
   $im->stringUp(gdMediumBoldFont,20,72,"(in.)",$black);
}

if ($Tplot eq "1") {
   # Y-axis right legend
   $ileft = 0;
   $intleft3 = $ymx3 - $ymn3;
   $label = $ymn3;
   $intminor3 = $intleft3 / ($intleft/$intminor) ;
   while ($ileft <= ($intleft3 + .01)) {
      while (length($label) < 3) { $label = " " .$label; }
      if (length($label) > 5) { $label = substr($label, 0, 5); }
      $im->string(gdTinyFont,440,223 - int($ileft*$yratio3),"$label",$black);
      $ileft = $ileft + $intminor3;
      $label = $label + $intminor3;
   }

   # Y-axis right title
   $im->stringUp(gdMediumBoldFont,475,190,"Temperature (F)",$red);
}

# Y-axis Grid
$ileft = 0;
#$intleft = $ymx - $ymn;
while ($ileft < $intleft) {
   $im->line(437,228 - int($ileft*$yratio),70,228 - int($ileft*$yratio),$dgray);
   $ileft = $ileft + $intminor;
   $label = $label + $intminor;
}

# X-axis title                                             
$im->string(gdMediumBoldFont,210,250,"Day of Year",$black);

# Compute x-grid (date) minor
$xminor = 1;
if ($days > 7) { $xminor = 2; }
if ($days > 14) { $xminor = 5; }
if ($days > 31) { $xminor = 10; }
if ($days > 65) { $xminor = 15; }
if ($days > 185) { $xminor = 32; }
if ($days > 365) { $xminor = 92; }
if ($days > 730) { $xminor = 184; }
if ($days > 900) { $xminor = 365; }

# X-axis legend
$cdate = $sdate;
while ($cdate < $edate) {
   $cdateyr = substr($cdate, 0, 4);
   $cdatemo = substr($cdate, 4, 2);
   $cdatedy = substr($cdate, 6, 2);
   $label = "$cdatemo/$cdatedy";
   $num = length($label);
   $offset = int($num*2.5);
   $label2 = $cdateyr;
   $num2 = length($label2);
   $offset2 = int($num2*2.5);
   $cdays = &ndays($sdate, $cdate);
#  print ("$cdays $cdate\n");
   if ($dtype eq "H") {
      $count = $cdays * 24;
   } else {
      $count = $cdays;
   }
   $im->string(gdTinyFont, 70 + int((($count - $xmn) * $xratio) + 0.0) - $offset, 232, $label, $black);
   $im->string(gdTinyFont, 70 + int((($count - $xmn) * $xratio) + 0.0) - $offset2, 240, $label2, $black);

# X-axis Grid
   $im->dashedLine(70 + int((($count - $xmn) * $xratio) + 0.0),38,70 + int((($count - $xmn) * $xratio) + 0.0),228,$dgray);
   # increment date for next label   
   $cdate = &incdate($cdatedy, $cdatemo, $cdateyr, $xminor);
}

# print end date label
   $cdateyr = substr($edate, 0, 4);
   $cdatemo = substr($edate, 4, 2);
   $cdatedy = substr($edate, 6, 2);
   $label = "$cdatemo/$cdatedy";
   $num = length($label);
   $offset = int($num*2.5);
   $label2 = $cdateyr;
   $num2 = length($label2);
   $offset2 = int($num2*2.5);
   $cdays = &ndays($sdate, $edate);
   if ($dtype eq "H") {
      $count = $cdays * 24;
   } else {
      $count = $cdays;
   }
   $im->string(gdTinyFont, 70 + int((($count - $xmn) * $xratio) + 0.0) - $offset, 232, $label, $black);
   $im->string(gdTinyFont, 70 + int((($count - $xmn) * $xratio) + 0.0) - $offset2, 240, $label2, $black);

# Legends Box
$im->arc(80,273,10,10,180,270,$black);
$im->line(80,268,425,268,$black);
$im->arc(425,273,10,10,270,360,$black);
$im->line(430,273,430,280,$black);
$im->arc(425,280,10,10,0,90,$black);
$im->line(425,285,80,285,$black);
$im->arc(80,280,10,10,90,180,$black);
$im->line(75,280,75,273,$black);
#$im->fill(100,270,$gray);

# Write Legends
if ($Pplot eq "1") {
   $im->string(gdSmallFont,94,271,"Acc. Precipitation",$black);
   $im->line(87,274,83,278,$green);
   $im->line(83,274,87,278,$green);
}
if ($Splot eq "1") {
   $im->string(gdSmallFont,225,271,"Snow Water Content",$black);
   $im->line(216,276,220,276,$blue);
   $im->line(218,274,218,278,$blue);
}
if ($Tplot eq "1") {
   $im->string(gdSmallFont,355,271,"Temp. Mx/Mn",$black);
   $im->line(340,276,350,276,$red);
}

# Write WRCC id
$im->string(gdTinyFont,460,250,"Western",$black);
$im->string(gdTinyFont,458,260,"Regional",$black);
$im->string(gdTinyFont,460,270,"Climate",$black);
$im->string(gdTinyFont,463,280,"Center",$black);

# Write USDA/NRCS Acknowledgement
$im->string(gdTinyFont,4,250,"Provisional",$black);
$im->string(gdTinyFont,15,260,"data",$black);
$im->string(gdTinyFont,4,270,"provided by",$black);
$im->string(gdTinyFont,10,280,"USDA/NRCS",$black);

# Titles
$title1 = substr($snotelmeta, 19, 21);
$title1 =~ s/[\t ]+$//;
$title1 = "\L$title1\E";
$title1 =~ s/\b[a-z]/\u$&/g;
$title1 = $title1 . ", " . $State . "   Snotel Site";
$num = length($title1);
$im->string(gdLargeFont,253-($num * 4),8,$title1,$black);
$title2 = "Lat ";
$title2 = $title2 . substr($snotelmeta, 40, 2) . " " . substr($snotelmeta, 42, 2) ."' N";
$title2 = $title2 . "   Lon ";
$title2 = $title2 . substr($snotelmeta, 45, 3) . " " . substr($snotelmeta, 48, 2) ."' W";
$title2 = $title2 . "   Elev ";
$title2 = $title2 . substr($snotelmeta, 51, 5) . " feet";
$num = length($title2);
#$im->string(gdSmallFont,253-int($num*3),28,$title2,$black);
$im->string(gdTinyFont,55,25,$title2,$black);
$title3 = "NRCS ID " . substr($snotelmeta, 6, 6) . "   NWS ID " . substr($snotelmeta, 0, 5);
$num = length($title3);
$im->string(gdTinyFont,310,25,$title3,$black); 

# Display data
# Display points
if ($dtype eq "H") {
   $npoints = $days * 24;
} else {
   $npoints = $days;
}
for ($i=0; $i<=$npoints; $i++) {
   # PCPN  (x symbol)
   if (($pcpn[$i] ne "-999") && ($Pplot eq "1")) {
         $im->setPixel(71 + int((($i - $xmn) * $xratio) + 0.5),228 - int((($pcpn[$i] - $ymn)*$yratio)+0.5) ,$green);
         $im->setPixel(72 + int((($i - $xmn) * $xratio) + 0.5),229 - int((($pcpn[$i] - $ymn)*$yratio)+0.5) ,$green);
         $im->setPixel(70 + int((($i - $xmn) * $xratio) + 0.5),227 - int((($pcpn[$i] - $ymn)*$yratio)+0.5) ,$green);
         $im->setPixel(70 + int((($i - $xmn) * $xratio) + 0.5),229 - int((($pcpn[$i] - $ymn)*$yratio)+0.5) ,$green);
         $im->setPixel(72 + int((($i - $xmn) * $xratio) + 0.5),227 - int((($pcpn[$i] - $ymn)*$yratio)+0.5) ,$green);
   }
   # SWE   (+ symbol)
   if (($swec[$i] ne "-999") && ($Splot eq "1")) {
         $im->setPixel(71 + int((($i - $xmn) * $xratio) + 0.5),228 - int((($swec[$i] - $ymn2)*$yratio2)+0.5) ,$blue);
         $im->setPixel(72 + int((($i - $xmn) * $xratio) + 0.5),228 - int((($swec[$i] - $ymn2)*$yratio2)+0.5) ,$blue);
         $im->setPixel(70 + int((($i - $xmn) * $xratio) + 0.5),228 - int((($swec[$i] - $ymn2)*$yratio2)+0.5) ,$blue);
         $im->setPixel(71 + int((($i - $xmn) * $xratio) + 0.5),229 - int((($swec[$i] - $ymn2)*$yratio2)+0.5) ,$blue);
         $im->setPixel(71 + int((($i - $xmn) * $xratio) + 0.5),227 - int((($swec[$i] - $ymn2)*$yratio2)+0.5) ,$blue);
   }
   # Daily max and min temp (vertical bar)
   if (($mxat[$i] ne "-999") && ($mnat[$i] ne "-999") && ($Tplot eq "1")) {
         $im->line(71 + int((($i - $xmn) * $xratio) + 0.5),228 - int((($mxat[$i] - $ymn3)*$yratio3)+0.5) ,71 + int(((($i) - $xmn) * $xratio) + 0.5),228 - int((($mnat[$i] - $ymn3)*$yratio3)+0.5) ,$red);
   }
   # Hourly temperature
   if (($obat[$i] ne "-999") && ($dtype eq "H") && ($Tplot eq "1")) {
         $im->setPixel(71 + int((($i - $xmn) * $xratio) + 0.5),228 - int((($obat[$i] - $ymn3)*$yratio3)+0.5) ,$red);
   }
}
# Display connecting lines for pcpn and swe if plotting daily data
if ($dtype ne "H") {
   for ($i=0; $i<$days; $i++) {
      # PCPN
      if (($pcpn[$i] ne "-999") && ($pcpn[$i+1] ne "-999") && ($Pplot eq "1")) {
            $im->line(71 + int((($i - $xmn) * $xratio) + 0.5),228 - int((($pcpn[$i] - $ymn)*$yratio)+0.5) ,71 + int(((($i+1) - $xmn) * $xratio) + 0.5),228 - int((($pcpn[$i+1] - $ymn)*$yratio)+0.5) ,$green);
      }
      # SWE
      if (($swec[$i] ne "-999") && ($swec[$i+1] ne "-999") && ($Splot eq "1")) {
            $im->line(71 + int((($i - $xmn) * $xratio) + 0.5),228 - int((($swec[$i] - $ymn2)*$yratio2)+0.5) ,71 + int(((($i+1) - $xmn) * $xratio) + 0.5),228 - int((($swec[$i+1] - $ymn2)*$yratio2)+0.5) ,$blue);
      }
   }
}

# print the image
if ($param == 1) {
   $gif_data = $im->gif;
   open (GIFFILE,">image.gif") || die;
   print GIFFILE $gif_data;
   close GIFFILE;
} else {
   print "Content-type: image/gif\n\n";
   print $im->gif;
}

#----------- ndays subroutine -------------
# number of days from sdate to edate (yyyymmdd)
# begin day and end day are both counted
sub ndays {
   local ($sdate, $edate) = @_;
   $ndays = 0;
   # find start date day of year equivalent
   $day = substr($sdate, 6, 2);
   $mon = substr($sdate, 4, 2);
   $syr = substr($sdate, 0, 4);
   $sjday = &jday($mon, $day, $syr);
   # find end date day of year equivalent
   $day = substr($edate, 6, 2);
   $mon = substr($edate, 4, 2);
   $eyr = substr($edate, 0, 4);
   $ejday = &jday($mon, $day, $eyr);
#  print ("$sjday $ejday\n");
   if ($syr eq $eyr) {
      $ndays = $ejday - $sjday + 1;
   } else {
      $ndays = $ejday;
      for ($i=$syr; $i<$eyr; $i++ ) {
         if (($i/4) == (int($i/4))) {
            $ndays = $ndays + (366 - $sjday);
            $sjday = 0;
         } else {
            $ndays = $ndays + (365 - $sjday);
            $sjday = 0;
         }
      }
      $ndays = $ndays + 1;
   } 
}
 
#----------- jday subroutine --------------
sub jday {
   local ($mon, $day, $yr) = @_;
#  print ("$mon $day \n");
   # leap year (every fourth (2000 not applied))
   if (($yr/4) == (int($yr/4))) { 
      if ($mon == " 01") {
         $jday = $day;
      } elsif ($mon == " 02") {
         $jday = 31 + $day;
      } elsif ($mon == " 03") {
         $jday = 60 + $day;
      } elsif ($mon == " 04") {
         $jday = 91 + $day;
      } elsif ($mon == " 05") {
         $jday = 121 + $day;
      } elsif ($mon == " 06") {
         $jday = 152 + $day;
      } elsif ($mon == " 07") {
         $jday = 182 + $day;
      } elsif ($mon == " 08") {
         $jday = 213 + $day;
      } elsif ($mon == " 09") {
         $jday = 244 + $day;
      } elsif ($mon == " 10") {
         $jday = 274 + $day;
      } elsif ($mon == " 11") {
         $jday = 305 + $day;
      } elsif ($mon == " 12") {
         $jday = 335 + $day;
      } else {
         print ("julian day conversion error.\n");
      }
   } else {
      if ($mon == " 01") {
         $jday = $day;
      } elsif ($mon == " 02") {
         $jday = 31 + $day;
      } elsif ($mon == " 03") {
         $jday = 59 + $day;
      } elsif ($mon == " 04") {
         $jday = 90 + $day;
      } elsif ($mon == " 05") {
         $jday = 120 + $day;
      } elsif ($mon == " 06") {
         $jday = 151 + $day;
      } elsif ($mon == " 07") {
         $jday = 181 + $day;
      } elsif ($mon == " 08") {
         $jday = 212 + $day;
      } elsif ($mon == " 09") {
         $jday = 243 + $day;
      } elsif ($mon == " 10") {
         $jday = 273 + $day;
      } elsif ($mon == " 11") {
         $jday = 304 + $day;
      } elsif ($mon == " 12") {
         $jday = 334 + $day;
      } else {
         print ("julian day conversion error.\n");
      }
   }
}

#----------- incdate subroutine --------------
sub incdate {
   local ($dy, $mo, $yr, $inc) = @_;
   if ($inc > 31) {
      if (($inc/30) > 12) {
         $dy = "01";
         $mo = "01";
         $yr = $yr +1;
      } else {
         $dy = "01";
         $mo = $mo + int($inc/30);
         if ($mo > 12) {
            $yr = $yr + 1; 
            $mo = $mo - 12;
         }
      }
   } else {   
      $dy = $dy + $inc;
      if (($dy + ($inc/2)) > 31) { $dy = 32; }
      if ($dy >31) {
         $dy = "01";
         $mo = $mo + 1;
      } elsif (($dy > 30) && (($mo eq "04") || ($mo eq "06") || ($mo eq "09") || ($mo eq "11"))) {
         $dy = "01";
         $mo = $mo + 1;
      } elsif (($dy > 29) && ($mo eq "02") && (int($yr/4) == ($yr/4))) {
         $dy = "01";
         $mo = "03";
      } elsif (($dy > 28) && ($mo eq "02")) {
         $dy = "01";
         $mo = "03";
      }
   }
   if (length($dy) < 2) { $dy = "0" . $dy; }
   if (length($mo) < 2) { $mo = "0" . $mo; }
   if ($mo > 12) { 
      $yr = $yr + 1; 
      $mo = "01";
   }
   $incdate = $yr . $mo . $dy;
}