#!/usr/bin/perl -w

# Norman Kuring		 8-Feb-2012	Initial development
# Norman Kuring		15-May-2018	Allow user to specify projection.
# Norman Kuring		10-Oct-2018	Pass GMT -J and -R arguments verbatim.
# Norman Kuring		20-Dec-2018	Add -S option to mapproject.

use strict;

my $hdp = '/usr/bin/hdp';
my $gmtset = '/usr/bin/gmt set';
my $mapproject = '/usr/bin/gmt mapproject -S';
my $nearneighbor = '/usr/bin/gmt nearneighbor';
my $grdimage = '/usr/bin/gmt grdimage';
my $gs = '/usr/bin/gs';
my $convert = '/usr/bin/convert';

my ($jarg,$rarg) = @ARGV;
defined $rarg or die <<"--";
Usage:
$0 GMT_J_option GMT_R_option
--

### Set up a work space ####################################################
my $cleanup = &workspace;

my ($cmd,$d);

doit("$gmtset PROJ_LENGTH_UNIT inch");

my $prjcmd = "$mapproject $jarg $rarg";

my ($wid,$hgt) = split ' ',`$prjcmd -W`;
my ($w100,$h100) = map int(100*$_ + 0.5), $wid,$hgt;
($wid,$hgt) = map $_/100, $w100,$h100;

my $grdcmd = "$prjcmd|$nearneighbor -I0.01 -R0/$wid/0/$hgt -N4 -S0.04 -r";

$cmd = "$grdcmd -Gred.grd";
open RP,'|-',$cmd or die "Error opening pipe to $cmd ($!)";
$cmd = "$grdcmd -Ggrn.grd";
open GP,'|-',$cmd or die "Error opening pipe to $cmd ($!)";
$cmd = "$grdcmd -Gblu.grd";
open BP,'|-',$cmd or die "Error opening pipe to $cmd ($!)";

while(my $path = <STDIN>){
  chop $path;
  my $f = (split /\//,$path)[-1];

  print STDERR "[",scalar localtime,"] Reading $path\n";

  $cmd = "$hdp dumpsds -n longitude -b -d $path > lon";
  system($cmd) == 0 or die "Error executing $cmd ($!)";
  $cmd = "$hdp dumpsds -n latitude -b -d $path > lat";
  system($cmd) == 0 or die "Error executing $cmd ($!)";
  $cmd = "$hdp dumpgr -pd -b $path > pal";
  system($cmd) == 0 or die "Error executing $cmd ($!)";
  $cmd = "$hdp dumpgr -i 0 -d -b $path > img";
  system($cmd) == 0 or die "Error executing $cmd ($!)";

  open F,'<',"pal" or die "Could not open pal ($!)";
  my @lut;
  for(my $i=0;$i<256;$i++){
    read(F,$d,3) == 3 or die "Error reading palette ($!)";
    @{ $lut[$i] } = unpack "CCC",$d;
  }
  close F;
  unlink "pal";

  open I,'<',"img" or die "Could not open img ($!)";
  open L,'<',"lon" or die "Could not open lon ($!)";
  open P,'<',"lat" or die "Could not open lat ($!)";
  while(read I,$d,1){
    my ($lon,$lat);
    read L,$lon,2;
    read P,$lat,2;
    $lon = unpack "s",$lon;
    $lat = unpack "s",$lat;
    $d   = unpack "C",$d;
    $lon /= 180;
    $lat /= 360;
    print RP "$lon $lat $lut[$d][0]\n";
    print GP "$lon $lat $lut[$d][1]\n";
    print BP "$lon $lat $lut[$d][2]\n";
  }
  close I;
}
close RP;
close GP;
close BP;

doit("
$grdimage red.grd grn.grd blu.grd -JX$wid/$hgt -P -X0 -Y0 -Q
| $gs -sDEVICE=pngalpha -sOutputFile=- -q -r100 -g${w100}x$h100 -
| $convert - -transparent '#808080' rgb.png
");

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

### Clean up the work space ################################################
&$cleanup;



sub vet_number{
  my $str = shift;
  $str =~ /^\s*([+\-]?(?:\d+\.?\d*|\.\d+))\s*$/;
  $1;
}

# Create a workspace and go there then return a subroutine
# that will get rid of the workspace.
sub workspace{
  my $tmpdir = join '','/dev/shm/',(split /\//,$0)[-1],'.',time,'.',$$;
  mkdir $tmpdir,0755 or die "Could not create directory, $tmpdir ($!)";
  chdir $tmpdir or die "Could not cd to $tmpdir ($!)";
  return sub{
    chdir;
    system("/bin/rm -r $tmpdir") == 0 or die "Could not remove $tmpdir ($!)";
  };
}

sub doit{
  my $cmd = join ' ', split ' ', shift;
  system($cmd) == 0 or die "Error executing $cmd ($!)";
}
