#!/usr/bin/env perl

# Norman Kuring         26-Jun-2018     Initial development
# Norman Kuring          1-Oct-2018

use strict;
use Time::Local qw(timegm);

my ($gmt_J,$gmt_R) = @ARGV;
$gmt_J =~ /^-J/ and $gmt_R =~ /^-R/ or die <<"--";
Usage:
$0 GMT_J_option GMT_R_option
See https://gmt.soest.hawaii.edu/doc/5.4.2/mapproject.html#required-arguments .
Note that this scriptuses a restricted subset of the full GMT -R and -J
specifications.
--

# Read scene identifiers -- one per line -- from standard input.
my @scenes;
while(<STDIN>){
  if(/^(V\d{13})/){
    push @scenes, $1;
  }
  else{
    warn "Unrecognized scene identifier: $_";
  }
}

# Sort the scenes and remove inadvertent duplicates.
my %scenes;
grep $scenes{$_}++ && 0, @scenes;
@scenes = sort keys %scenes;

# Group the scenes by orbit by assuming that a gap of more than
# 40 minutes (2400 seconds) indicates a new dayside portion of the orbit.
my @orbits;
my $i = -1;
my $prev_time = 0;
foreach my $scene (@scenes){
  my ($y,$d,$h,$m,$s) = $scene =~
     /^V(\d{4})(\d{3})(\d\d)(\d\d)(\d\d)/;
  my $t = timegm($s,$m,$h,1,0,$y-1900) + 86400*($d-1);
  if($t - $prev_time > 2400){
    $i++;
  }
  $prev_time = $t;
  push @{$orbits[$i]},$scene;
}

# Get the dimensions of the resulting images.
# Stipulate a display resolution of 100 pixels per inch.
doit("gmt set PROJ_LENGTH_UNIT inch");
my ($wid_in,$hgt_in) = split ' ',`gmt mapproject -W $gmt_J $gmt_R`;
my ($wid_px,$hgt_px) = map int(100*$_ + 0.5), $wid_in,$hgt_in;
   ($wid_in,$hgt_in) = map     $_/100       , $wid_px,$hgt_px;


# Download and process each scene to level-2.  Make red/green/blue rhos
# products plus rhos_862 for help with land/water discrimination.
# Also generate chlor_a product.
foreach my $scene (@scenes){
  doit("
  wget -q --load-cookies ~/.urs_cookies --save-cookies ~/.urs_cookies --auth-no-challenge=on --content-disposition https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/$scene.L1A_SNPP.nc;
  wget -q --load-cookies ~/.urs_cookies --save-cookies ~/.urs_cookies --auth-no-challenge=on --content-disposition https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/$scene.GEO-M_SNPP.nc;
  getanc.py $scene.L1A_SNPP.nc;
  l2gen ifile=$scene.L1A_SNPP.nc
        geofile=$scene.GEO-M_SNPP.nc
        par=$scene.L1A_SNPP.nc.anc
        ofile=$scene.L2.nc
        l2prod=rhos_443,rhos_551,rhos_671,rhos_862,chlor_a;
  rm $scene.L1A_SNPP.nc $scene.GEO-M_SNPP.nc
  ");
}

# For each orbit, extract the longitude, latitude, and reflectance for
# each band, then project the data and grid them to the output raster form.
my $tmppfx = (split/\//,$0)[-1] . ".$^T.$$";
my $nodata = pack "f",-32767;
foreach my $orb (@orbits){
  my $pfx = $orb->[0];
  foreach my $scene (@$orb){
    foreach my $p (qw(longitude latitude)){
      doit("
      h5dump -d /navigation_data/$p -b -o $tmppfx $scene.L2.nc;
      cat $tmppfx >>$pfx.$p
      ");
    }
    foreach my $p (qw(rhos_443 rhos_551 rhos_671 rhos_862)){
      doit("
      h5dump -d /geophysical_data/$p -b -o $tmppfx $scene.L2.nc;
      cat $tmppfx >>$pfx.$p
      ");
    }
  }
  foreach my $p (qw(rhos_443 rhos_551 rhos_671 rhos_862)){
    open X,'<',"$pfx.longitude" or die "Could not open $pfx.longitude ($!)";
    open Y,'<',"$pfx.latitude"  or die "Could not open $pfx.latitude ($!)";
    open Z,'<',"$pfx.$p"        or die "Could not open $pfx.$p ($!)";
    my $cmd = join ' ', split ' ', <<"    --";
    gmt mapproject $gmt_J $gmt_R -bi3f -bo3f |
    gmt nearneighbor -G$pfx.$p.grd -I0.01 -R0/$wid_in/0/$hgt_in
                     -bif -N4 -S0.06 -r
    --
    open P,'|-',$cmd or die "Error opening pipe to $cmd ($!)";
    my ($x,$y,$z);
    while(read Z,$z,4){
      read X,$x,4;
      read Y,$y,4;
      next if $z eq $nodata;
      print P $x,$y,$z;
    }
    close P or die "Error executing command:\n$cmd ($!)";
    close X; close Y; close Z;
    doit(qq(
    gmt grdmath $pfx.$p.grd
    0.005 MAX 0.005 DIV LOG 2 0.005 DIV LOG DIV 1 MIN 65535 MUL = $tmppfx.grd;
    gmt grd2xyz $tmppfx.grd -do0 -ZTLHw >$tmppfx.dat;
    /bin/echo -e "P5\\n$wid_px $hgt_px\\n65535"|cat - $tmppfx.dat > $pfx.$p.pgm
    ));
  }
  doit(qq(
  convert $pfx.rhos_671.pgm $pfx.rhos_551.pgm $pfx.rhos_443.pgm
  -combine $pfx.rhos_671_551_443.tif
  ));
}

# Clean up a bit.
doit("/bin/rm $tmppfx*");

# Use these minima and maxima to log-scale the chlor_a and PIC data.
my %minmax = (
chlor_a => [0.01   ,50   ]
);

# For chlor_a, just throw all of the data into the hopper at once
# without bothering to sparate according to orbit.
my $pfx = substr $scenes[0],0,8;
foreach my $scene (@scenes){
  foreach my $p (qw(chlor_a)){
    doit("
    h5dump -d /geophysical_data/$p -b -o $tmppfx $scene.L2.nc;
    cat $tmppfx >>$pfx.$p
    ");
  }
}
# Concatenate the previously extracted longitude and latitude data.
doit("cat *.longitude >$pfx.lon; cat *.latitude  >$pfx.lat");
# Project and grid the products.
foreach my $p (qw(chlor_a)){
  open X,'<',"$pfx.lon" or die "Could not open $pfx.lon ($!)";
  open Y,'<',"$pfx.lat" or die "Could not open $pfx.lat ($!)";
  open Z,'<',"$pfx.$p"  or die "Could not open $pfx.$p ($!)";
  my $cmd = join ' ', split ' ', <<"  --";
  gmt mapproject $gmt_J $gmt_R -bi3f -bo3f |
  gmt nearneighbor -G$pfx.$p.grd -I0.01 -R0/$wid_in/0/$hgt_in
                   -bif -N4 -S0.06 -r
  --
  open P,'|-',$cmd or die "Error opening pipe to $cmd ($!)";
  my ($x,$y,$z);
  while(read Z,$z,4){
    read X,$x,4;
    read Y,$y,4;
    next if $z eq $nodata;
    print P $x,$y,$z;
  }
  close P or die "Error executing command:\n$cmd ($!)";
  close X; close Y; close Z;
  doit(qq(
  gmt grdmath $pfx.$p.grd
  $minmax{$p}[0] MAX $minmax{$p}[0] DIV LOG
  $minmax{$p}[1] $minmax{$p}[0] DIV LOG DIV
  1 MIN 65535 MUL = $tmppfx.grd;
  gmt grd2xyz $tmppfx.grd -do0 -ZTLHw >$tmppfx.dat;
  /bin/echo -e "P5\\n$wid_px $hgt_px\\n65535"|cat - $tmppfx.dat
  | convert - $pfx.$p.tif
  ));
}




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