#! /bin/csh

if ($#argv != 13) then
    echo "usage: `basename $0` img imgIncInSec arcticGrd ofile W E S N incInSec minArcticLat maxPredLat sidOfArctic oDir"
    echo "  example: `basename $0` ../img/topo_9.1b.img 60 ../predTile/w140n40.grd -141 -99 -11 41 30 79 80 65535 /tmp"
    echo "  don't forget to have a 1 deg pad"
    echo
    echo " increment must be numerical seconds, i.e. 75 -NOT- 1.25m"
    echo
    exit
endif

# need lots of sig figs to keep surface from moaning

gmtset D_FORMAT = %.16lg

# basic idea is to blockmedian and surface the predicted and then 
# cut out bad stuff about 70 N and paste in arctic bathymetry.


set img         = $1;   shift
set imgIncInSec = $1;   shift
set arcticGrd   = $1;   shift
set ofile       = $1;   shift
set w           = $1;   shift
set e           = $1;   shift
set s           = $1;   shift
set n           = $1;   shift
set incInSec    = $1;   shift
set minArctic   = $1;   shift
set maxPred     = $1;   shift
set sidOfArctic = $1;   shift
set oDir        = $1;   shift

set sidOfPred   = 0

set stem        = `basename -s .grd $ofile`
set dir         = `dirname  $ofile`

# fiddle with increment; we need both number and string

set inc = `echo $incInSec | awk '{ printf "-I%sc\n", $1 }'`

# much screwing around with WESN stuff...
#
# add 1 degree pad assumed to be in xyzi data
# slightly tricky, and have to mess around with 180 rollover
# csh only does integer math so need integer values for lat and lon.

set padWidth = 1
set w = `echo $w|awk '{ tmp = $1-'$padWidth'; printf "%.16lg\n", tmp }'`
set e = `echo $e|awk '{ tmp = $1+'$padWidth'; printf "%.16lg\n", tmp }'`
set s = `echo $s|awk '{ tmp = $1-'$padWidth'; if (tmp < -90) tmp = -90; printf "%.16lg\n", tmp }'`
set n = `echo $n|awk '{ tmp = $1+'$padWidth'; if (tmp > +90) tmp = +90; printf "%.16lg\n", tmp }'`

set w360 = `echo $w|awk '{ tmp = $1; if (tmp < 0) tmp = tmp+360; printf "%.16lg\n", tmp }'`
set e360 = `echo $e|awk '{ tmp = $1; if (tmp < 0) tmp = tmp+360; printf "%.16lg\n", tmp }'`

set region = "-R/$w/$e/$s/$n"


# and a few tmp files...

set tmpFile0 = $oDir/tmp/`basename $0`.$stem.$incInSec.mercator
set tmpFile1 = $oDir/tmp/`basename $0`.$stem.$incInSec.mercator.plusMinus180

# whew, back to the problem at hand


# Surface img. 
#
# If we are in arctic, we have to limit predicted to it's region of good data, with arctic
# grid above that.
#
# Odd thing here is that img2web returns mercator PIXEL projection. 
# Remove that with a surface onto a geographic NODE grid.

# First get predicted bathymetry out of PIXEL registered 0-360 binary file, then convert
# into +/-180 that GMT understands better, then blockemedian and surface to remove mercator.

# note: $maxPred is inefficient in call to img2web. we only need up to $n
set img2webNorth = `echo $n $maxPred  | awk '{ if ($1 < $2) print $1; else print $2; }'`

if ($w360 < $e360) then # NOT 360 rollover, but might cross +/- 180
    img2web $img $s $img2webNorth $w360 $e360 1 >!  $tmpFile0.xyz
else # 360 roll over
    img2web $img $s $img2webNorth $w360 360   1 >!  $tmpFile0.xyz
    img2web $img $s $img2webNorth 0     $e360 1 >>! $tmpFile0.xyz
endif

# deal with +/-180 rollover

if ( $w > 0 ) then
    # nothing needs to be done, although e140 might seem like a special case.
else # west edge is in west hemisphere
    if ( $e > 0) then
        cat $tmpFile0.xyz|\
            awk '{ tmp = $1; if (tmp > 180) tmp = tmp-360; printf "%.16lg %.16lg %.16lg\n", tmp, $2, $3 }' >!  $tmpFile1.xyz
        mv $tmpFile1.xyz $tmpFile0.xyz
    else
        cat $tmpFile0.xyz|\
            awk '{ printf "%.16lg %.16lg %.16lg\n", $1-360, $2, $3 }' >!  $tmpFile1.xyz
        mv $tmpFile1.xyz $tmpFile0.xyz
    endif
endif


# At this point we have a PIXEL registered MERCATOR grid of predicted with -180<lon<180. 
#
# If we are north of +70 get arctic bathymetry and block median with predicted along 
# 1 deg wide transition.
#
# FIXME: arctic grd does NOT include -180, just plus 180. This leaves a line along 180

if ($n >= $minArctic) then
    grd2xyz $arcticGrd -R/$w/$e/$minArctic/$n -S >>! $tmpFile0.xyz
endif

# At this point we have PIXEL mercator projected "predicted" & NODE registered arctic xyz.

# Mixing these two would be a problem, but we know that the predicted has no wavelengths
# shorter than 1 min, so shifting the grid by 30c or less is harmless. This isn't true in
# general of course.

# Block median into final NODE registered -and- QUICK grid or surface complains about 
# multiple hits in single node. Need to surface to fill in holes from mercator projection
# and so we can resample everything to PIXEL. grdsample does NOT handle NaN well...

cat $tmpFile0.xyz|\
    blockmedian -Q $inc $region|\
        surface    $inc $region -G$tmpFile0.grd 

# Convert to PIXEL, and we are done with depths!
# grdsample $tmpFile0.grd -G$dir/$stem.grd -Ql -F $region
mv $tmpFile0.grd $dir/$stem.grd

# Make sid. It's easy, just use grdmath to make a constant matrix by multiplying depth 
# grid by zero and adding sid.
#
# Note: we can -NOT- filter sid in anyway!

grdmath $dir/$stem.grd 0 MUL $sidOfPred ADD = $dir/$stem.sid.grd


# -IF- we used arctic grid, then we need to account for that...

if ($n > $maxPred) then
    
    # Cut off arctic data, set it to arctic sid, then paste it back with the predicted.
    
    set tmpFile2 = $oDir/tmp/`basename $0`.$stem.$incInSec.sid.a
    set tmpFile3 = $oDir/tmp/`basename $0`.$stem.$incInSec.sid.b
    set tmpFile4 = $oDir/tmp/`basename $0`.$stem.$incInSec.sid.c

	if ($s < $maxPred) then
    	grdcut   $dir/$stem.sid.grd -R/$w/$e/$maxPred/$n -G$tmpFile2.grd     
    	grdmath  $tmpFile2.grd ISNAN 0 MUL $sidOfArctic ADD = $tmpFile4.grd
    	grdcut   $dir/$stem.sid.grd -R/$w/$e/$s/$maxPred -G$tmpFile3.grd
    	grdpaste $tmpFile3.grd $tmpFile4.grd -G$dir/$stem.sid.grd
	else 
    grdcut   $dir/$stem.sid.grd -R/$w/$e/$s/$n -G$tmpFile2.grd     
    grdmath  $tmpFile2.grd ISNAN 0 MUL $sidOfArctic ADD = $tmpFile4.grd
    mv $tmpFile4.grd $dir/$stem.sid.grd	
    endif
endif

# all done

/bin/rm -rf $oDir/tmp/`basename $0`.$stem* 

