#!/bin/sh
#	DRFT06RR GRIDDING OF MULTIBEAM DATA
#
# This will grid ascii data for a given box id.  The box areas
# are defined in the file ../drft06rr_boxes.lis which contains the
# coordinates of the box center, and we assume each box is 6x6
# degrees.  Preliminary analysis of the raw xyz data suggest we
# should use no smaller grid spacing that ~100m.  Such as spacing
# is roughly 0.001 degrees (which is ~110 m).  We will first
# blockmean the data down to 0.001 degrees then run nearneighbor
# on the binned averages with a 0.003 degree seach radius (~330m).
# We will automatically determine the actual -R extent of the data 
# inside the box and use that -R for the gridding.  Once the multi-
# beam data has been gridded we extract the Sandwell/Smith predicted
# bathymetry for the same area, resample at desired grid spacing,
# and use grdblend to blend the two grids into a final box grid.
# We will use gridline registration to facilitate blending with
# Easter/Salas y Gomez grids.
# We also make a 0.003 degree grid
#
# Paul Wessel, Nov. 10, 2001 on the Revelle
#
# Dec 13,2001: There was one day of missing MB data that we had to pull
# out from the Simrad system.  Its format is slightly different.  The
# raw xyz file is called 20011211-All.xyz and is play xyz which I then
# converted once and for all to 20011211-All.xyz.bin and then bzip2'ed the
# xyz file for arcdhive purposes.  The rest of the script looks for all
# *.bin files so both the huge MB file and this missing day will be used.
if [ $# -eq 0 ]; then
	echo "usage: drft06rr_gridxyz box [-p] [-f]" >&2
	echo "  -p will plot the grid and generate a raster for xv"
	echo "  -f will force regridding if grid already exists"
	exit
fi
cd gridding
# Define processing parameters:
Dfine=0.001
Dcourse=0.003
Srad=0.003
PAD=15k
GEOWARE2M=5
quiet=1
do_view=0
force=0
if [ $# -ge 2 ] && [ $2 = "-p" ]; then
	do_view=1
	if [ $# -eq 3 ] && [ $2 = "-f" ]; then
		force=1
	fi
elif [ $# -eq 2 ] && [ $2 = "-f" ]; then
	force=1
fi
if [ $quiet -eq 0 ]; then
	babble="-V"
#	babble=
fi
box=$1
if [ -f grids/drft06rr.box$box.mb+ss.${Dfine}.grd ] && [ $force -eq 0 ]; then
	echo "drft06rr_gridxyz: Grid for box $box exists, use -f to force regridding" >&2
	exit
fi
HW=3
Wbox=`grep ^$box ../drft06rr_boxes.lis | awk '{print $2}'`
Ebox=`grep ^$box ../drft06rr_boxes.lis | awk '{print $3}'`
Sbox=`grep ^$box ../drft06rr_boxes.lis | awk '{print $4}'`
Nbox=`grep ^$box ../drft06rr_boxes.lis | awk '{print $5}'`
Rbox="$Wbox/$Ebox/$Sbox/$Nbox"

# Need output format to have decimal seconds since 0.001 degree = 3.6 arc secs
gmtset OUTPUT_DEGREE_FORMAT ddd:mm:ss.xF GRIDFILE_SHORTHAND true

# 1. Convert the ascii xyz files to more efficient native double precision binary files
#    Note that data were passed to me as a *.zip file (on revell in /Work.Disk.08/DRFT06RRxyz)
#    which I unzip -a to deal with the CR/LF issue.  The records contain comma-separater x,y,z
#    triplets: Here is one example:
#     080-10-09.946W,019-55-08.472S,3005.330
#    which explains why we use tr to replace the - with : before letting GMT see it.

echo "drft06rr_gridxyz: 1. Convert xyz->bin..." >&2
cd raw
nzip=`(ls *.zip 2> /dev/null) | wc -l`
if [ $nzip -gt 0 ]; then
	for i in *.zip; do
		echo "  unzip -a $i"
		unzip -a $i
		rm -f $i
	done
fi
ntxt=`(ls *.txt 2> /dev/null) | wc -l`
if [ $ntxt -gt 0 ]; then
	for i in *.txt; do
		if [ ! -f $i.bin ]; then
			echo "  convert $i to GMT double precision binary file $i.bin"
			tr '-' ':' < $i | gmtconvert -bod -f0x,1y > $i.bin
			echo "  bzip2 $i to $i.bz2"
			rm -f $i.bz2
			bzip2 -9 $i
		fi
	done
fi
nbin=`(ls *.bin 2> /dev/null) | wc -l`
if [ $nbin -eq 0 ]; then
	echo no data available yet
	exit
fi
cd ..

# 2. Keep only the points inside the current 6x6 box

echo "drft06rr_gridxyz: 2. Extract subset..." >&2
gmtselect -R$Rbox raw/*.bin -bi3d -bod -f0x,1y $babble > use.bb
size=`ls -l use.bb | awk '{print $5}'`
if [ $size -eq 0 ]; then
	echo no data available for box $box yet
	exit
fi

# Determine region to the nearest 2' since we must coregister with predicted bathy grids

info=`minmax use.bb -I2m -bi3d -f0x,1y -C`
# Add 20m in all directions for padding
# This gives an -R useful for pulling out 2m grids from grdraster
echo $info | awk '{print $1}' > tmp
W2m=`gmtmath tmp -f0x -Ca 20 60 DIV SUB =`
echo $info | awk '{print $2}' > tmp
E2m=`gmtmath tmp -f0x -Ca 20 60 DIV ADD =`
echo $info | awk '{print $3}' > tmp
S2m=`gmtmath tmp -f0y -Ca 20 60 DIV SUB =`
echo $info | awk '{print $4}' > tmp
N2m=`gmtmath tmp -f0y -Ca 20 60 DIV ADD =`
R2m="-R$W2m/$E2m/$S2m/$N2m"
info=`minmax use.bb -I${Dcourse} -bi3d -f0x,1y -C`
# Add 4.5m in all directions for padding when gridding, but make sure we do not exceed the 6x6 box
# This gives a -R useful for gridding at 0.001 degrees
echo $info | awk '{print $1}' > tmp
W=`gmtmath tmp -f0x -Ca 4.5 60 DIV SUB $Wbox MAX $Dcourse DIV FLOOR $Dcourse MUL =`
echo $info | awk '{print $2}' > tmp
E=`gmtmath tmp -f0x -Ca 4.5 60 DIV ADD $Ebox MIN $Dcourse DIV CEIL $Dcourse MUL =`
echo $info | awk '{print $3}' > tmp
S=`gmtmath tmp -f0y -Ca 4.5 60 DIV SUB $Sbox MAX $Dcourse DIV FLOOR $Dcourse MUL =`
echo $info | awk '{print $4}' > tmp
N=`gmtmath tmp -f0y -Ca 4.5 60 DIV ADD $Nbox MIN $Dcourse DIV CEIL $Dcourse MUL =`
R="-R$W/$E/$S/$N"
# Add 5.4m in all directions for padding when blockmean.
# This gives a -R useful for blockmeaning at 0.001 degrees since we want nearneighbor to have values
# outside its -R for isotropic solution there too.
echo $info | awk '{print $1}' > tmp
Wb=`gmtmath tmp -f0x -Ca 5.4 60 DIV SUB =`
echo $info | awk '{print $2}' > tmp
Eb=`gmtmath tmp -f0x -Ca 5.4 60 DIV ADD =`
echo $info | awk '{print $3}' > tmp
Sb=`gmtmath tmp -f0y -Ca 5.4 60 DIV SUB =`
echo $info | awk '{print $4}' > tmp
Nb=`gmtmath tmp -f0y -Ca 5.4 60 DIV ADD =`
Rb="-R$Wb/$Eb/$Sb/$Nb"
if [ $quiet -eq 0 ]; then
	echo "  ${Dfine} x ${Dfine} gridding in region $R"
	echo "  ${Dfine} x ${Dfine} blockmean in region $Rb"
	echo "  2m x 2m grdraster in region $R2m"
fi

# 3. Get predicted bathymetry, make positive down, and resample at final grid resolution

echo "drft06rr_gridxyz: 3. Extract and prep SS background..." >&2
grdraster $GEOWARE2M $R2m -Gss.i2 $babble
grdsample ss.i2 -Gss.i2 -T
# Flip sign temporarily to match MB data, then mask out data far from MB
grdmath $babble ss.i2 NEG = ss.i2
grdsample ss.i2 -Gss${Dfine}.i2 -I${Dfine} $Rb $babble -f0x,1y
rm -f ss.i2 tmp

# 4. Blockmean and grid the mb data using nearnest neighbor method

echo "drft06rr_gridxyz: 4. Block average..." >&2
blockmean $Rb -I${Dfine} -bi3d use.bb -bod -f0x,1y $babble > use.b${Dfine}
rm -f use.bb

# 5. Use nearneighbor to grid the mb data

echo "drft06rr_gridxyz: 5. Grid data..." >&2
nearneighbor $R -I${Dfine} -S$Srad -Ggrids/drft06rr.box$box.mb.${Dfine}.grd -bi3d -f0x,1y $babble use.b${Dfine}

# 6. Blend the mb and the predicted data to yield final grid

echo "drft06rr_gridxyz: 6. Blend mb and SS..." >&2
cat << EOF > blend.d
grids/drft06rr.box$box.mb.${Dfine}.grd  $R  1
ss${Dfine}.i2         $R  0.01
EOF
Rblend=$R
if [ $box = F ]; then	# blend in old 1993 Easter grid
	echo "../data/F.0.001.i2	-R-105.001/-103/-28/-24.5	0.1" >> blend.d
	Rblend="-R-105/-99/-28.001/-24.5"
	grdblend -Gt.grd -I${Dfine} $Rblend -f0x,1y $babble << EOF
grids/drft06rr.box$box.mb.${Dfine}.grd  $R  1
EOF
	\mv -f t.grd grids/drft06rr.box$box.mb.${Dfine}.grd	
	grdraster $GEOWARE2M -R-105/-99/-28/-24.5 -Gss2.i2 $babble
	grdsample ss2.i2 -Gss2.i2 -T
	# Flip sign temporarily to match MB data, then mask out data far from MB
	grdmath $babble ss2.i2 NEG = ss2.i2
	grdsample ss2.i2 -Gss2${Dfine}.i2 -I${Dfine} -R-105/-99/-28/-24.5 $babble -f0x,1y
	rm -f ss2.i2 tmp
	echo "ss2${Dfine}.i2 -R-105/-99/-28/-24.5 0.01" >> blend.d
fi
if [ $box = G ]; then	# blend in old 1993 Easter grid
	echo "../data/G.0.001.i2	-R-110.001/-104.001/-29/-24.5	0.1" >> blend.d
	Rblend="-R-110.001/-104.001/-29/-24.5"
	grdblend -Gt.grd -I${Dfine} $Rblend -f0x,1y $babble << EOF
grids/drft06rr.box$box.mb.${Dfine}.grd  $R  1
EOF
	\mv -f t.grd grids/drft06rr.box$box.mb.${Dfine}.grd	
	grdraster $GEOWARE2M -R-110/-104/-29/-24.5 -Gss2.i2 $babble
	grdsample ss2.i2 -Gss2.i2 -T
	# Flip sign temporarily to match MB data, then mask out data far from MB
	grdmath $babble ss2.i2 NEG = ss2.i2
	grdsample ss2.i2 -Gss2${Dfine}.i2 -I${Dfine} -R-110/-104/-29/-24.5 $babble -f0x,1y
	rm -f ss2.i2 tmp
	echo "ss2${Dfine}.i2 -R-110/-104/-29/-24.5 0.01" >> blend.d
fi
grdblend -Ggrids/drft06rr.box$box.mb+ss.${Dfine}.grd -I${Dfine} $Rblend -Z-1 -f0x,1y $babble blend.d
rm -f use.b${Dfine} use.bb ss.i2 ss${Dfine}.i2 blend.d

# 7. Downgrade to ${Dcourse} degree grids using grdfilter

echo "drft06rr_gridxyz: 7. Create lower-resolution grid..." >&2
grdfilter grids/drft06rr.box$box.mb+ss.${Dfine}.grd -Ggrids/drft06rr.box$box.mb+ss.${Dcourse}.grd -Fg0.006 -D0 -I${Dcourse} $babble 
grdfilter grids/drft06rr.box$box.mb.${Dfine}.grd -Ggrids/drft06rr.box$box.mb.${Dcourse}.grd -Fg0.006 -D0 -I${Dcourse} $babble 

# 8. Encode data origins by letting all SS depths be in integer meters and all
#    MB soundings be rounded to nearest 5 cm, so z = N * 0.05 m + 0.025 m
#    We can then use is_MB = (abs(fmod(z,1)) > 0.02) to pull out the MB data.

echo "drft06rr_gridxyz: 8. Encode origins of data..." >&2
grdmath grids/drft06rr.box$box.mb.${Dcourse}.grd ISNAN = mb_0_ss_1.bit
grdmath grids/drft06rr.box$box.mb+ss.${Dcourse}.grd 20 MUL FLOOR 20 DIV 0.025 ADD mb_0_ss_1.bit 1 SUB NEG MUL grids/drft06rr.box$box.mb+ss.${Dcourse}.grd RINT mb_0_ss_1.bit MUL ADD = grids/drft06rr.box$box.mb+ss.${Dcourse}.grd
grdedit grids/drft06rr.box$box.mb+ss.${Dcourse}.grd -Ddegree/degree/meter/1/0/"DRFT06RR [Nazca Ridge 2001, PI: D. Naar] Box A Merged Bathymetry (downgraded)"/"For multibeam only: grdmath drft06rr.box$box.mb+ss.${Dcourse}.grd 1 FMOD ABS 0.02 GT 0 NAN drft06rr.box$box.mb+ss.${Dcourse}.grd MUL = mb.grd"
grdmath grids/drft06rr.box$box.mb.${Dfine}.grd ISNAN = mb_0_ss_1.bit
grdmath grids/drft06rr.box$box.mb+ss.${Dfine}.grd 20 MUL FLOOR 20 DIV 0.025 ADD mb_0_ss_1.bit 1 SUB NEG MUL grids/drft06rr.box$box.mb+ss.${Dfine}.grd RINT mb_0_ss_1.bit MUL ADD = grids/drft06rr.box$box.mb+ss.${Dfine}.grd
grdedit grids/drft06rr.box$box.mb+ss.${Dfine}.grd -Ddegree/degree/meter/1/0/"DRFT06RR [Nazca Ridge 2001, PI: D. Naar] Box A Final Merged Bathymetry"/"For multibeam only: grdmath drft06rr.box$box.mb+ss.${Dfine}.grd 1 FMOD ABS 0.02 GT 0 NAN drft06rr.box$box.mb+ss.${Dfine}.grd MUL = mb.grd"
rm -f mb_0_ss_1.bit

# 9. Get intensity grid for various grdimage plotting

echo "drft06rr_gridxyz: 9. Calculate intensity grids..." >&2
grdgradient grids/drft06rr.box$box.mb+ss.${Dfine}.grd -Ggrids/drft06rr.box$box.mb+ss.${Dfine}.int.grd -Nt1.2 -A45 $babble 
grdgradient grids/drft06rr.box$box.mb+ss.${Dcourse}.grd -Ggrids/drft06rr.box$box.mb+ss.${Dcourse}.int.grd -Nt1.2 -A45 $babble 

if [ $do_view -eq 1 ]; then
# 10. Make a bathymetric plot
	echo "drft06rr_gridxyz: 10. Make plot..." >&2
	echo "  Make PostScript plot" >&2
	height=`echo $Ebox $Nbox | mapproject -R$Rbox -JM6 | cut -f2`
	y=`gmtmath -Q $height 1.0 ADD =`
	makecpt -Crainbow -T-6000/0/500 -Z > t.cpt
	grdimage -R$Rbox -JM6 -B1g1WSne grids/drft06rr.box$box.mb+ss.${Dcourse}.grd -Igrids/drft06rr.box$box.mb+ss.${Dcourse}.int.grd -Ct.cpt -P -Y0.5 $babble > t.ps
	rm -f t.cpt
	echo "  Convert to Sun rasterfile" >&2
	ps2ras -w 7.1 -h $y -d 24 -r 300 t.ps t.ras
	echo "  View rasterfile" >&2
	xv t.ras &
fi

# We can now delete the mb only grids

rm -f grids/drft06rr.box$box.mb.${Dfine}.grd grids/drft06rr.box$box.mb.${Dcourse}.grd
echo "drft06rr_gridxyz: Done" >&2
echo "  Grid is grids/drft06rr.box$box.mb+ss.${Dfine}.grd" >&2
echo "  Intensity is grids/drft06rr.box$box.mb+ss.${Dfine}.int.grd)" >&2
cd ..