#! /bin/csh -f

# script to postscript plot of grd file.
# that's right: 140+ lines to plot a grid. 
# GMT is a time saver.
#
# magic and hidden files: Another of the wonders of GMT

rm -f .gmt*

# default color space is rgb. For some reason that is fucked up.
# default longitude is 0-360, make it +/-180

gmtset	COLOR_MODEL		hsv
gmtset	PLOT_DEGREE_FORMAT	ddd:mmF

gmtset PAPER_MEDIA archD
gmtset PAPER_MEDIA a4

if ($#argv<8 || $#argv>9) then
	echo "usage: $0 name W E S N scale pings pingInc options"
	echo "  example: $0 arabia 20 60 -010 40 .25 e020n40.xyzi 30c -P"
	echo " -P for portrait is optional"
	exit
endif

set stem	= $1;   shift
set w		= $1;   shift
set e		= $1;   shift
set s		= $1;   shift
set n		= $1;   shift
set inchPerDeg	= $1;   shift
set pings	= $1;   shift
set pingInc	= -I$1;	shift
set options	= "$*"

set ofile	= ../plots/$stem.ps
set tmpStem	= /tmp/`basename $0`.$stem

set w = `echo $w | awk '{ if ($1 > 180) print $1-180.0; else print $1 }'`

set region = "-Rg/$w/$e/$s/$n"
set coastlineColor = 40/40/40
set colorTable = /sw/share/gmt/cpt/GMT_topo.cpt
set colorTable = ../plotHitMap/topo.cpt
set projection = "-Jm$inchPerDeg"


set intervalsFile = ../com/intervals_none
set intervalsFile = ../com/intervals_0
set intervalsFile = ../com/intervals_3000


set azimuth = -A315
set azimuth = '-A0/270 -N1.0'
set azimuth = '-E315/15 -Nt1.0'
set azimuth = '-E315/90 -Nt0.6'
set azimuth = '-A0/270 -N10.0'

#set vertExhaggeration1 = 20000.
#set vertExhaggeration2 = 8000.


# grdfilter "distance" connotes pixel type and units for filter width
# grdfilter width set width of filter, like how many pixels get averaged.
# grdfilter type is gaussian, like mean, median, gaussian weighted filter.
#
# "D2" means grid is in degrees, filter width is in km 
# filter type is Gaussian, with width of 0.25 km
#
# grdfilter: CONTOUR ONLY is done at 2 minute spacing = 1/30 of a degree.

set filterDistance = -D2
set filterType = g
set filterWidth = 0.25
set grdFilterRule = -F$filterType$filterWidth
set outputInc = -I.03333333333333333

# annotate only W and S sides but draw east and north sides too
# Why is this always such a pain in the ass?

set tickInfo = -B1g1WSen
set tickInfo = "-Ba2g1/a1g1:."$stem":WS"
set tickInfo = -Ba10f1g0:."$stem":WSen

# left bottom corner of plot 1.5" right and 1.00" above bottom/left of page.

set xPlotOrigin = -X1.5
set yPlotOrigin = -Y1.00




# grdimage has many bugs, not least is crashing 
# when faced with mercator projections from a grd 
# that includes pole even if region given excludes pole.
#
# cut out region of interest and spoon feed grdimage.

grdcut ../grd/$stem.Bathymetry.srtm.grd -G$tmpStem.cut.grd $region -V


# shade topography 
# image dirt with more exaggeration than bathymetry? No.

# plan a) use grdgradient as it was intended.

grdgradient $tmpStem.cut.grd $azimuth -G$tmpStem.grad.grd

# plan b) equalize gradient explictly

#grdhisteq -V $tmpStem.grad.grd -G$tmpStem.grad.grd -N1.0

# plan c) and then there is Dave's way

#grdmath -V $tmpStem.grad.grd $vertExhaggeration1 DIV = $tmpStem.grad.grd



# make image

grdimage $tmpStem.cut.grd -I$tmpStem.grad.grd -C$colorTable \
	$projection $tickInfo $xPlotOrigin $yPlotOrigin $options \
	-fg -K >! $tmpStem.ps


# turn on clipping, and use pscoast to "clip dry areas", 
# render land (presumable with more vertical exaggeration),
# then turn off clipping.
#
#pscoast $projection $region -Dh -Gc -O -K >> $tmpStem.ps
#grdimage $tmpStem.cut.grd -I$tmpStem.grad2.grd -C$colorTable \
#	$projection $options -O -K >> $tmpStem.ps
#pscoast $projection $region -Q -O -K >> $tmpStem.ps


# image is used at full resolution, but contour is made using 
# decimated image so contour is smooth(er)

grdfilter $tmpStem.cut.grd $region $filterDistance $grdFilterRule \
		-G$tmpStem.filtered.grd $outputInc



# use -Q option so xy returned is xy at median z

../bin/medianId $pings $pingInc $region -F -Q -V -fg \
	|  awk '{print $1, $2, $3}' \
	|  sort -n -k1,2 \
	>! $tmpStem.pings.xy

set pingSymbolSize = 0.01c

psxy $tmpStem.pings.xy $projection $region -K -O -Sa$pingSymbolSize -fg >> $tmpStem.ps



grdcontour $tmpStem.filtered.grd -C$intervalsFile $projection -O \
		-Wthinnest/$coastlineColor \
			>> $tmpStem.ps

# All done

/bin/mv -f $tmpStem.ps $ofile
/bin/rm -f $tmpStem.cut.grd $tmpStem.grad1.grd $tmpStem.grad2.grd $tmpStem.filtered.grd $tmpStem.grad.grd $tmpStem.pings.xy