#! /bin/csh

if ($#argv != 7) then
	echo "usage: `basename $0` xyziFile W E S N incInSec oDir"
	echo "  example: `basename $0` ../cmTile/w140n40.xyzi -141 -99 -11 41 3 .."
	echo "  don't forget to have a 1 deg pad"
	echo
	exit
endif

# ignoring all details here is a "flow chart"
#
# note "polishing", i.e. two paths leading to one really cool file.
#
#	use img2web,
# 	then separately surface [z-pred] and pred, 
#	then add surfaced outputs.
#
#
#   "pixels from cm"        "predicted" (pixel grid)        
#       |                                   |
#   xyzi (random ascii)                     |
#       |                                   |
#   blockmedian (pixel ascii)               |
#   grdtrack (pixel ascii) < ---------------|
#   z-predicted (pixel ascii)               |
#   surface (pixel grid)                    |
#       |                                   |
#       --->grdmath "add" (pixel grid) <---
#           grdcut  (pixel grid)
#             |
#             |
#	"polished cmdata in a pixel GTOP30 tile"

# need lots of sig figs to keep surface from moaning

gmtset D_FORMAT = %.16lg


set ifile	= $1;   shift
set w		= $1;   shift
set e		= $1;   shift
set s		= $1;   shift
set n		= $1;   shift
set incInSec = $1;  shift
set oDir 	= $1;   shift
set opts	= "$*"

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

set stem = `basename -s .xyzi $ifile`
set predicted = $oDir/predTile/$stem.grd
set ping      = $oDir/tmp/`basename $0`.$stem.$incInSec.ping
set landmask  = $oDir/tmp/`basename $0`.$stem.$incInSec.landmask
set polished  = $oDir/tmp/`basename $0`.$stem.$incInSec.polished
set sidTmp    = $oDir/tmp/`basename $0`.$stem.$incInSec.sidTmp

# deal with csh nonsense of not groking plus sign (+21 must be 21)

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

# Tension in surface

set tension = -T0.65
set tension = -T0.25
set tension = -T1.0

# Prepare to subtract 1 degree pad assumed to be in xyzi data

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

# define working (padded) and final (trimmed) regions

set region = "-R$w/$e/$s/$n"
set finalRegion = "-R/$finalW/$finalE/$finalS/$finalN"

# set lt0 = `echo $s`
# set ltf = `echo $n`
# @ ltm = ( $lt0 + $ltf ) / 2 

# whew, back to the job at hand

echo "working on the $stem tile"


# blockmedian ping data
#
# NOTE use hacked version of blockmedian that puts sid of median
# ping in column 4
#
# sid must be NODE -AND- we can -NOT- filter sid in anyway
#
# use -Q option so xy returned is xy of median z, and avoid hard to debug beating!!!
#
# place 0 everywhere SRTM has land to keep interpolation by surface reasonable near land

# srtmTileMaker/makeArbitrarySRTM.sh $w $e $s $n `echo $incInSec | awk '{ printf "%sc\n", $1 }'` 
srtmTileMaker/makeArbitrarySRTM.sh $w $e $s $n `echo $incInSec | awk '{ printf "%sc\n", $1 }'` $oDir/srtm/grd/$stem.grd
exit
grdmath   $oDir/srtm/grd/$stem.grd 0 MUL = $landmask.grd
grdcut $finalRegion $oDir/srtm/grd/$stem.grd \
	-G$oDir/srtm/grd/$stem.grdcut.grd
mv $oDir/srtm/grd/$stem.grdcut.grd $oDir/srtm/grd/$stem.grd

grd2xyz   $landmask.grd -S |\
	awk '{ printf "%.16lg %.16lg %.16lg %.16lg\n", $1, $2, $3, 0 }' \
			>! $ping.xyzi

# interpolate predicted to location of median ping; no nans here. grdtrack will keep sid 
# and add 5th column to our data.
#
# each row is xy of median, median, sid of median, predicted at xy of median
#
# calculate x, y, delta = median(ping)-predicted , sid
#
# normally we need diff of actual and predicted. if on land, put 0 in for depth and 
# predicted, so that surface won't generate wild excusions over land... 
#
# Surface ignores sid but why waste time adding it to file?
#
# interpolate delta (surface at full res grid)

if ($w < -180) then
	cat $ifile | awk '{ tmp=$1; if (tmp< +180 && tmp>='$w360') printf "%.16lg %.16lg %.16lg %.16lg\n",tmp-360,$2,$3,$4}' \
	>>! $ping.xyzi

	cat $ifile | awk '{ tmp=$1; if (tmp>=-180 && tmp<=  0) print $0 }' \
	>>! $ping.xyzi
else
	cat $ifile \
	>>! $ping.xyzi
endif

../bin/medianId -Q $inc $region $ping.xyzi | tee $ping.MedianId.xyzi |\
grdtrack -G$predicted $region -S |\
awk '{ tmp=$1; if (tmp>=+180) tmp=tmp-360; printf "%.16lg %.16lg %.16lg\n", tmp, $2, $3-$5}' |\
surface $region $inc $tension -G$ping.xyd.surface.grd	


cp $ping.MedianId.xyzi /tmp/ping.MedianId.xyzi
echo "ping count of `wc -l $ping.MedianId.xyzi`"

# add (actual - predicted) to predicted to get actual, at least where we have actual...

grdmath \
	$ping.xyd.surface.grd \
	$predicted \
	ADD = $polished.grd
	
grdcut $finalRegion $polished.grd \
	-G$polished.grdcut.grd
grdcut $finalRegion $landmask.grd \
	-G$landmask.grdcut.grd

# make sure that water is always negative, don't care about land because next step is to
# overlay SRTM land in mergeBathymetryTileWithSRTM.csh

# FIXME: if the -R is NOT an even multiple of "inc" then -R gets tweaked in polished, so tweak landmask to match
# grdedit -R$polished.grdcut.grd $landmask.grd
grdmath \
	$landmask.grdcut.grd \
	ISNAN 0 NAN -2 MUL \
	$polished.grdcut.grd \
	MIN = $oDir/grd/$stem.Bathymetry.grd

# create --->>>sid<<<---- grid. 
#
#   Use sid from above where there is a sid, 
#	but leave it Nan otherwise. (at least for the moment)
#
# Note we keep 1 deg pad in sid grid to keep xyz2grd from moaning,
# 	but we have to add/trim arctic latter, 
#		so this isn't a big deal.

# FIXME: xyz2grd reports multiple hits in some nodes, should not happen after blockmedian
cat $ping.MedianId.xyzi | \
	awk '{ printf "%.16lg %.16lg %05d\n", $1, $2, $4 }' | \
		xyz2grd $region $inc \
			-G$sidTmp.NaN.grd

# use grdmath to find sid with NaN, and it replace it with sid from
# predicted ("pred") sid grid. Predicted sid is either the arctic sid,
# or the plain old predicted sid.

grdmath \
	$sidTmp.NaN.grd \
	$oDir/predTile/$stem.sid.grd \
	AND = \
	$sidTmp.untrimmed.grd

grdcut $sidTmp.untrimmed.grd \
	-G$oDir/grd/$stem.sid.grd \
	$finalRegion


# all done 

rm $oDir/tmp/`basename $0`.$stem*
