#! /bin/csh

if ($#argv < 8) then
	echo "usage: $0 name W E S N scale incInSec oDir options "
	echo "  example: $0 e020n40 20 60 -010 40 .25 -P .."
	exit
endif

# 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

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 incInSec = $1;	shift
set oDir	= $1;   shift
set options	= "$*"

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

set w = `echo $w | awk '{ if ($1 > 180) print $1-180.0; else print $1 }'`
set s = `echo $s | awk '{ if ($1<-89) tmp = -89; else tmp = $1; printf "%.16lg\n", tmp }'`
set n = `echo $n | awk '{ if ($1>89) tmp = 89; else tmp = $1; printf "%.16lg\n", tmp }'`

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


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


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 at $incInSec':WSen
set tickInfo = -Ba10f1g0:."$incInSec sec":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 $oDir/grd/$stem.Bathymetry.srtm.grd -G$tmpStem.cut.grd $region -V

# debug that darn New Siberian Island issue...
#grdcut $oDir/grd/$stem.Bathymetry.grd -G$tmpStem.cut.grd $region -V
#grdcut $oDir/srtm/grd/$stem.grd -G$tmpStem.cut.grd $region -V
#grdcut $oDir/predTile/$stem.grd -G$tmpStem.cut.grd $region -V
#grdcut $oDir/tmp/jjMedian.e140n90.median.xyd.surface.grd -G$tmpStem.cut.grd $region -V
#grdcut $oDir/tmp/foo.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
# 
# grdcontour $tmpStem.filtered.grd $projection -C$intervalsFile \
# 		-Wthinnest/$coastlineColor -O \
# 			>> $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

