Listing 3. plotdays.pl
#!/usr/bin/perl
#
#
# plotdays.pl
#
# This guy uses the PGPLOT module, which in turn is based
# on the pgplot FORTRAN libraries... It's a long story, but
# it works. The goal is to produce some GIF format files
# of interesting statistics for the last X days.
#
use PGPLOT; # Load PGPLOT module
require "timelocal.pl";
$DATADIR = "/home/weather/data";
$twentyfourhours = 86400; #seconds in a day
@daysofweek = ('Sunday','Monday','Tuesday','Wednesday','Thursday','Friday','Saturday');
# Get the number of days of interest
if( $#ARGV >= 0 ) { $days=$ARGV[0]; }
else { $days = 7; }
$maxouttemp = 0; $minouttemp = 100;
$maxbp = 0; $minbp = 100;
$maxraind = 0; $minraind = 100;
$maxrainrate = 0; $minrainrate = 100;
$maxwinds = 0; $minwinds = 100;
# Loop through the number of days we are interested in.
# My loop index _decreases_, day X is further back
# in history than X-1.
#
# In true FORTRAN style, we are keeping parallel arrays for
# each statistic and label. And $j holds the total number
# of data points for each statistic's array.
$j = 0;
for($i=$days; $i >= 0; $i-- )
{
# Caculate file name for that day in the past.
$dtime = time();
$ttime = $dtime-($twentyfourhours * $i);
@time = localtime($ttime);
$DAYFILE = sprintf "%s/%02d%02d%02d",
$DATADIR, $time[4]+1, $time[3], $time[5];
$daylabel[$i] = sprintf "%s", $daysofweek[$time[6]];
$daydate[$i] = sprintf "%2d/%02d/%02d",$time[4]+1,$time[3],$time[5];
if( -r "$DAYFILE" )
{
# Open that day's file and read the data into local arrays
# We are interested in outtemp, humidity,
# pressure, rainfall, and windspeed.
open(IN, "$DAYFILE") or die "Cannot open file $DAYFILE";
while ()
{
($time, $date, $wdir, $winds, $aux, $intemp,
$outtemp, $hum, $bp, $raind, $rainm, $rainrate)
= split;
# Just in case a Min or Max slips in somehow
if( $time =~ /[M]/ ) { next; }
($hour,$min) = split(/:/,$time);
# Make basetime be the time in seconds at
# first of the day on that recent day in question
# We want to end up with a Unix-style long integer
# time value for each data point.
$time[0] = $time[1] = $time[2] = 0;
$basetime[$i] = timelocal(@time);
$today_time = (((60*$hour)+$min)*60);
$readingtime = $basetime[$i] + $today_time;
$readingtime = $readingtime - $basetime[$days];
# This is the X value for each plot
$xval[$j] = $readingtime;
$outtemp[$j] = $outtemp;
if( $outtemp > $maxouttemp ) { $maxouttemp = $outtemp; }
if( $outtemp < $minouttemp ) { $minouttemp = $outtemp; }
$dp[$j] = &dewpoint($outtemp,$hum);
$bp[$j] = $bp;
if( $bp > $maxbp ) { $maxbp = $bp; }
if( $bp < $minbp ) { $minbp = $bp; }
$raind[$j] = $raind;
if( $raind > $maxraind ) { $maxraind = $raind; }
if( $raind < $minraind ) { $minraind = $raind; }
$rainrate[$j] = $rainrate;
if( $rainrate > $maxrainrate )
{ $maxrainrate = $rainrate; }
if( $rainrate < $minrainrate )
{ $minrainrate = $rainrate; }
$winds[$j] = $winds;
if( $winds > $maxwinds ) { $maxwinds = $winds; }
if( $winds < $minwinds ) { $minwinds = $winds; }
$j++;
}
close IN;
}
}
# Then I just create the plots using a whole bunch
# of calls to pgplot routines. Not real pretty, but
# it gets the job done.
#
# Temp and dewpoint plot
#
$dev = "temp.gif/GIF";
$dev = "?" unless defined $dev;
pgbegin(0,$dev,1,1); # Open plot device
pgpap(7.0,0.55);
pgscr(1,0.75,0.75,0.75);
pgscr(5,0.40,0.70,0.40);
pgscr(0,1.0,1.0,1.0);
pgscf(1); # Set character font
pgslw(1); # Set line width
pgsch(1.2); # Set character height
pglabel("","",""); # Labels
$top = (int(($maxouttemp)/10)+1) * 10 ;
$bottom = (int(($minouttemp)/10)) * 10 ;
$left = 0;
$right = ($basetime[0] + $twentyfourhours) - $basetime[$days];
pgswin($left,$right,$bottom,$top);
pgbox( 'ABCGI',$twentyfourhours,0,'ABCGINMVS',10,2);
pgsci(0); # Change colour
pgdraw($xval[0],$outtemp[0]);
pgsci(2); # Change colour
for( $i=0; $i < $j; $i++)
{
pgdraw($xval[$i],$outtemp[$i]);
}
pgsci(2); # Change colour
pgptxt(0, $top+3.5 ,0,0,"Temperature (F)");
pgsci(4); # Change colour
pgptxt(0, $top+1 ,0,0,"Dewpoint (F)");
pgsci(5); # Change colour
pgptxt(int(($right-$left)/2), $top+3.5 ,0,0, "$time $daylabel[0] $date ");
pgsci(0); # Change colour
pgdraw($xval[0],$dp[0]);
pgsci(4); # Change colour
for( $i=0; $i < $j; $i++)
{
pgdraw($xval[$i],$dp[$i]);
}
pgsci(5); # Change colour
pgsch(1.20); # Set character height
for( $i=$days; $i >= 0; $i-- )
{
$offset = (($days-$i) * $twentyfourhours) + int($twentyfourhours/2);
pgptxt($offset, $bottom-2 ,0,0.5,$daylabel[$i]);
pgptxt($offset, $bottom-4 ,0,0.5,$daydate[$i]);
}
pgend; # Close plot
###
#
# Atmospheric Pressure Plot
#
$dev = "pressure.gif/GIF";
pgbegin(0,$dev,1,1); # Open plot device
pgpap(7.0,0.55);
pgscr(1,0.75,0.75,0.75);
pgscr(5,0.40,0.70,0.40);
pgscr(0,1.0,1.0,1.0);
pgscf(1); # Set character font
pgslw(1); # Set line width
pgsch(1.2); # Set character height
pglabel("","",""); # Labels
$top = (int($maxbp * 2) + 1) / 2;
$bottom = (int($minbp * 2)) / 2;
$left = 0;
$right = ($basetime[0] + $twentyfourhours) - $basetime[$days];
pgswin($left,$right,$bottom,$top);
pgbox( 'ABCGI',$twentyfourhours,0,'ABCGINMVS',.5,2);
pgsci(0); # Change colour
pgdraw($xval[0],$bp[0]);
pgsci(2); # Change colour
for( $i=0; $i < $j; $i++)
{
pgdraw($xval[$i],$bp[$i]);
}
pgsci(2); # Change colour
pgptxt(0, $top+.05 ,0,0,"Pressure (inches of mercury)");
pgsci(5); # Change colour
pgptxt(int(($right-$left)/2), $top+.05 ,0,0, "$time $daylabel[0] $date ");
pgsci(5); # Change colour
pgsch(1.20); # Set character height
for( $i=$days; $i >= 0; $i-- )
{
$offset = (($days-$i) * $twentyfourhours) + int($twentyfourhours/2);
pgptxt($offset, $bottom-.05 ,0,0.5,$daylabel[$i]);
pgptxt($offset, $bottom-.10 ,0,0.5,$daydate[$i]);
}
pgend; # Close plot
###
#
# Daily Rainfall and Rainfall Rate
#
$dev = "raind.gif/GIF";
pgbegin(0,$dev,1,1); # Open plot device
pgpap(7.0,0.55);
pgscr(1,0.75,0.75,0.75);
pgscr(5,0.40,0.70,0.40);
pgscr(0,1.0,1.0,1.0);
pgscf(1); # Set character font
pgslw(1); # Set line width
pgsch(1.2); # Set character height
pglabel("","",""); # Labels
if( $maxraind > $maxrainrate ) { $maxrain = $maxraind; }
else { $maxrain = $maxraind; }
if( $minraind > $minrainrate ) { $minrain = $minraind; }
else { $minrain = $minraind; }
$top = (int($maxrain * 2) + 1) / 2;
$bottom = (int($minrain * 2)) / 2;
$left = 0;
$right = ($basetime[0] + $twentyfourhours) - $basetime[$days];
pgswin($left,$right,$bottom,$top);
pgbox( 'ABCGI',$twentyfourhours,0,'ABCGINMVS',.5,2);
pgsci(0); # Change colour
pgdraw($xval[0],$raind[0]);
pgsci(2); # Change colour
for( $i=0; $i < $j; $i++)
{
pgdraw($xval[$i],$raind[$i]);
}
pgsci(4); # Change colour
for( $i=0; $i < $j; $i++)
{
# Just print positive data points for rainrate
if( $rainrate[$i] > 0 )
{
$y[0] = $rainrate[$i];
$x[0] = $xval[$i];
pgpt(1,\@x,\@y,20);
}
}
pgsci(2); # Change colour
pgptxt(0, $top+.08 ,0,0,"Daily Rainfall (inches)");
pgsci(4); # Change colour
pgptxt(0, $top+.04 ,0,0,"Rainfall Rate (inches per hour)");
pgsci(5); # Change colour
pgptxt(int(($right-$left)/2), $top+.05 ,0,0, "$time $daylabel[0] $date ");
pgsci(5); # Change colour
pgsch(1.20); # Set character height
for( $i=$days; $i >= 0; $i-- )
{
$offset = (($days-$i) * $twentyfourhours) + int($twentyfourhours/2);
pgptxt($offset, $bottom-.05 ,0,0.5,$daylabel[$i]);
pgptxt($offset, $bottom-.10 ,0,0.5,$daydate[$i]);
}
pgend; # Close plot
###
#
# Windspeed Plot
#
$dev = "winds.gif/GIF";
pgbegin(0,$dev,1,1); # Open plot device
pgpap(7.0,0.55);
pgscr(1,0.75,0.75,0.75);
pgscr(5,0.40,0.70,0.40);
pgscr(0,1.0,1.0,1.0);
pgscf(1); # Set character font
pgslw(1); # Set line width
pgsch(1.2); # Set character height
pglabel("","",""); # Labels
$top = (int(($maxwinds)/5)+1) * 5 ;
$bottom = (int(($minwinds)/5)) * 5 ;
$left = 0;
$right = ($basetime[0] + $twentyfourhours) - $basetime[$days];
pgswin($left,$right,$bottom,$top);
pgbox( 'ABCGI',$twentyfourhours,0,'ABCGINMVS',5,0);
pgsci(0); # Change colour
pgdraw($xval[0],$winds[0]);
pgsci(2); # Change colour
for( $i=0; $i < $j; $i++)
{
pgdraw($xval[$i],$winds[$i]);
}
pgsci(2); # Change colour
pgptxt(0, $top+1 ,0,0,"Windspeed (MPH)");
pgsci(5); # Change colour
pgptxt(int(($right-$left)/2), $top+1 ,0,0, "$time $daylabel[0] $date ");
pgsci(5); # Change colour
pgsch(1.20); # Set character height
for( $i=$days; $i >= 0; $i-- )
{
$offset = (($days-$i) * $twentyfourhours) + int($twentyfourhours/2);
pgptxt($offset, $bottom-.70 ,0,0.5,$daylabel[$i]);
pgptxt($offset, $bottom-1.4 ,0,0.5,$daydate[$i]);
}
pgend; # Close plot
exit;
# Subroutines to do things like calculate dewpoint
# and corrected atmospheric pressure.
sub cbp
{
local($in) = @_;
return(sprintf("%2.2f", $in + $CBP_CORRECTION));
}
sub ctof
{
local ($C) = @_;
$F = (9/5 * $C) + 32;
return $F;
}
sub ftoc
{
local ($F) = @_;
$C = 5/9 * ($F - 32);
return $C;
}
sub dewpoint
{
local ($F, $H) = @_;
$rh = $H/100; # RH. Fraction, [0..1] not %.
$t = ftoc($F); #Celsius Temperature
$es = 6.112 * exp(17.67 * $t / ($t + 243.5));
$e = $rh * $es;
$loge = log($e/6.112);
$dp = 243.5 * $loge/(17.67 - $loge);
return ctof($dp);
}
