#!/usr/bin/perl -w

# Norman Kuring		25-Jun-2010	Original development
# Norman Kuring		13-Aug-2010	Emphasize grid lines corresponding to
#					time annotations.
# Norman Kuring         17-Nov-2010	Increase output image height to
#					accomodate legends containing many
#					parameters.
# Norman Kuring		27-Aug-2010	Plot parameters down to the fractional
#					second and add 1-second gridlines to
#					plots that span a minute or less.
# Norman Kuring		 9-Sep-2010	Allow the user to specify the minimum
#					and maximum values plotted.  Allow the
#					y-scale for each units grouping to be
#					constrained independently.
# Norman Kuring		21-Apr-2011	Allow comments in parameter index files.
# Norman Kuring		 6-May-2011	Allow plotting of SAC-D service
#					platform telemetry as well as the
#					Aquarius instrument telemetry.
# Norman Kuring		23-Jun-2011	Implement leap second correction when
#					converting GPS time to UTC.
# Norman Kuring          1-Jul-2011     Large data values cause a problem.
#					It turns out that the problem is
#					twofold. First, all of the values
#					for this parameter in the recently
#					corrupted data are equal to
#					39251829815707348593623868279812521984 
#					. To get the minimum and maximum
#					value for my plot I use the
#					maximum and minimum for the
#					data. When these are equal, as in
#					this case, I add and subtract 0.5
#					to give the plot a finite height.
#					This leads to the second part of
#					the problem. Adding 0.5 to such
#					a large number is no different
#					from adding zero in the discrete
#					world of computer arithmetic, so
#					I end up requesting a plot with
#					zero height which fails.
#					This has been fixed.
# Norman Kuring		11-Jul-2011	Read condensed versions of the data for
#					longer time periods.
# Norman Kuring		12-Jul-2011	Make plot title smaller and allow it to
#					wrap to accommodate some of the longer
#					descriptive strings in the metadata.
# Norman Kuring		13-Jul-2011	Handle a no-data condition a little
#					more gracefully.
# Norman Kuring		20-Jul-2011	Look back far enough to include 3-day
#					stored telemetry files from CONAE.
# Norman Kuring		25-Aug-2011	If the min and max value of a time
#					series were the same negative value,
#					the plot would fail. Fixed.
# Norman Kuring		18-Nov-2011	Use modular versions of read_units and
#					get_param_specs subroutines.
# Norman Kuring		11-Jan-2012	Add capability to read condensed daily
#					data.
# Norman Kuring		12-Jan-2012	Move read_condensed subroutine to
#					ReadTelemMeta module and fix a plot
#					boundary truncation problem.
# Norman Kuring		13-Jun-2012	Plot the mean instead of the median
#					for aggregated statistics.

use strict;
use lib '/web/aquarius/scripts';
use ReadTelemMeta qw(get_param_specs read_units read_condensed);
use GpsToUtc qw(&read_leap_secs &gps2utc);
use Utility qw(isleap);
use Getopt::Long;
use Time::Local;
use POSIX;

$ENV{PATH} = '/bin:/usr/local/gmt/bin';

my $psbasemap = '/usr/local/gmt/bin/psbasemap';
my $psxy      = '/usr/local/gmt/bin/psxy';
my $pslegend  = '/usr/local/gmt/bin/pslegend';
my $pstext    = '/usr/local/gmt/bin/pstext';
my $gs        = '/usr/bin/gs';

# Aggregation thresholds
my $aggthresh1 =     86_400;	# 1 day
my $aggthresh2 =  6_912_000;	# 80 days

my @colors = qw(
0/175/0    140/0/70  0/0/255   128/0/150
0/80/0     50/30/0  0/130/255 255/0/255
128/50/0  255/50/0  80/0/255  230/110/0
140/100/60 80/80/0  0/102/102 160/60/60
);

# Make lighter, reduced-saturation versions of the above colors
# for the plotting of maxima and minima.  I do this in hopes that
# the more muted colors will not distract from the median values
# when plotting aggregated data sets.
my @muted =
map {
  my ($r,$g,$b) = split /\//;
  my ($h,$s,$v) = rgb2hsv($r,$g,$b);
  $v += 0.7*(1 - $v);
  $s *= 0.3;
  join '/',hsv2rgb($h,$s,$v);
} @colors;

# Use these symbols
#    a - star      i - inverted triangle
#    c - circle    n - pentagon
#    d - diamond   s - square
#    g - octagon   t - triangle
#    h - hexagon   x - cross
my @symbols = qw(c a d g h i n s t x);

# I use the following hash to determine the width of the plot legend.
# This only works as long as I use a 10-point Helvetica-Bold font in
# my pslegend call.
#
# I generated the following numbers as follows.
#   perl -e '$q=chr 39;for($i=33;$i<127;$i++){' \
#        -e '$c=chr $i;' \
#        -e '$cmd="pstext -R0/150/-0.5/0.5 -JX150/1 -P -X0 -Y0"' \
#        -e '   . "| gs -sDEVICE=pgmraw -sOutputFile=- -q -r100"' \
#        -e '   . " -g15000x100 -dTextAlphaBits=4 - | pnmcrop>nak.pgm";' \
#        -e 'open P,"|-",$cmd or die;' \
#        -e 'print P "0 0 10 0 1 BL ",$c x 1000,"\n";' \
#        -e 'close P or die;' \
#        -e '$cmd="pnmfile nak.pgm";' \
#        -e 'open P,"-|",$cmd or die;' \
#        -e '$w=(split " ",<P>)[3]/1000;' \
#        -e 'close P or die;' \
#        -e 'printf "$q$c$q => %6.3f, ",$w}' \
#        -e 'print "\n"'
# I then manually modified the value for the at-sign
# and added a value for the space and the degree symbol.
my $degcode = chr 217;
my %cwid = (
'!' =>  4.657, '"' =>  6.632, '#' =>  7.781, '$' =>  7.781, '%' => 12.444,
'&' => 10.104, "'" =>  3.889, '(' =>  4.660, ')' =>  4.659, '*' =>  5.445,
'+' =>  8.176, ',' =>  3.889, '-' =>  4.660, '.' =>  3.889, '/' =>  3.891,
'0' =>  7.780, '1' =>  7.777, '2' =>  7.780, '3' =>  7.780, '4' =>  7.781,
'5' =>  7.780, '6' =>  7.780, '7' =>  7.781, '8' =>  7.780, '9' =>  7.780,
':' =>  4.657, ';' =>  4.657, '<' =>  8.176, '=' =>  8.176, '>' =>  8.176,
'?' =>  8.549, '@' => 13.646, 'A' => 10.105, 'B' => 10.103, 'C' => 10.104,
'D' => 10.103, 'E' =>  9.335, 'F' =>  8.550, 'G' => 10.889, 'H' => 10.103,
'I' =>  3.889, 'J' =>  7.780, 'K' => 10.104, 'L' =>  8.549, 'M' => 11.658,
'N' => 10.103, 'O' => 10.889, 'P' =>  9.335, 'Q' => 10.889, 'R' => 10.104,
'S' =>  9.336, 'T' =>  8.551, 'U' => 10.103, 'V' =>  9.336, 'W' => 13.215,
'X' =>  9.336, 'Y' =>  9.336, 'Z' =>  8.550, '[' =>  4.659, '\\'=>  1.945,
']' =>  4.659, '^' =>  8.175, '_' =>  7.781, '`' =>  3.889, 'a' =>  7.781,
'b' =>  8.549, 'c' =>  7.781, 'd' =>  8.550, 'e' =>  7.781, 'f' =>  4.660,
'g' =>  8.550, 'h' =>  8.549, 'i' =>  3.889, 'j' =>  3.890, 'k' =>  7.780,
'l' =>  3.889, 'm' => 12.444, 'n' =>  8.549, 'o' =>  8.550, 'p' =>  8.549,
'q' =>  8.550, 'r' =>  5.445, 's' =>  7.780, 't' =>  4.660, 'u' =>  8.549,
'v' =>  7.781, 'w' => 10.891, 'x' =>  7.781, 'y' =>  7.781, 'z' =>  7.000,
'{' =>  5.445, '|' =>  3.915, '}' =>  5.444, '~' =>  8.176, ' ' =>  3.888,
"$degcode"=>  8.440 );

my %units = &read_units;

# Make a table for the usage message.
my $unit_table = "Unit key   Meaning\n" . '-' x 8 . '   ' . '-' x 36 . "\n";
foreach my $u (sort { $units{$a}[0] cmp $units{$b}[0] } keys %units){
  $units{$u}[0] or next;
  $unit_table .= sprintf "%-8s   %s\n",$u,$units{$u}[0];
}

my $plot_wid = 10;
my $plot_hgt = 6;

my ($min, $max, $show_alarms, $show_lines, $abscissa);
GetOptions(
  'min=f%'     => \$min,
  'max=f%'     => \$max,
  'alarms'     => \$show_alarms,
  'lines'      => \$show_lines,
  'abscissa=s' => \$abscissa,
);

my ($start_time,$end_time,@params) = @ARGV;

die <<"--" unless @params;
Usage:
$0 [options] start_time end_time parameter_id [parameter_id [parameter_id ...]]
The start_time and end_time arguments should be specified as
YYYY-MM-DDThh:mm:ss (e.g. 2010-05-19T14:27:31 for May 19, 2010
at 2:27:31 PM) with times given in UTC.  If you leave parts off
of the end of the string, start times will be the beginning of
the specified period and end times will get set to the end of the
specified period.  For example: a start time of 2010-05 equates to
2010-05-01T00:00:00; an end time of 2010-05-19T14 equates to
2010-05-19T14:59:59.

Options:
  -min key=value	Use value as the minimum value to be plotted
			instead of the minimum value of the data; key
			indicates which unit the minimum applies to
			(see the table below for available keys).
  -max key=value	Use value as the maximum value to be plotted
			instead of the maximum value of the data; key
			indicates which unit the maximum applies to
			(see the table below for available keys).
  -alarms		Highlight the yellow and red alarm thresholds
			of those parameters that have them
  -lines		Draw lines connecting the plotted points of
			a time series
  -abscissa parameter	Use the specified parameter as the abscissa
			of the plot

$unit_table
--

my $st = parse_time($start_time,0);
my $et = parse_time($end_time,1);

my($bp,$bs,$ai,$gi,$gd) = time_axis($et - $st);

# Get parameter information as a hash of arrays.
#                          [    0   ,     1      ,  2  , 3, 4, 5, 6,    7      ]
# $pindex{TLMTYP}{PARAM_ID}[GROUP_ID,GROUP_COLUMN,UNITS,RL,YL,YH,RH,DESCRIPTION]
my %pindex;
foreach my $telemtyp (qw/AQ SP/){
  %{$pindex{$telemtyp}} = get_param_specs($telemtyp);
}

# Initialize the GPS to UTC conversion function.
&read_leap_secs;

my @timeseries;
my %minmax;
my %unitcount;
foreach my $p (@params){
  my $telemtyp = telemtyp($p);
  if(not defined $telemtyp){
    $p = undef;
    next;
  }

  my ($st,$et,@ts) =
  read_param($telemtyp,@{ $pindex{$telemtyp}{$p} }[0,1],$st,$et);

  push @timeseries, [ @ts ];
  if( @ts ){
    my $units = $pindex{$telemtyp}{$p}[2];
    $unitcount{$units}++;
    push @{ $minmax{$units} }, [ minmax( \@ts ) ];
  }

  # The read_param subroutine might modify the $st and $et values,
  # so I turn them back to strings here.
  $start_time = tstring($st);
  $end_time   = tstring($et);
}
# Remove any parameters that were not found in the metadata.
for(my $i=$#params; $i>=0; $i--){
  splice(@params,$i,1) if not defined $params[$i];
}

my (%min, %max, %inc);
foreach my $units (keys %unitcount){
  $min{$units} = defined $min->{$units} ?
                 $min->{$units}         :
                 ( sort {$a <=> $b} map $_->[0], @{ $minmax{$units} } )[0];
  $max{$units} = defined $max->{$units} ?
                 $max->{$units}         :
                 ( sort {$b <=> $a} map $_->[1], @{ $minmax{$units} } )[0];
  my $range = $max{$units} - $min{$units};
  if($range > 0){
    $inc{$units} = nice_increment($range);
    $min{$units} -= 0.02*$range unless defined $min->{$units};
    $max{$units} += 0.02*$range unless defined $max->{$units};
  }
  elsif($range < 0){
    if(defined $min->{$units}){
      $max{$units} = $min->{$units} + 1;
    }
    elsif(defined $max->{$units}){
      $min{$units} = $max->{$units} - 1;
    }
    $inc{$units} = 0.1;
  }
  else{
    my $factor = $min{$units} == 0 ? 1 : abs($min{$units});
    $inc{$units}  = 0.1*$factor;
    $min{$units} -= 0.5*$factor;
    $max{$units} += 0.5*$factor;
  }
}
my $unitcount = keys %unitcount;

my $top_unit;
if($unitcount == 0){
  $top_unit = $pindex{AQ}{$params[0]}    ?
              $pindex{AQ}{$params[0]}[2] :
              $pindex{SP}{$params[0]}[2];
  $unitcount++;
  $unitcount{$top_unit} = 0;
  $min{$top_unit} = defined $min->{$top_unit} ?
                    $min->{$top_unit}         :
                    defined $max->{$top_unit} ?
                    $max->{$top_unit} - 1     :
                    0;
  $max{$top_unit} = defined $max->{$top_unit} ?
                    $max->{$top_unit}         :
                    defined $min->{$top_unit} ?
                    $min->{$top_unit} + 1     :
                    1;
  $inc{$top_unit} = nice_increment($max{$top_unit} - $min{$top_unit});
}
else{
  $top_unit = (sort { $unitcount{$b} <=> $unitcount{$a} } keys %unitcount)[0];
}

# Make a temporary working directory and go there.
my $tmpdir = '/tmpfs01/' . (split /\//,$0)[-1] . ".$^T.$$";
mkdir $tmpdir,0755 or die "Could not make $tmpdir ($!)";
chdir $tmpdir or die "Could not cd to $tmpdir ($!)";

# Draw the initial basemap frame.
doit( qq(
$psbasemap
-R$start_time/$end_time/$min{$top_unit}/$max{$top_unit}
-JX${plot_wid}T/$plot_hgt
-P -X$unitcount -Y1
-Bp$bp/0Sn
-Bs$bs/0
$gd
--LABEL_OFFSET=0.01i
--ANNOT_FONT_PRIMARY=Helvetica
--ANNOT_FONT_SIZE_PRIMARY=14p
--ANNOT_OFFSET_PRIMARY=0.075i
--BASEMAP_FRAME_RGB=+0/0/0
-K
>basemap.ps
) );

# Draw the Y-axis separately so I can set font separately.
doit( qq(
$psbasemap
-R$start_time/$end_time/$min{$top_unit}/$max{$top_unit}
-JX${plot_wid}T/$plot_hgt -P
-Bp0/$inc{$top_unit}:"$units{$top_unit}[1]":We
--ANNOT_FONT_PRIMARY=Helvetica-Bold
--ANNOT_FONT_SIZE_PRIMARY=10p
--ANNOT_OFFSET_PRIMARY=0.075i
--LABEL_FONT=Helvetica-Bold
--LABEL_FONT_SIZE=16p
--LABEL_OFFSET=0.01i
--BASEMAP_FRAME_RGB=+0/0/0
-O -K
>>basemap.ps
) );

# Keep track of overall dimensions.
my $final_wid = $plot_wid + $unitcount + 0.2;
my $final_hgt = $plot_hgt + 2;

if( @timeseries > 1 ){
  # Open a file in which to write legend information.
  open F,'>',"legendspec.txt" or die "Could not open legendspec.txt ($!)";
}

my $legendrows = 0;
my $sym = 0;
my $clr = 0;
my $maxwid = 0;

for(my $i=0; $i<@params; $i++){

  if( @{ $timeseries[$i] } ){

    my $telemtyp = telemtyp($params[$i]);

    my $units = $pindex{$telemtyp}{$params[$i]}[2];

    if($show_alarms){
      # Highlight the alarm ranges for this parameter if they are available.
      my($rl,$yl,$yh,$rh) = @{ $pindex{$telemtyp}{$params[$i]} }[3..6];

      my $do_rl = $rl =~ /\d/ && $min{$units} <= $rl;
      my $do_rh = $rh =~ /\d/ && $max{$units} >= $rh;
      my $do_yl = $yl =~ /\d/ && $min{$units} <= $yl;
      my $do_yh = $yh =~ /\d/ && $max{$units} >= $yh;

      if($do_rl or $do_rh or $do_yl or $do_yh){

        my $cmd = join ' ',split ' ',<<"        --";
        $psxy
        -R$start_time/$end_time/$min{$units}/$max{$units}
        -JX${plot_wid}T/$plot_hgt
        -P -M -O -K
        >>basemap.ps
        --

        open P,'|-',$cmd or die "Error opening pipe to $cmd ($!)";

        my $yellow = '255/255/80';
        my $red    = '255/150/150';	# More pink really, but I want a light
					# enough color not to overwhelm the
					# plotted data.

        if($do_yl){
          print P "> -G$yellow\n";
          print P "$start_time $min{$units}\n";
          print P "$end_time $min{$units}\n";
          print P "$end_time $yl\n";
          print P "$start_time $yl\n";
        }
        if($do_yh){
          print P "> -G$yellow\n";
          print P "$start_time $max{$units}\n";
          print P "$end_time $max{$units}\n";
          print P "$end_time $yh\n";
          print P "$start_time $yh\n";
        }
        if($do_rl){
          print P "> -G$red\n";
          print P "$start_time $min{$units}\n";
          print P "$end_time $min{$units}\n";
          print P "$end_time $rl\n";
          print P "$start_time $rl\n";
        }
        if($do_rh){
          print P "> -G$red\n";
          print P "$start_time $max{$units}\n";
          print P "$end_time $max{$units}\n";
          print P "$end_time $rh\n";
          print P "$start_time $rh\n";
        }
        close P or die "Error executing $cmd ($!)";
      }
    }
  }
}
my %telemtyp;
my $span = $et - $st;
my $iscondensed = $span > $aggthresh1 ? 1 : 0;

if($iscondensed){
  # Make custom symbols to plot minima and maxima.
  open CS,'>',"min.def" or die "Could not create min.def ($!)";
  print CS "-0.5 0.32087 M\n",
          "0 0.55 1.1 204.62 335.38 A\n",
          "0 0.55 1.1 335.38 204.62 A\n";
  close CS;
  open CS,'>',"max.def" or die "Could not create max.def ($!)";
  print CS "0.5 -0.32087 M\n",
          "0 -0.55 1.1 24.62 155.38 A\n",
          "0 -0.55 1.1 155.38 24.62 A\n";
  close CS;
}

my $data_plotted = 0;
for(my $i=0; $i<@params; $i++){

  if( @{ $timeseries[$i] } ){

    my $v = $iscondensed ? 2 : 1;

    my $telemtyp = telemtyp($params[$i]);

    # Keep track of the telemetry types so I can assign
    # an appropriate title to multi-parameter plots.
    $telemtyp{$telemtyp}++;

    my $units = $pindex{$telemtyp}{$params[$i]}[2];

    if($iscondensed){

      # Plot the minimum-to-maximum range for each aggregate statistic.

      my $delta_t = $span > $aggthresh2 ? 86_400 : 300;
      my $delt = $delta_t/2;

      my $cmd = join ' ', split ' ', <<"      --";
      $psxy
      -R$start_time/$end_time/$min{$units}/$max{$units}
      -JX${plot_wid}T/$plot_hgt
      -P -G180 -M -O -K      >>basemap.ps
      --

      open P,'|-',$cmd or die "Error opening pipe to $cmd ($!)";
      foreach my $p (@{ $timeseries[$i] }){
        print P ">\n",
        tstring($p->[0]-$delt)," $p->[4]\n",
        tstring($p->[0]-$delt)," $p->[5]\n",
        tstring($p->[0]+$delt)," $p->[5]\n",
        tstring($p->[0]+$delt)," $p->[4]\n";
      }
      close P or die "Error executing pipe $cmd ($!)";

      # Cap each range with the custom symbols defined above and give
      # them the same color (but muted) as the corresponding plotted parameter.
      $cmd = join ' ', split ' ', <<"      --";
      $psxy
      -R$start_time/$end_time/$min{$units}/$max{$units}
      -JX${plot_wid}T/$plot_hgt
      -P -S0.08 -W3,$muted[$clr] -O -K >>basemap1.ps
      --

      open P,'|-',$cmd or die "Error opening pipe to $cmd ($!)";
      foreach my $p (@{ $timeseries[$i] }){
        print P
        tstring($p->[0])," $p->[4] kmin\n",
        tstring($p->[0])," $p->[5] kmax\n";
      }
      close P or die "Error executing pipe $cmd ($!)";

      # Put a label on the plot indicating that
      # it was generated from aggregated data.
      $cmd = join ' ', split ' ', <<"      --";
      $pstext
      -R0/$plot_wid/0/$plot_hgt -JX$plot_wid/$plot_hgt -P -N
      -G0/0/255 -O -K >>basemap1.ps
      --

      my $ai = $delta_t==300 ? "5-minute" : "daily";	# Aggregating interval

      open P,'|-',$cmd or die "Error opening pipe to $cmd ($!)";
      print P "$plot_wid @{[$plot_hgt+0.15]} 12 0 21 BR ",
              "Means, minima, and maxima plotted ",
              "from data aggregated over $ai intervals\n";
      close P or die "Error executing pipe $cmd ($!)";
      
    }

    # Plot the parameter time series.
    if($show_lines){
      my $cmd = join ' ', split ' ', <<"      --";
      $psxy
      -R$start_time/$end_time/$min{$units}/$max{$units}
      -JX${plot_wid}T/$plot_hgt
      -P -W5,$colors[$clr] -O -K
      >>basemap2.ps
      --

      open P,'|-',$cmd or die "Error opening pipe to $cmd ($!)";
      foreach my $p (@{ $timeseries[$i] }){
        print P tstring($p->[0])," $p->[$v]\n";
      }
      close P or die "Error executing pipe $cmd ($!)";
    }

    my $cmd = join ' ', split ' ', <<"    --";
    $psxy
    -R$start_time/$end_time/$min{$units}/$max{$units}
    -JX${plot_wid}T/$plot_hgt
    -P -S$symbols[$sym]0.08 -W3,$colors[$clr] -O -K
    >>basemap2.ps
    --

    open P,'|-',$cmd or die "Error opening pipe to $cmd ($!)";
    foreach my $p (@{ $timeseries[$i] }){
      print P tstring($p->[0])," $p->[$v]\n";
    }
    close P or die "Error executing pipe $cmd ($!)";

    $data_plotted = 1;

    if( @timeseries > 1){

      my $ustr = $unitcount > 1 ? " ($units{$units}[1])" : '';
      my $legendtxt = "$params[$i]$ustr";

      # Get the width of the legendtxt in pixels.
      my $wid = 0;
      my $txt = $legendtxt;
      $txt =~ s/\\217/$degcode/g;
      grep { $wid += $cwid{$_}; 0 } split //, $txt;

      $maxwid = $wid if $wid > $maxwid;

      print F
      "S 0.1i $symbols[$sym] 0.15i - 2p,$colors[$clr] 0.24i $legendtxt\n";

      $legendrows++;
    }
  }

  $sym = ($sym + 1)%@symbols;
  $clr = ($clr + 1)%@colors;
}

if($data_plotted){
  # Tack the plotted points onto the rest of the graph, so far.
  # I do it this way so that all of the plotted points lie on top
  # of all of the plotted minimum-to-maximum ranges when condensed
  # data are being displayed.
  my $bm1 = $iscondensed ? 'basemap1.ps' : '';
  doit( qq(/bin/cat $bm1 basemap2.ps >>basemap.ps) );
}

# Overlay a grid on the plot.
doit( qq(
$psbasemap
-R$start_time/$end_time/$min{$top_unit}/$max{$top_unit}
-JX${plot_wid}T/$plot_hgt
-P -Bg$gi/g$inc{$top_unit} -O -K
--GRID_CROSS_SIZE=0
>>basemap.ps
) );

if($ai ne $gi){
  # Emphasize annotation grid lines on the time axis.
  doit( qq(
  $psbasemap
  -R$start_time/$end_time/$min{$top_unit}/$max{$top_unit}
  -JX${plot_wid}T/$plot_hgt
  -P -Bg$ai/0 -O -K
  --GRID_CROSS_SIZE=0
  >>basemap.ps
  ) );
}


if($unitcount > 1){
  # Draw additional Y-axes on the plot if parameters having more than
  # one set of units are plotted.
  foreach my $unit (keys %unitcount){
    next if $unit eq $top_unit;		# This one's already drawn.
    doit( qq(
    $psbasemap
    -R0/1/$min{$unit}/$max{$unit}
    -JX1/$plot_hgt
    -P -B0/$inc{$unit}:"$units{$unit}[1]":W -X-1 -O -K
    --ANNOT_FONT_SIZE_PRIMARY=8p
    --ANNOT_OFFSET_PRIMARY=0.075i
    --LABEL_FONT=Helvetica-Bold
    --LABEL_FONT_SIZE=14p
    --LABEL_OFFSET=0.01i
    >>basemap.ps
    ) );
  }
}

# Draw a legend on the plot if needed.
my $origin = '-X-1 -Y-1';
if( @timeseries > 1 ){
  close F;

  if($legendrows){

    my $wid = 0.4 + $maxwid/100;
    # The following height is based on the 10-point ANNOT_FONT_SIZE_PRIMARY
    # and the 6-point line-spacing (-L1.6).
    my $hgt = $legendrows*10/72 + ($legendrows + 1)*6/72;
    $final_wid += 0.25 + $wid;
    my $x = $final_wid - 0.25;
    my $y = $final_hgt > $hgt ? ($final_hgt + $hgt)/2 : $hgt;
    $final_hgt = $hgt + 6/72 if $hgt >= $final_hgt;

    doit( qq(
    $pslegend legendspec.txt
    -R0/$final_wid/0/$final_hgt
    -JX$final_wid/$final_hgt -P $origin
    -Dx$x/$y/$wid/$hgt/TR -F
    --ANNOT_FONT_PRIMARY=Helvetica-Bold
    --ANNOT_FONT_SIZE_PRIMARY=10p
    -L1.6 -O -K
    >>basemap.ps
    ) );

    $origin = '';
  }
}

# Put a title on the plot and close the PostScript.
my $title;
if( @params == 1 ){
  my $telemtyp = (each %telemtyp)[0];
  $title = $pindex{$telemtyp}{$params[0]}[7];
}
else{
  my %longname = (SP => 'Service Platform', AQ => 'Aquarius');
  $title = join(' and ',map $longname{$_}, keys %telemtyp);
  $title = $title ? "$title Telemetry" : 'No Data Found';
}
my $title_x = $final_wid/2;
my $title_y = $final_hgt - 0.2;
doit( qq(
/bin/echo -e "> $title_x $title_y 12 0 1 CT 12p $final_wid c\\n$title"
| $pstext -R0/$final_wid/0/$final_hgt
  -JX$final_wid/$final_hgt -P -m -O $origin
>> basemap.ps
) );

$final_wid = int(100*$final_wid);
$final_hgt = int(100*$final_hgt);

doit( qq(
$gs -sDEVICE=pngalpha -sOutputFile=- -q -r100
-dGraphicsAlphaBits=4 -dTextAlphaBits=4 -g${final_wid}x$final_hgt -
<basemap.ps >basemap.png
) );

open F,'<',"basemap.png" or die "Could not open $tmpdir/basemap.png ($!)";
local($/) = undef;
my $png = <F>;
close F;

print $png;

chdir;
doit("/bin/rm -rf $tmpdir");

sub parse_time{
  my ($tstr,$wantend) = @_;

  my ($y,$mo,$d,$h,$m,$s);

  if(
    $tstr =~ /
    ^
    (\d{4})		# $1: year
    (?:
     \D 
      (\d\d?)		# $2: month		(1-12)
      (?:
        \D
        (\d\d?)		# $3: day of the month	(1-31)
        (?:
          \D
          (\d\d?)	# $4: hour		(0-23)
          (?:
            \D
            (\d\d?)	# $5: minute		(0-59)
            (?:
              \D
              (\d\d?)	# $6: second		(0-59)
            )?
          )?
        )?
      )?
    )?
    $
    /x
  ){
    $y = $1;
    if(defined $2){
      $mo = $2;
      if(defined $3){
        $d = $3;
        if(defined $4){
          $h = $4;
          if(defined $5){
            $m = $5;
            if(defined $6){
              $s = $6;
            }
            else{
              $s = $wantend ? 59 : 0;
            }
          }
          else{
            ($m,$s) = $wantend ? (59,59) : (0,0);
          }
        }
        else{
          ($h,$m,$s) = $wantend ? (23,59,59) : (0,0,0);
        }
      }
      else{
        ($d,$h,$m,$s) = $wantend ? (31,23,59,59) : (1,0,0,0);
      }
    }
    else{
      ($mo,$d,$h,$m,$s) = $wantend ? (12,31,23,59,59) : (1,1,0,0,0);
    }
  }
  else{
    die "Unrecognized time string ($tstr)\n";
  }
  timegm($s,$m,$h,$d,$mo-1,$y-1900);
}

sub tstring{
  my $t = shift;
  my $fs = sprintf "%f",$t - int($t);	# fractional second
  $fs =~ s/^0//;
  $fs == 0 and $fs = '';
  my ($s,$m,$h,$d,$mo,$y) = (gmtime $t)[0..5];
  sprintf "%04d-%02d-%02dT%02d:%02d:%02d$fs",$y+1900,$mo+1,$d,$h,$m,$s;
}

sub time_axis{
  my $timespan = shift;

  my ($bp,$bs);	# Primary and secondary boundary annotations (GMT -B option)
  my $gi;	# Interval between Vertical grid lines
  my $ai;	# Annotation interval
  my %gd = (	# GMT default parameter settings
    ANNOT_FONT_PRIMARY        => 'Helvetica',
    ANNOT_FONT_SECONDARY      => 'Helvetica',
    ANNOT_FONT_SIZE_PRIMARY   => '16p',
    ANNOT_FONT_SIZE_SECONDARY => '18p',
  );

  # The GMT settings in this function are for a plot that is 10 inches
  # (1000 pixels at a 100dpi screen resolution) wide, i.e. -JX10 .

  if   ($timespan <= 12){		# 12 seconds
    $bp = '1Cf1c';		$bs = '1D';	$gi = '1c';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm:ss';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 60){		# one minute
    $bp = '10Cf1c';		$bs = '1D';	$gi = '1c';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm:ss';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 90){		# one and a half minutes
    $bp = '10Cf1c';		$bs = '1D';	$gi = '10c';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm:ss';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 180){		# three minutes
    $bp = '20Cf1c';		$bs = '1D';	$gi = '10c';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm:ss';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 300){		# five minutes
    $bp = '30Cf2c';		$bs = '1D';	$gi = '10c';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm:ss';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 900){		# fifteen minutes
    $bp = '1Mf10c';		$bs = '1D';	$gi = '10c';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 1800){		# half an hour
    $bp = '2Mf10c';		$bs = '1D';	$gi = '1m';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 3600){		# one hour
    $bp = '5Mf20c';		$bs = '1D';	$gi = '1m';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 7200){		# two hours
    $bp = '10Mf1m';		$bs = '1D';	$gi = '2m';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 14_400){		# four hours
    $bp = '20Mf1m';		$bs = '1D';	$gi = '5m';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 21_600){		# six hours
    $bp = '30Mf2m';		$bs = '1D';	$gi = '10m';
    $gd{PLOT_CLOCK_FORMAT} = '-hh:mm';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 86400){		# one day
    $bp = '1Hf5m';		$bs = '1D';	$gi = '20m';
    $gd{PLOT_CLOCK_FORMAT} = '-hh';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 172_800){		# two days
    $bp = '2Hf30m';		$bs = '1D';	$gi = '1h';
    $gd{PLOT_CLOCK_FORMAT} = '-hh';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/3/86400;
  }
  elsif($timespan <= 345_600){		# four days
    $bp = '3Hf30m';		$bs = '1D';	$gi = '1h';
    $gd{PLOT_CLOCK_FORMAT} = '-hh';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'a';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/5/86400;
  }
  elsif($timespan <= 432_000){		# five days
    $bp = '4Hf1h';		$bs = '1D';	$gi = '2h';
    $gd{PLOT_CLOCK_FORMAT} = '-hh';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'a';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/5/86400;
  }
  elsif($timespan <= 691_200){		# eight days
    $bp = '6Hf1h';		$bs = '1D';	$gi = '3h';
    $gd{PLOT_CLOCK_FORMAT} = '-hh';
    $gd{PLOT_DATE_FORMAT} = '-dd o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'a';
    $gd{TIME_INTERVAL_FRACTION} = $timespan/5/86400;
    $gd{ANNOT_FONT_SIZE_SECONDARY} = '14p';
  }
  elsif($timespan <= 2_073_600){	# 24 days
    $bp = '1df3h';		$bs = '1O';	$gi = '6h';
    $gd{PLOT_DATE_FORMAT} = '-o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = '0.6';
  }
  elsif($timespan <= 3_456_000){	# 40 days
    $bp = '1df6h';		$bs = '1O';	$gi = '12h';
    $gd{PLOT_DATE_FORMAT} = '-o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = '0.6';
    $gd{ANNOT_FONT_SIZE_PRIMARY} = '12p';
  }
  elsif($timespan <= 4_320_000){	# 50 days
    $bp = '2df12h';		$bs = '1O';	$gi = '1d';
    $gd{PLOT_DATE_FORMAT} = '-o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = '0.6';
  }
  elsif($timespan <= 5_184_000){	# 60 days
    $bp = '3df12h';		$bs = '1O';	$gi = '1d';
    $gd{PLOT_DATE_FORMAT} = '-o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = '0.6';
  }
  elsif($timespan <= 6_912_000){	# 80 days
    $bp = '4df12h';		$bs = '1O';	$gi = '1d';
    $gd{PLOT_DATE_FORMAT} = '-o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = '0.6';
  }
  elsif($timespan <= 8_640_000){	# 100 days
    $bp = '7df1d';		$bs = '1O';	$gi = '1d';
    $gd{PLOT_DATE_FORMAT} = '-o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = '0.6';
  }
  elsif($timespan <= 17_280_000){	# 200 days
    $bp = '11df1d';		$bs = '1O';	$gi = '1d';
    $gd{PLOT_DATE_FORMAT} = '-o yyyy';
    $gd{TIME_FORMAT_SECONDARY} = 'a';
    $gd{TIME_INTERVAL_FRACTION} = '0.6';
  }
  elsif($timespan <= 21_600_000){	# 250 days
    $bp = '1O';		$bs = '1Y';		$gi = '1o';
    $gd{PLOT_DATE_FORMAT} = 'o';
    $gd{TIME_FORMAT_PRIMARY} = 'f';
    $gd{TIME_INTERVAL_FRACTION} = '0.8';
  }
  elsif($timespan <= 51_840_000){	# 600 days
    $bp = '1O';		$bs = '1Y';		$gi = '1o';
    $gd{PLOT_DATE_FORMAT} = 'o';
    $gd{TIME_FORMAT_PRIMARY} = 'a';
    $gd{TIME_INTERVAL_FRACTION} = '0.8';
  }
  elsif($timespan <= 94_694_400){	# three years
    $bp = '1O';		$bs = '1Y';		$gi = '1o';
    $gd{PLOT_DATE_FORMAT} = 'o';
    $gd{TIME_FORMAT_PRIMARY} = 'c';
    $gd{TIME_INTERVAL_FRACTION} = '0.8';
  }
  elsif($timespan <= 189_302_400){	# six years
    $bp = '1O';		$bs = '1Y';		$gi = '1o';
    $gd{PLOT_DATE_FORMAT} = 'o';
    $gd{TIME_FORMAT_PRIMARY} = 'c';
    $gd{TIME_INTERVAL_FRACTION} = '0.8';
    $gd{ANNOT_FONT_SIZE_PRIMARY} = '12p';
  }
  elsif($timespan <= 511_233_120){	# 16.2 years
    $bp = '1Y';		$bs = '';		$gi = '1o';
  }
  else{
    my $mult = int($timespan/511_233_120) + 1;
    $bp = "${mult}Yf1y";	$bs = '';	$gi = '1y';
  }

  $ai = $bp;
  $ai =~ s/f.*//;

  # Format the GMTDefaults arguments for psbasemap.
  my @gd;
  foreach my $k (sort keys %gd){
    my $v = $gd{$k} =~ /\s/ ? qq("$gd{$k}") : $gd{$k};
    push @gd, "--$k=$v";
  }
  my $gd = join ' ',@gd;
 
  ($bp,$bs,$ai,$gi,$gd);
}

sub read_param{
  my ($typ,$grp,$col,$st,$et) = @_;
  my @ts;

  my $span = $et - $st;

  if($span > $aggthresh1){
    # If the requested time span is greater than the threshold
    # then read the data at reduced temporal resolution.

    my $dir = '/web/aquarius/metadata/condensed_hkt';

    my $delta_t = 300;	# five minutes
    my $fnammin = sprintf "%dmin",$delta_t/60;

    if($span > $aggthresh2){
      $delta_t = 86_400;	# one day
      $fnammin = "day";
    }

    # Adjust the start and end times to the coarser time interval.
    $st -= $st % $delta_t;
    my $mod_et = $et % $delta_t;
    $et += $delta_t - $mod_et if $mod_et;

    my $numstats = 5;	# population size, mean, median, minimum, and maximum

    if($span > $aggthresh2){
      my $sy = (gmtime $st)[5] + 1900;
      my $ey = (gmtime $et)[5] + 1900;
      for(my $y=$sy; $y<=$ey; $y++){
        my $numgroups = isleap($y) ? 366 : 365;
        my $begin = timegm(0,0,0,1,0,$y-1900);
        my $f = "$dir/$y/Q$y.$fnammin.${typ}_$grp.txt";

        push @ts,
        read_condensed($f,$begin,$delta_t,$st,$et,$numgroups,$numstats,$col);
      }
    }
    else{
      my $numgroups = 86400/$delta_t;
      my $sd = int($st/86400);
      my $ed = int($et/86400);
      for(my $d=$sd; $d<=$ed; $d++){
        my ($yr,$doy) = (gmtime 86400*$d)[5,7];
        $yr  = sprintf "%04d", $yr  + 1900;
        $doy = sprintf "%03d", $doy + 1;
        my $begin = 86400*$d;
        my $f = "$dir/$yr/$doy/Q$yr$doy.$fnammin.${typ}_$grp.txt";

        push @ts,
        read_condensed($f,$begin,$delta_t,$st,$et,$numgroups,$numstats,$col);
      }
    }
  }
  else{
    # Otherwise read the data at the original temporal resolution.

    my %dir = (
    AQ => '/web/data/Aquarius/DOWNLINKS/L1',
    SP => '/web/data/Aquarius/HKT',
    );

    # Expand the time window a bit when deriving directory names because
    # stored telemetry files from CONAE can contain 3 days' worth of data,
    # and they are placed in directories according to their start times.
    my $st2 = 86400*(int($st/86400) - 3);
    for(my $t=$st2; $t<=$et; $t+=86400){
      my ($d,$m,$y) = (gmtime $t)[3..5];
      my $dir = sprintf "$dir{$typ}/%04d/%02d%02d",
                $y+1900,$m+1,$d;
      if(opendir D,$dir){
        my @f = grep /\.${typ}_$grp.txt$/,readdir D;
        closedir D;
        foreach my $f (@f){
          open F,'<',"$dir/$f" or die "Could not open $dir/$f ($!)";
          while(<F>){
            my ($ts,$val) = (split ' ')[0,$col];
            next unless defined $val and $val =~ /\d/;
            $ts = gps2utc($ts);
            next if $ts < $st;
            next if $ts > $et;
            push @ts,[$ts,$val];
          }
          close F;
        }
      }
    }
    # Deliver the time series in chronological order.
    @ts = sort { $a->[0] <=> $b->[0] } @ts;
  }
  # The strt and end times may have been modified, so I return them as well.
  ($st,$et,@ts);
}

sub minmax{
  my $ts = shift;
  my ($min,$max);
  if(@{ $ts->[0] } == 2){
    # This is a full-resolution time series.
    # Each row contains a time stamp [0] and a value [1].
    # I only care about the values here.
    ($min,$max) = map $_->[1], (sort { $a->[1] <=> $b->[1] } @$ts)[0,-1];
  }
  else{
    # This is a condensed time series.
    # Each row contains six elements:
    #    [0] time stamp
    #    [1] population size
    #    [2] mean value
    #    [3] median value
    #    [4] minimum value
    #    [5] maximum value
    # (See condense_hkt.pl around line 278.)
    # I only care about the minima and maxima here.
    $min = (sort { $a->[4] <=> $b->[4] } @$ts)[0][4];
    $max = (sort { $b->[5] <=> $a->[5] } @$ts)[0][5];
  }
  ($min,$max);
}

sub nice_increment{
  my $range = shift || 100;
  my $factor = 10**floor(log10($range));
  $factor*ceil($range/$factor)/20;
}

sub telemtyp{
  my $p = shift;
  my $t = $pindex{AQ}{$p} ? 'AQ' : ($pindex{SP}{$p} ? 'SP' : undef);
  unless($t){
    warn "Unrecognized parameter mnemonic '$p'\n";
  }
  $t;
}

# The following two subroutines (rgb2hsv and hsv2rgb) are based
# on information I found at http://en.wikipedia.org/wiki/HSL_and_HSV .

sub rgb2hsv{
  my ($r,$g,$b) = @_;
  $r /= 255;
  $g /= 255;
  $b /= 255;
  my $min = $g<=$b ? ($r<=$g ? $r : $g) : ($r<=$b ? $r : $b);
  my $max = $g>=$b ? ($r>=$g ? $r : $g) : ($r>=$b ? $r : $b);
  my $value = $max;
  my $chroma = $max - $min;
  my $hue;
  if($chroma == 0){
    $hue = 0;
  }
  elsif($max == $r){
    $hue = fmod( ($g - $b)/$chroma , 6 );
  }
  elsif($max == $g){
    $hue = ($b - $r)/$chroma + 2;
  }
  else{
    $hue = ($r - $g)/$chroma + 4;
  }
  $hue *= 60;
  my $sat = $chroma == 0 || $value == 0 ? 0 : $chroma/$value;
  ($hue,$sat,$value);
}

sub hsv2rgb{
  my ($h,$s,$v) = @_;
  my $chroma = $v*$s;
  $h /= 60;
  my $x = $chroma*(1 - abs(fmod($h,2) - 1));
  my ($r,$g,$b) = $h < 1 ? ($chroma,$x,0) :
                ( $h < 2 ? ($x,$chroma,0) :
                ( $h < 3 ? (0,$chroma,$x) :
                ( $h < 4 ? (0,$x,$chroma) :
                ( $h < 5 ? ($x,0,$chroma) :
                           ($chroma,0,$x)  ))));
  my $m = $v - $chroma;

  map int(255*($_ + $m) + 0.5), $r,$g,$b;
}

sub doit{
  my $cmd = join ' ', split ' ', shift;
  #print STDERR "[",scalar localtime,"] $cmd\n";
  system($cmd) == 0 or die "Error executing $cmd ($!)";
}
