#!/usr/bin/perl 

# $Id: topo2grd 431 2011-09-22 16:29:56Z cbworden $ USGS

# A script to convert DEM grids into a topo correction map 
# for ShakeMap using the slope-to-vs30 calculation from Wald
# and Allen.

use strict;
use FindBin;
use lib "$FindBin::Bin/../lib";
use ShakeConfig;
use File::Copy;
use POSIX;

#---------------------------------------------
# Global variables

my $topo_dir = $DEMdir;
my $VS30_OCEAN = 600;
my $VS30_MINIMUM = 180;
my $coast_res = 'f';
my ($XMIN,$XMAX,$YMIN,$YMAX);
my $VERBOSE = 1;
my $v = $VERBOSE ? '-V' : '';

#---------------------------------------------

print (join ' ',($0,@ARGV),"\n") if $VERBOSE;

my ($evid,$bounds,$dxdy,$regime, $demonly) = @ARGV;
my ($dx,$dy) = split "/", $dxdy;

die "Usage: $0 evid wesnbounds dx/dy [regime=active|stable]\n" 
  unless ($evid and $bounds and $dxdy);

my $tmp = "tmp.topo2grd.$$";

my $outputdir = ($demonly) ? "$shake_home/data/$evid/mapping" : "$shake_home/data/$evid/output";
my $grd_topo = "$outputdir/topo_data.grd";
my $qtm_file = "$outputdir/sitecorr.grd";

die "$outputdir does not exist" unless (-d $outputdir);

#my $res = '30c';
my $res = convert_res($dx) <= convert_res($dy) ? $dx : $dy;
print "Using $topo_dir. Final resolution fixed at $res.\n" if $VERBOSE;

my $command;
my $results = dem2grd($bounds, $tmp, $res);
print $results if $VERBOSE;

my ($mask_bounds) = ($results =~ /FINAL BOUNDS: (.*?)\n/);
die "Could not get final bounds (need newer version of dem2grd.pl?)" 
  unless ($mask_bounds);


unless (-s "$topo_dir/$tmp.grd") {
  print "Could not find $topo_dir/$tmp.grd, aborting.\n";
  exit;
}

print `mv $topo_dir/$tmp.grd $grd_topo`;
exit if ($demonly);
print `$gmt_bin/grdgradient $grd_topo -G$grd_topo.i -S$grd_topo.grad -M -D $v`;

#-----------------------------------------------------------------
# Add ocean mask here and set ocean gridpoints to -9999
#-----------------------------------------------------------------

my $map_bounds = $mask_bounds;
$map_bounds =~ s/^-R//;
my ($west, $east, $south, $north) = split /\//,$map_bounds;
my $pad = convert_res($res)*25;
$map_bounds = sprintf("-R%f/%f/%f/%f", $west-$pad, $east+$pad, $south-$pad, $north+$pad);
my $landmask = "$grd_topo.landmask";
$command = "$gmt_bin/grdlandmask $map_bounds -I$res "
  ." -G$landmask -D$coast_res -N-9999/NaN -F";
print $command."\n" if $VERBOSE;
system($command);
die "Could not create file $landmask" unless (-s $landmask);
print("$gmt_bin/grdsample $landmask -G$landmask  $mask_bounds -I$res\n") if $VERBOSE;
system("$gmt_bin/grdsample $landmask -G$landmask  $mask_bounds -I$res");

$command = "$gmt_bin/grdmath $landmask $grd_topo.grad AND = $grd_topo.grad.1";
print $command."\n" if $VERBOSE;
system($command);
die "Could not use landmask on $grd_topo.grad" unless (-s "$grd_topo.grad.1");
copy("$grd_topo.grad.1","$grd_topo.grad");

#-----------------------------------------------------------------
# Figure out which coefficients to use. In normal (and now
# deprecated) usage, use get_slope function to get average
# slope over the map area, then compare it to an arbitrary
# cutoff point to select active or stable coefficients.
# 
# Now we support the regime flag which takes 'active' and 
# 'stable'
#-----------------------------------------------------------------

if (defined $regime) {
  ($regime) = ($regime =~ /regime=(.*)$/);
  chomp $regime;

  die "Could not parse regime flag (need regime=active|stable)"
    unless ($regime);

}
else {
  my $slope = get_slope("$grd_topo.grad");
  print "Mean slope for this map is $slope.\n" if $VERBOSE;
  $regime = ($slope > 0.05) ? 'active' : 'stable';
}

my $slopefunc = slope2vs30($regime);

#-----------------------------------------------------------------
# End coefficents block, start GMT processing
#-----------------------------------------------------------------

system ("$gmt_bin/grd2xyz $grd_topo.grad > $grd_topo.xyz");
convert("$grd_topo.xyz","$grd_topo.vs30");

# Now .vs30 is an xyz file with the vs30 values.
# Convert to the right resolution:

system ("$gmt_bin/nearneighbor $grd_topo.vs30 -I$dxdy -R$bounds -N1/1 -S$res -G$outputdir/sitecorr.grd.$res");
system ("$gmt_bin/grdsample $outputdir/sitecorr.grd.$res -G$outputdir/sitecorr.grd.gmt -I$dxdy -Qn ");

system("/bin/rm $outputdir/sitecorr.grd.$res $grd_topo*");

#-----------------------------------------------------------------
# Begin subroutines
#-----------------------------------------------------------------

sub dem2grd {
  my ($coords, $name, $res) = @_;

  my @coords = split /\//,$coords;

  my $raw_res = "30c";
  $#coords == 3 or die "Usage: $0 lon0/lon1/lat0/lat1 event [res]\n";

  my (@lon,@lat);
  ($lon[0],$lon[1],$lat[0],$lat[1]) = @coords;
  $lon[1] = 180 if ($lon[1] > 180);
  $lon[0] = -180 if ($lon[0] < -180);

  my @lonbounds = ( "w180 -180 -140",
		  "w140 -140 -100",
		  "w100 -100 -60",
		  "w060 -60 -20",
		  "w020 -20 20",
		  "e020 20 60",
		  "e060 60 100",
		  "e100 100 140",
		  "e140 140 180",
		);

  my @latbounds = ( "n90 40 90",
		  "n40 -10 40",
		  "s10 -60 -10",
		);

  # Adjust lat/lon to match raw grid resolution

  my $raw_res1 = convert_res($raw_res);

  # Iterate through available files

  my ($try,$bound1,$bound2);
  my @list;

  my ($w,$e) = @lon;
  my ($s,$n) = @lat;

  my $fullgrd = "$topo_dir/$name.grd";
  unlink $fullgrd if (-e $fullgrd);

  my $n_extract = 0;

  foreach my $lonbound (@lonbounds) {
    my ($lon,$wbound,$ebound) = split /\s+/,$lonbound;

    my $rowgrd = "$topo_dir/tmp.$$.$lon";
    unlink $rowgrd if (-e $rowgrd);

    foreach my $latbound (@latbounds) {
      my ($lat,$sbound,$nbound) = split /\s+/,$latbound;
 
      my $file = "$topo_dir/DATA/$lon$lat.grd"; 
      if (not -e $file) {
        print "DEM grid file: $file NOT FOUND. Why?\n" if $VERBOSE;
        next;
      }

      if ($sbound > $n) {
        print "Skipping $file: $sbound > $n\n" if $VERBOSE;
        next;
      }
      if ($nbound < $s) {
        print "Skipping $file: $nbound < $s\n" if $VERBOSE;
        next;
      }
      if ($wbound > $e) {
        print "Skipping $file: $wbound > $e\n" if $VERBOSE;
        next;
      }
      if ($ebound < $w) {
        print "Skipping $file: $ebound < $w\n" if $VERBOSE;
        next;
      }
      print "Proceeding to extract $file\n" if $VERBOSE;
      $n_extract++;

      my $tmpgrd = "$topo_dir/tmp.$$.$lon$lat";
      extract($w,$e,$s,$n,$wbound,$ebound,$sbound,$nbound,
        $file,$tmpgrd,$raw_res1);
      paste($rowgrd,$tmpgrd);
      unlink $tmpgrd;
    }
    # At this point, $rowgrd has the grdfile for this row
    if (-e $rowgrd) {
      paste($fullgrd,$rowgrd);
      unlink $rowgrd;
    }
  }

  if ($n_extract == 0) {
    print "Error in topo2grd: Nothing extracted.\n"
        . "Set \$VERBOSE = 1 and examine the output.\n"
        . "It is likely you have something misconfigured in your DEM data.\n";
    die "Quitting topo2grd.";
  }

  print "Done with pasting.\n" if ($VERBOSE);
  my $maxbounds = "-R$XMIN/$XMAX/$YMIN/$YMAX";

  print "FINAL BOUNDS: $maxbounds\n" if $VERBOSE;

  my $res1 = convert_res($res);
  if (!$res1 or $res1 == $raw_res1) {
    return "FINAL BOUNDS: $maxbounds\n";
  }

  my $rflag = "-R$lon[0]/$lon[1]/$lat[0]/$lat[1]";
  my $tmp = "$fullgrd.1";
  command("cp $fullgrd $tmp");
  command("$gmt_bin/grdsample $tmp -G$fullgrd -I$res");
  unlink $tmp;

  return "FINAL BOUNDS: $maxbounds\n";

}

sub get_slope {
  my $file = shift;
  my $result;
  
  my $ghd = new GMThd($file, $gmt_bin);
  ($result, undef, undef) = $ghd->mean();

  return $result;
}

sub slope2vs30 {
  my $regime = shift;
  if ($regime eq 'active') {
    print "Using active tectonic factors.\n" if $VERBOSE;
    return sub {
      my $s = shift;

      return $VS30_OCEAN if ($s == -9999 or $s == 0);
      my $vs30 = ($s < 0.000032) ? $VS30_MINIMUM :
              ($s < 0.0022)   ? 210 :
              ($s < 0.0063)   ? 270 :
              ($s < 0.018)    ? 330 :
              ($s < 0.05)     ? 425 :
              ($s < 0.1)      ? 555 :
              ($s < 0.138)    ? 690 :
                                800;
      return $vs30;
    };
  }
  elsif ($regime eq 'stable') {
    print "Using stable continent factors.\n" if $VERBOSE;
    return sub {
      my $s = shift;

      return $VS30_OCEAN if ($s == -9999 or $s == 0);
      my $vs30 = ($s < 0.000001) ? $VS30_MINIMUM :
              ($s < 0.002)   ? 210 :
              ($s < 0.004)   ? 270 :
              ($s < 0.0072)   ? 330 :
              ($s < 0.013)    ? 425 :
              ($s < 0.018)    ? 555 :
              ($s < 0.025)    ? 690 :
                                760;
      return $vs30;
    };
  }
  die "ERROR: Unknown tectonic regime $regime";
}
 
sub convert {
  my ($in,$out) = @_;
  open IN,"<$in";
  open OUT,">$out";


  while (<IN>) {
    my ($x,$y,$slope) = split /\s+/;
    my $vs30 = $slopefunc->($slope);
    print OUT $x,' ',$y,' ',$vs30,"\n";
  }  
  close IN;
  close OUT;
}

sub extract {
  my ($w,$e,$s,$n,$wbound,$ebound,$sbound,$nbound,$grdfile,$name,
      $raw_res1) = @_;

  my $w1 = ($w > $wbound) ? $w : $wbound;
  my $e1 = ($e < $ebound) ? $e : $ebound;
  my $n1 = ($n < $nbound) ? $n : $nbound;
  my $s1 = ($s > $sbound) ? $s : $sbound;

  $w1 = offset($w1,$raw_res1,-1);
  $e1 = offset($e1,$raw_res1,1);
  $s1 = offset($s1,$raw_res1,-1);
  $n1 = offset($n1,$raw_res1,1);
 
  my $w1 = ($w1 > $wbound) ? $w1 : $wbound;
  my $e1 = ($e1 < $ebound) ? $e1 : $ebound;
  my $n1 = ($n1 < $nbound) ? $n1 : $nbound;
  my $s1 = ($s1 > $sbound) ? $s1 : $sbound;

  print "EXTRACTING $grdfile W:$w1 E:$e1 S:$s1 N:$n1\n" if $VERBOSE;

  my $command = "$gmt_bin/grdcut -R$w1/$e1/$s1/$n1 $grdfile -G$name";
  print "$command\n" if $VERBOSE;
  command($command);
  die "Could not extract $grdfile\n" unless (-e $name);

  $XMIN = $w1 unless (defined $XMIN and $XMIN < $w1);
  $XMAX = $e1 unless (defined $XMAX and $XMAX > $e1);
  $YMIN = $s1 unless (defined $YMIN and $YMIN < $s1);
  $YMAX = $n1 unless (defined $YMAX and $YMAX > $n1);

}

sub paste {
  my ($f1,$f2) = @_;

  if (!-e $f1) {
    command("mv $f2 $f1");
    return;
  }
  command("$gmt_bin/grdpaste $f1 $f2 -G$f1.tmp");
  command("mv $f1.tmp $f1");  
} 

sub convert_res {
  my $res = shift;
  my $res1;

  return $res unless ($res =~ /m|c/);
  if ($res =~ /^([0-9]*\.?[0-9]+)([m|c])$/) {
    $res1 = $1/60 if ($2 eq 'm');
    $res1 = $1/3600 if ($2 eq 'c');
  }
  elsif ($res =~ /^[\d\.]$/) {
    return $res;
  }
  else {
    die "Invalid res $res";
  }

  die "Could not convert $res" unless ($res1 > 0);
  return $res1;
}

sub command {
  my $command = shift;
  print $command."\n" if ($VERBOSE);
  return `$command`;
}

sub offset {
  my ($deg,$res,$dir) = @_;

  my $val = $deg/$res;
  return $deg if ($val == int $val);

  $val = int $val + $dir;
  $val *= $res;
  
  return $val;
}