#!/bin/bash
# Written by Collin Ferebee & Scott T. Marshall 3/25/2009
# Updated 4/9/2012 by Scott T. Marshall
#
# This script plots the a map of Southern California
# with a shaded DEM from the USGS NED dataset.
# On top of the map is plotted earthquakes from the 
# Lin et al (2007) relocated seismicity catalog.
# Data courtesy of the SCEC datacenter and was processed using getEQloc.pl
# 
#

#---------------------------------------------------------------------------------------------#
#---------------------------------  DEFINE SCRIPT VARIABLES  ---------------------------------#
#---------------------------------------------------------------------------------------------#

#Should I plot the DEM? 1 for yes 0 for no.
plotDEM=1
#Name the output file
map_file=SoCal_3arc.eps
#Define the W/E/S/N edges in degrees
area=-121.0/-114.0/32.0/37.0   
#The horizontal map width (c=cm i=in p=pixels)
scale=46.5i                      
#Define the lighting angle
azi=050
#How much should I move the plot over to center it on the page?
xShift=1.5i
yShift=1.25i

#If the .eps file exists, remove it
if [ -f $map_file ]
then
	rm $map_file
fi

#-----------------------------------------------------------------------------------------#
#---------------------------------  GMT SETUP VARIABLES  ---------------------------------#
#-----------------------------------------------------------------------------------------#

#set the GMT variable for eps paper size
gmtset PAPER_MEDIA Custom_49ix42i
#set the GMT variable for dots per inch (resolution)
gmtset DOTS_PR_INCH 300
#set the GMT variable for font size
gmtset ANNOT_FONT_SIZE_PRIMARY 20
#set the GMT variable for map scale height
gmtset MAP_SCALE_HEIGHT 0.25i


#----------------------------------------------------------------------------------------#
#---------------------------------  PROCESS GRID FILES  ---------------------------------#
#----------------------------------------------------------------------------------------#

#convert from the binary .bil file to a GMT .grd file. GMT will only plot .grd files
#Processing northwest corner
#xyz2grd NED_50238428.bil -G1arc/1.grd -R-121.0/-117.5/34.5/37.0 -I1c  -ZTLh -V -F
#give info for the newly created grid
#grdinfo 1arc/1.grd

#Processing northeast corner
#xyz2grd NED_33677582.bil -G1arc/2.grd -R-117.5/-114.0/34.5/37.0 -I1c  -ZTLh -V -F
#give info for the newly created grid
#grdinfo 1arc/2.grd

#Processing southwest corner
#xyz2grd NED_30379128.bil -G1arc/3.grd -R-121.0/-117.5/32.0/34.5 -I1c  -ZTLh -V -F
#give info for the newly created grid
#grdinfo 1arc/3.grd

#Processing southeast corner
#xyz2grd NED_32169894.bil -G1arc/4.grd -R-117.5/-114.0/32.0/34.5 -I1c  -ZTLh -V -F
#give info for the newly created grid
#grdinfo 1arc/4.grd

#downsample the .grd files to 3 arc seconds
#grdsample 1arc/1.grd -G3arc/1.grd -I3c
#grdsample 1arc/2.grd -G3arc/2.grd -I3c
#grdsample 1arc/3.grd -G3arc/3.grd -I3c
#grdsample 1arc/4.grd -G3arc/4.grd -I3c

#make intensity files for to figure out the range of sigma values for all grd files
#grdgradient 3arc/1.grd -A$azi -G3arc/1_int.grd -Nt -V
#grdgradient 3arc/2.grd -A$azi -G3arc/2_int.grd -Nt -V
#grdgradient 3arc/3.grd -A$azi -G3arc/3_int.grd -Nt -V
#grdgradient 3arc/4.grd -A$azi -G3arc/4_int.grd -Nt -V

#Now that I know the sigmas from these grids, apply the avg sigma value using the -Nt[amp/sigma/offset] option.
#This should give uniform shading for all grid tiles.
#grdgradient 3arc/1.grd -A$azi -G3arc/1_int.grd -Nt1/14000/0
#grdgradient 3arc/2.grd -A$azi -G3arc/2_int.grd -Nt1/14000/0
#grdgradient 3arc/3.grd -A$azi -G3arc/3_int.grd -Nt1/14000/0
#grdgradient 3arc/4.grd -A$azi -G3arc/4_int.grd -Nt1/14000/0


#--------------------------------------------------------------------------------------------#
#---------------------------------  START PLOTTING THE MAP  ---------------------------------#
#--------------------------------------------------------------------------------------------#

if [ $plotDEM -eq 1 ]
then
	#plot topography and print a progress report
	echo printing 1arc/1.grd...
	grdimage 3arc/1.grd -X$xShift -Y$yShift -R$area -JM$scale -CCaliColorMap.cpt -I3arc/1_int.grd -P -K >> $map_file 
	echo printing 1arc/2.grd...
	grdimage 3arc/2.grd -R -JM -CCaliColorMap.cpt -I3arc/2_int.grd -O -K >> $map_file 
	echo printing 1arc/3.grd...
	grdimage 3arc/3.grd -R -JM -CCaliColorMap.cpt -I3arc/3_int.grd -O -K >> $map_file 
	echo printing 1arc/4.grd...
	grdimage 3arc/4.grd -R -JM -CCaliColorMap.cpt -I3arc/4_int.grd -O -K >> $map_file 
	
	#plot basemap mark every degree
	psbasemap -R -JM -Bf60ma60m/f60ma60m:horizontal:WSne -O -K --ANNOT_FONT_SIZE_PRIMARY=28 --FRAME_WIDTH=0.20i --ANNOT_OFFSET_PRIMARY=0.25i >> $map_file
	
else
	#plot basemap mark every degree
	psbasemap -X$xShift -Y$yShift -R$area -JM$scale -Bf60ma60m/f60ma60m:horizontal:WSne -P -K --ANNOT_FONT_SIZE_PRIMARY=28 --FRAME_WIDTH=0.20i --ANNOT_OFFSET_PRIMARY=0.25i >> $map_file

fi

#Plot the coast and make the water blue and the land white
pscoast -R -JM -Na/4 -W4 -S102/153/204 -Df -O -K >> $map_file

#plot seismicity as red circles
echo plotting earthquakes...
psxy Lin_et_al.mag -R$area -JM$scale -Sc0.20c -W1 -G190/30/45 -H1 -O -K >> $map_file

#plot faults
psxy fault.lin -R -JM -H0 -m -W8/0/0/0 -O -K >> $map_file

#plot beachballs
psmeca HardebeckShearerFocalMech.txt -R -JM -Sa1/20/100 -W1p -H1 -O -K >> $map_file
#psmeca neicSoCal.mt2 -R -JM -Sa1/20/0 -W1p -H1 -O -K >> $map_file
#psmeca HarvardMT.txt -R -JM -Sm1/20/0 -W1p -H1 -O -K >> $map_file

#Plot significant historical earthquakes as stars
psxy historicalEQ.txt -R -JM -H1 -W4 -G210/172/043 -Sa0.5i -O -K >> $map_file
pstext historicalEQ.txt -R -JM -H1 -D0/0.35i -G210/172/043 -S4 -O -K >> $map_file

#plot a legend 
pslegend -R -JM -Dx0.5i/0.5i/13.75i/8.6i/BL -G255/255/255 -F -O -K << END >> $map_file
G 0.25i
H 36 1 Southern California Earthquakes > M@-w@-0  (1981-2005)
G 0.15i
H 28 1 405,712 Total Events in 25 Years; 16,229 Events/Average Year
G 0.25i
H 24 1 Compiled by Collin L. Ferebee and Scott T. Marshall
G 0.15i
H 24 1 Department of Geology, Appalachian State University
G 0.25i
S 1.0i c 0.2c 190/30/45 0.5 1.5i @:19:Earthquake Epicenters from Lin et al. (2007)@::
G 0.25i
S 1.0i c 0.0001i 255/255/255 0.0 1.5i @:19:Focal Mechanisms > M@-W@-5.0 from Hardebeck and Shearer (2003)@::
G 0.25i
S 1.0i a 0.5i 210/172/43 4.0 1.5i @:19:Large Historical Epicenters from Stover and Coffman (1993)@::
G 0.25i
S 1.0i y 0.4i 0/0/0 8.0 1.5i @:19:Fault Traces from Jennings (1992)@::
G 0.50i
M 242.5 34.5 100k+u f 
G 0.35i
M 242.5 34.5 100m+u f 
G 0.35i
B CaliColorMap.cpt 0.5i 0.5i -B:"Elevation Above Sea Level (meters)": --LABEL_FONT_SIZE=20 --LABEL_FONT=0 --ANNOT_FONT_SIZE_PRIMARY=20
END

#Plot the North Arrow and Directional Rose
pscoast -R -JM -W0 -Df -Tfx12.25i/5.25i/2.25i/2 -O -K --HEADER_FONT_SIZE=26 --HEADER_OFFSET=0.2c >> $map_file

#plot a beachball for the legend
psmeca -R -JM -Sa1/20/100 -W1p -O -K <<END>> $map_file
-120.767 32.714 24.0 320 90 180 5.0 0 0
END

#This is what you need to finish off an .eps file correctly
echo showpage >> $map_file


#------------------------------------------------------------------------------------------------------------------#
#---------------------------------  DISPLAY THE MAP AND CONVERT TO OTHER FORMATS  ---------------------------------#
#------------------------------------------------------------------------------------------------------------------#

#display $map_file in a program that can read .eps files.
#gv $map_file &
#acrobat $map_file &
#illustrator $map_file &

#make a pdf, jpg, and transparent png from the eps files
echo converting to pdf
ps2raster $map_file -Tf
#Now open the pdf with acrobat
acrobat SoCal_3arc.pdf &
echo converting to jpg
ps2raster $map_file -Tj
echo converting to png
ps2raster $map_file -TG

#print a done message and ring the system bell
echo -e "Done! \a" 
exit

#--------------------------------------------------------#
#--------------------------------------------------------#