#!/bin/bash


gmt gmtset FORMAT_GEO_MAP DD
gmt gmtset MAP_FRAME_TYPE plain
gmt gmtset FONT_TITLE 24p
gmt gmtset MAP_TITLE_OFFSET -0.4c 
gmt gmtset FONT_LABEL 16p
gmt gmtset MAP_LABEL_OFFSET 0.3c
gmt gmtset FONT_ANNOT_PRIMARY 12p

#Set variables for calculations and plotting
bigregion="-R318490/320370/3769236/3771478"
region="-R319200/319650/3769950/3770750"
insetRegion="-R-116/-98/30/41"
insetProj="-JM?i"

zoominc="-I0.25/0.25"
biginc="-I10/10"
zoomproj="-Jx1:4000/1:4000"
bigproj="-Jx1:24000/1:24000"

#keep_open="-K -V"

#add="-O -K -Vd"

#close="-O -V"

#psFile="BCDFile.ps"

#pdfFile="BCDFile.pdf"

#path2ElevFile="/Users/spri902/Desktop/SubTER_BCD/BCD_DEM"

verb="-Vd"

#view="-p120/25"

coastmisc="-Ggrey -Sblue -Na/1p"
xymisc="-Sa0.2c -W0.1p -G255/255/0 -:"

# blockmean the xyz file to preprocess data for surface function
gmt blockmean bcd_dem.xyz $bigregion $biginc > bcd_dem_blk.grd

# call surface function to grid the table data using adjustable tension continuous curvature splines
gmt surface bcd_dem_blk.grd -Gbcd_grid.nc $bigregion $biginc 
gmt surface bcd_dem_blk.grd -Gbcd_gridcut.nc $region $zoominc 

#Cut the grid file to the size i need
#gmt grdcut bcd_dem_blk.grd -Gbcd_dem_blkcut.grd $region

#Set the sun shadow to illuminate the data (sun is 90degrees north and at 50deg angle from horizon)
gradientmisc="-A20" 
gmt grdgradient bcd_grid.nc -Gbcd_topo.nc $gradientmisc
gmt grdgradient bcd_gridcut.nc -Gbcd_topocut.nc $gradientmisc

gmt grdhisteq bcd_topo.nc -Gbcd_topohist.nc -N 
gmt grdhisteq bcd_topocut.nc -Gbcd_topohistcut.nc -N 
gmt grdmath bcd_topohist.nc 5.5 DIV = bcd_toponorm.nc
gmt grdmath bcd_topohistcut.nc 5.5 DIV = bcd_toponormcut.nc

# Create color pallet 
#colormisc = "-Chaxby -D -V -T1524/2126/1"
#gmt makecpt $colormisc > bcd_topo.cpt
gmt makecpt -Cdem2 -D -V -T1524/2126/0.1 > bcd_topo.cpt
gmt makecpt -Cterra -D -V -T100/2200/0.1 > region.cpt
gmt makecpt -Cgrey -D -V -T1524/2126/0.1 >bcd_grey.cpt
# make subplots w top fig of large scale DEM w inset and bottom fig zoomed in on BCD


gmt begin bcdmap
	gmt set MAP_ANNOT_OBLIQUE 34 FORMAT_GEO_MAP ddd:mmmF
	gmt subplot begin 2x2 -T"Blue Canyon Dome Site Figure" -Fs5c,20c/10c,10c -M1c -Ba -BWSne $verb
		#upper figure column1
		gmt subplot set 0,0
			gmt grdimage @earth_relief_03s -R-111/111/-111/111+uk -JA107:14W/34:04N/? -Cbcd_topo.cpt
			gmt inset begin -DjTL+w2.0c+o0.2c -M0 $verb
				gmt basemap $insetRegion $insetProj -Bwesn $verb
				gmt coast $insetRegion $insetProj $coastmisc $verb
				gmt plot $insetRegion $insetProj $xymisc <<-EOF 
				34.059 -106.956
				EOF
			gmt inset end
		# upper figure column2
		gmt subplot set 0,1
			gmt grdimage bcd_grid.nc $bigproj $bigregion -Ibcd_toponorm.nc -Cbcd_grey.cpt $verb
			
		#lower figure
		gmt subplot set 1,0
			gmt grdimage bcd_gridcut.nc $projection $region -Ibcd_toponormcut.nc -Cbcd_topo.cpt $verb

			# plot well head points
			wellmisc="-Sc0.1c -Gred"
			gmt plot wells.txt $region $projection $wellmisc
			# This adds the north arrow to the plot
			gmt text $region $projection <<-EOF 
			319700 3770750 36 90 34  BL \343
			EOF
	gmt subplot end
gmt end show

echo "Finished"
rm bcd_topohist.nc
rm bcd_toponorm.nc
rm bcd_topo.nc
rm bcd_grid.nc