#!/bin/bash
#
# EastUS.gmt
# Written by Julia T Irizarry and Scott T. Marshall 8/2011
# 
# This script plots a shaded DEM of the eastern US along with
# the entire NEIC earthquake epicenter catalog from 1973-2011 (8/22/2011)
# The script also plots two insets of the New Madrid and Middleton Place-Summerville
# seismic zones.
#
#

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

#Should I plot the DEM? 1 for yes 0 for no.
plotDEM=1
#define the plot file
file=EastUS.eps
#define map lon/lat range
area=-95/-65/24/48
#The horizontal map width (c=cm i=in p=pixels)
scale=40.1i
#define two water colors
waterColor=146/182/221
waterColor2=210/230/250
#waterColor2=185/206/234
#define EQ dot red color
redColor=225/27/82
#define EQ dot blue color
blueColor=15/10/170
#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 $file ]
then
	rm $file
fi

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

#set the GMT variable for eps paper size
gmtset PAPER_MEDIA=Custom_42.5ix42i
#change axis font size
gmtset ANNOT_FONT_SIZE_PRIMARY=32p
#change the axis label offset
gmtset ANNOT_OFFSET_PRIMARY=0.25i
#change map scale height
gmtset MAP_SCALE_HEIGHT=0.25i
#change frame width
gmtset FRAME_WIDTH=0.25i
#change tick length
gmtset TICK_LENGTH=0.2
#change to fancy basemap
gmtset BASEMAP_TYPE=fancy
#set the GMT variable for frame width of the legend
gmtset FRAME_PEN=3p
#make a color palette file
#makecpt -CCharlestonColor.cpt -T0/540/90 -Z > EasternUS.cpt
#makecpt -CCaliColorMap.cpt -T0/1500/250 -Z > ColorTopo.cpt


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

#---only run once (makes grid file)
#convert the USGS binary *.bil files to a *.grd file some files may require pixel node registration (-F)
#xyz2grd n36w077.bil -G1.grd -R-77/-76/36/37  -I3.0c -ZTLh -V

#---only run once (makes the intensity file)
#make file for shadows given a color pallette file and a lighting direction
#echo "Shading U.S. DEM Tile 1/4..."
#grdgradient EastUS_DEM/1_6arc.grd -A270 -GEastUS_DEM/1_6arc_int.grd -V -Nt1/1000/0
#echo "Shading U.S. DEM Tile 2/4..."
#grdgradient EastUS_DEM/2_6arc.grd -A270 -GEastUS_DEM/2_6arc_int.grd -V -Nt1/1000/0
#echo "Shading U.S. DEM Tile 3/4..."
#grdgradient EastUS_DEM/3_6arc.grd -A270 -GEastUS_DEM/3_6arc_int.grd -V -Nt1/1000/0
#echo "Shading U.S. DEM Tile 4/4..."
#grdgradient EastUS_DEM/4_6arc.grd -A270 -GEastUS_DEM/4_6arc_int.grd -V -Nt1/1000/0


#EAST US DEM DOWNSAMPLING...
#Downsample the grid files
#grdsample EastUS_DEM/1_6arc.grd -GEastUS_DEM/1_12arc.grd -I12c
#grdsample EastUS_DEM/2_6arc.grd -GEastUS_DEM/2_12arc.grd -I12c
#grdsample EastUS_DEM/3_6arc.grd -GEastUS_DEM/3_12arc.grd -I12c
#grdsample EastUS_DEM/4_6arc.grd -GEastUS_DEM/4_12arc.grd -I12c

#Make the intensity file for shading
#grdgradient EastUS_DEM/1_12arc.grd -A270 -GEastUS_DEM/1_12arc_int.grd -V -Nt1/1000/0
#grdgradient EastUS_DEM/2_12arc.grd -A270 -GEastUS_DEM/2_12arc_int.grd -V -Nt1/1000/0
#grdgradient EastUS_DEM/3_12arc.grd -A270 -GEastUS_DEM/3_12arc_int.grd -V -Nt1/1000/0
#grdgradient EastUS_DEM/4_12arc.grd -A270 -GEastUS_DEM/4_12arc_int.grd -V -Nt1/1000/0

#Downsample the grid files
#grdsample EastUS_DEM/r1_6arc.grd -GEastUS_DEM/r1_12arc.grd -I12c
#grdsample EastUS_DEM/r2_1_6arc.grd -GEastUS_DEM/r2_1_12arc.grd -I12c
#grdsample EastUS_DEM/r2_2_6arc.grd -GEastUS_DEM/r2_2_12arc.grd -I12c
#grdsample EastUS_DEM/r3_6arc.grd -GEastUS_DEM/r3_12arc.grd -I12c
#grdsample EastUS_DEM/r4_6arc.grd -GEastUS_DEM/r4_12arc.grd -I12c

#Make the intensity file for shading
#grdgradient EastUS_DEM/r1_12arc.grd -A270 -GEastUS_DEM/r1_12arc_int.grd -V -Nt1/1000/0
#grdgradient EastUS_DEM/r2_1_12arc.grd -A270 -GEastUS_DEM/r2_1_12arc_int.grd -V -Nt1/1000/0
#grdgradient EastUS_DEM/r2_2_12arc.grd -A270 -GEastUS_DEM/r2_2_12arc_int.grd -V -Nt1/1000/0
#grdgradient EastUS_DEM/r3_12arc.grd -A270 -GEastUS_DEM/r3_12arc_int.grd -V -Nt1/1000/0
#grdgradient EastUS_DEM/r4_12arc.grd -A270 -GEastUS_DEM/r4_12arc_int.grd -V -Nt1/1000/0


#NEW MADRID DEM DOWNSAMPLING...
#Downsample the grid files
#grdsample NewMadrid_DEM/1_3arc.grd -GNewMadrid_DEM/1_6arc.grd -I6c
#grdsample NewMadrid_DEM/2_3arc.grd -GNewMadrid_DEM/2_6arc.grd -I6c
#grdsample NewMadrid_DEM/3_3arc.grd -GNewMadrid_DEM/3_6arc.grd -I6c
#grdsample NewMadrid_DEM/4_3arc.grd -GNewMadrid_DEM/4_6arc.grd -I6c

#Make the intensity file for shading
#grdgradient NewMadrid_DEM/1_6arc.grd -A270 -GNewMadrid_DEM/1_6arc_int.grd -V -Nt1/1000/0
#grdgradient NewMadrid_DEM/2_6arc.grd -A270 -GNewMadrid_DEM/2_6arc_int.grd -V -Nt1/1000/0
#grdgradient NewMadrid_DEM/3_6arc.grd -A270 -GNewMadrid_DEM/3_6arc_int.grd -V -Nt1/1000/0
#grdgradient NewMadrid_DEM/4_6arc.grd -A270 -GNewMadrid_DEM/4_6arc_int.grd -V -Nt1/1000/0


#CHARLESTON DEM DOWNSAMPLING...
#Downsample the grid files
#grdsample Charleston_DEM/S1_3arc.grd -GCharleston_DEM/S1_6arc.grd -I6c
#grdsample Charleston_DEM/S2_3arc.grd -GCharleston_DEM/S2_6arc.grd -I6c

#Make the intensity file for shading
#grdgradient Charleston_DEM/S1_6arc.grd -A270 -GCharleston_DEM/S1_6arc_int.grd -V -Nt1/1000/0
#grdgradient Charleston_DEM/S2_6arc.grd -A270 -GCharleston_DEM/S2_6arc_int.grd -V -Nt1/1000/0


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

#print a message
echo -e "Plotting East US Map..."

if [ $plotDEM -eq 1 ]
then
	#plot topography (SRTM Tiles from getSRTM.pl)
	echo -e "plotting EastUS DEM tile 1/4..."
	grdimage EastUS_DEM/1_12arc.grd -X$xShift -Y$yShift -R$area -JM$scale -CColorTopo.cpt -IEastUS_DEM/1_12arc_int.grd -P -K >> $file 
	echo -e "plotting EastUS DEM tile 2/4..."
	grdimage EastUS_DEM/2_12arc.grd -R -JM -CColorTopo.cpt -IEastUS_DEM/2_12arc_int.grd -O -K >> $file 
	echo -e "plotting EastUS DEM tile 3/4..."
	grdimage EastUS_DEM/3_12arc.grd -R -JM -CColorTopo.cpt -IEastUS_DEM/3_12arc_int.grd -O -K >> $file 
	echo -e "plotting EastUS DEM tile 4/4..."
	grdimage EastUS_DEM/4_12arc.grd -R -JM -CColorTopo.cpt -IEastUS_DEM/4_12arc_int.grd -O -K >> $file 
	echo -e "";
	
	echo -e "filling EastUS DEM hole 1/5 (Iowa)..."
	grdimage EastUS_DEM/r1_12arc.grd -R -JM -CColorTopo.cpt -IEastUS_DEM/r1_12arc_int.grd -O -K >> $file 
	echo -e "filling EastUS DEM hole 2/5 (Florida/Georgia)..."
	grdimage EastUS_DEM/r2_1_12arc.grd -R -JM -CColorTopo.cpt -IEastUS_DEM/r2_1_12arc_int.grd -O -K >> $file 
	echo -e "filling EastUS DEM hole 3/5 (Florida/Georgia)..."
	grdimage EastUS_DEM/r2_2_12arc.grd -R -JM -CColorTopo.cpt -IEastUS_DEM/r2_2_12arc_int.grd -O -K >> $file
	echo -e "filling EastUS DEM hole 4/5 (Carolinas)..."
	grdimage EastUS_DEM/r3_12arc.grd -R -JM -CColorTopo.cpt -IEastUS_DEM/r3_12arc_int.grd -O -K >> $file
	echo -e "filling EastUS DEM hole 5/5 (Rhode Island/Massachusetts)..."
	grdimage EastUS_DEM/r4_12arc.grd -R -JM -CColorTopo.cpt -IEastUS_DEM/r4_12arc_int.grd -O -K >> $file
	echo -e "";
	
	#plot basemap
	psbasemap -R -JM -Bf5.0a5.0/f5.0a5.0:horizontal:WSen -O -K >> $file
		
else
	#plot basemap
	psbasemap -X$xShift -Y$yShift -R$area -JM$scale -Bf5.0a5.0/f5.0a5.0:horizontal:WSen -P -K >> $file

fi

#plot coastline
pscoast -R -JM -S$waterColor -Ir/1p,$waterColor -W0.5p -Df -N1/2p -N2/2p -O -K >> $file

#plot earthquake data WITHOUT magnitude scaling
#psxy NEIC_PDE_1973-2011-8-22.dat -R -JM -Sc0.2i -G$redColor -W1.25p -O -K >> $file
psxy ANSS_1972-2012.mag -R -JM -Sc0.16i -G$redColor -W1.25p -O -K >> $file
#psxy NEIC_SignificantUS_EQs.dat -R -JM -Sc0.2i -G25/10/175 -W1.25p -O -K >> $file
                 
#plot a star at the locations of significant Historical Events
psxy -R -JM -Sa0.5i -Ggold -W1.25p -O -K <<END>> $file
-70.300 42.700
-90.400 35.600
-89.600 36.300
-89.600 36.500
-80.000 32.900
-89.400 37.000
-80.700 37.300
-88.373 37.911
-77.9330 37.9360
END

#plot a text label for significant historical EQ's
pstext -R -JM -D-0.25i/-0.1i -Ggold -S2p -O -K <<END>> $file
-70.300 42.700 24 0 1 TR 1755 VIII Cape Ann
-90.400 35.600 24 0 1 TR 1811 M7.7 and M7.0 New Madrid
-89.600 36.300 24 0 1 TR 1812 M7.5 New Madrid
-89.600 36.500 24 0 1 BR 1812 M7.7 New Madrid
-80.000 32.900 24 0 1 TR 1886 M7.3 Charleston
-89.400 37.000 24 0 1 BR 1895 M6.6 Charleston
-80.700 37.300 24 0 1 TR 1897 M5.6 Giles County
-88.373 37.911 24 0 1 TR 1968 M5.4 Southern Illinois
-77.9330 37.9360 24 0 1 TR 2011 M5.8 Virginia
END

#plot the New Madrid rectangle to show where the subregion is
psxy -R -JM -W4p,gold -O -K <<END>> $file
-91 35
-91 38
-88 38
-88 35
-91 35
END

#plot the Charleston rectangle to show where the subregion is
psxy -R -JM -W4p,gold -O -K <<END>> $file
-81.5 32
-81.5 34
-78   34
-78   32
-81.5 32
END

#change axis font size for the legend
gmtset ANNOT_FONT_SIZE=20p

#plot a legend
pslegend -R -JM -Dx0.4i/0.4i/12.0i/6.7i/BL -F -Gwhite -O -K  <<END>> $file
G 0.25i
H 30 1 Eastern United States Earthquakes > M@-w@-0  (1972-2012)
G 0.15i
H 24 1 13,360 Total Events in 40 Years; 334 Events/Average Year
G 0.25i
H 20 1 Compiled by Julia T. Irizarry and Scott T. Marshall
G 0.14i
H 20 1 Department of Geology, Appalachian State University
G 0.14i
H 20 1 Earthquake Epicenters from the ANSS Catalog 
G 0.4i
M 275 37 200k+u f
G 0.25i
M 275 37 200m+u f 
G 0.25i
B ColorTopo.cpt 0.5i 0.5i -B:"Elevation Above Sea Level (meters)": --FRAME_PEN=2p --LABEL_FONT_SIZE=20 --LABEL_FONT=0
END


#------------------------------------------------------------------------------------------------#
#------------------------------------------------------------------------------------------------#
#----------------------------   NEW MADRID MAP INSET   ------------------------------------------#
#------------------------------------------------------------------------------------------------#
#------------------------------------------------------------------------------------------------#

#print a message
echo -e "Plotting New Madrid Inset..."

#define New Madrid Area and scale
NMarea=-91/-88/35/38
NMscale=11.5i
#Set the new XYshifts to move the insets over
xShift=27.25i
yShift=1.5i
#change frame width
gmtset FRAME_WIDTH=3p
#change axis font size
gmtset ANNOT_FONT_SIZE=30p
#change to simpl e basemap
gmtset BASEMAP_TYPE=plain

if [ $plotDEM -eq 1 ]
then
	#plot New Madrid topography
	echo -e "plotting New Madrid DEM tile 1/4..."
	grdimage NewMadrid_DEM/1_6arc.grd -X$xShift -Y$yShift -R$NMarea -JM$NMscale -CColorTopo.cpt -INewMadrid_DEM/1_6arc_int.grd -O -K >> $file 
	echo -e "plotting New Madrid DEM tile 2/4..."
	grdimage NewMadrid_DEM/2_6arc.grd -R -JM -CColorTopo.cpt -INewMadrid_DEM/2_6arc_int.grd -O -K >> $file 
	echo -e "plotting New Madrid DEM tile 3/4..."
	grdimage NewMadrid_DEM/3_6arc.grd -R -JM -CColorTopo.cpt -INewMadrid_DEM/3_6arc_int.grd -O -K >> $file 
	echo -e "plotting New Madrid DEM tile 4/4..."
	grdimage NewMadrid_DEM/4_6arc.grd -R -JM -CColorTopo.cpt -INewMadrid_DEM/4_6arc_int.grd -O -K >> $file
	echo -e "";
	
	#plot New Madrid basemap in SE corner of map
	psbasemap -R -JM -B0.5a1.0/0.5a1.0:horizontal:WSen -O -K >> $file
		
else
	#plot New Madrid basemap in SE corner of map
	psbasemap -X$xShift -Y$yShift -R$NMarea -JM$NMscale -B0.5a1.0/0.5a1.0:horizontal:WSen -O -K >> $file
	
fi

#plot coastline 
pscoast -R -JM -S$waterColor2 -Df -A1 -N1/8 -N2/8 -Ia/4,$waterColor2 -O -K >> $file

#plot earthquake data withOUT magnitude scale
#psxy NEIC_PDE_1973-2011-8-22.dat -R -JM -Sc0.13i -G$redColor -W1.25p -O -K >> $file
psxy ANSS_1972-2012.mag -R -JM -Sc0.08i -G$redColor -W1.25p -O -K >> $file

#plot a star at the locations of significant Historical Events
psxy -R -JM -Sa0.5i -Ggold -W1.25p -O -K <<END>> $file
-90.400 35.600
-89.600 36.300
-89.600 36.500
-89.400 37.000
-88.373 37.911
END

#plot a text label for significant historical EQ's
pstext -R -JM -D-0.25i/-0.1i -Ggold -S2p -O -K <<END>> $file
-89.600 36.300 24 0 1 TR 1812 M7.5 New Madrid
-89.600 36.500 24 0 1 BR 1812 M7.7 New Madrid
-89.400 37.000 24 0 1 BR 1895 M6.6 Charleston
-88.373 37.911 24 0 1 TR 1968 M5.4 Southern Illinois
END

#plot a text label for significant historical EQ's
pstext -R -JM -D0.25i/0.1i -Ggold -S2p -O -K <<END>> $file
-90.400 35.600 24 0 1 BL 1811 M7.7 and M7.0 New Madrid
END

#change map scale height               
gmtset MAP_SCALE_HEIGHT=0.15i
#change tick length 
gmtset TICK_LENGTH=0.1i
#change distance from scale to distance measurements
gmtset ANNOT_OFFSET_PRIMARY=0.1i
#change axis font size
gmtset ANNOT_FONT_SIZE=16p

#plot NM legend
pslegend -R -JM -Dx11.0i/0.5i/5.0i/2.25i/BR -F -Gwhite -O -K --FRAME_PEN=2p <<END>> $file
G 0.15i
H 22 1 New Madrid Seismic Zone
G 0.35i
M 270 36 50k+u f  
G 0.25i
M 270 36 50m+u f 
END


#------------------------------------------------------------------------------------------------#
#------------------------------------------------------------------------------------------------#
#----------------------------   CHARLESTON MAP INSET   ------------------------------------------#
#------------------------------------------------------------------------------------------------#
#------------------------------------------------------------------------------------------------#

#print a message
echo -e "Plotting Charleston Inset..."

#define Charleston Area and scale
CHarea=-81.5/-78/32/34
CHscale=10i
#Set the new XYshifts to move the insets over
xShift=1.5i
yShift=16.0i
#change axis font size
gmtset ANNOT_FONT_SIZE=30p

if [ $plotDEM -eq 1 ]
then
	#plot Charleston topography
	echo -e "plotting Charleston DEM tile 1/2..."
	grdimage Charleston_DEM/S1_6arc.grd -X$xShift -Y$yShift -R$CHarea -JM$CHscale -CColorTopo.cpt -ICharleston_DEM/S1_6arc_int.grd -O -K >> $file 
	echo -e "plotting Charleston DEM tile 2/2..."
	grdimage Charleston_DEM/S2_6arc.grd -R -JM -CColorTopo.cpt -ICharleston_DEM/S2_6arc_int.grd -O -K >> $file 
	echo -e "";
	
	#plot Charleston basemap-- in SE corner of map
	psbasemap -R -JM -Bf0.5a1.0/f0.5a1.0:horizontal:WSen -O -K >> $file
	
else
	#plot Charleston basemap-- in SE corner of map
	psbasemap -X$xShift -Y$yShift -R$CHarea -JM$CHscale  -Bf0.5a1.0/f0.5a1.0:horizontal:WSen -O -K >> $file
	 
fi

#plot coastline
pscoast -R -JM -S$waterColor2 -W0.5p -Df -A1 -N1/8 -N2/8 -Ia/4,$waterColor2 -O -K >> $file

#plot earthquake data withOUT magnitude scale
#psxy NEIC_PDE_1973-2011-8-22.dat -R -JM -Sc0.1i -G$redColor -W1.25p -O -K >> $file
psxy ANSS_1972-2012.mag -R -JM -Sc0.08i -G$redColor -W1.25p -O -K >> $file

#plot a star at the locations of significant Historical Events
psxy -R -JM -Sa0.5i -Ggold -W1.25p -O -K <<END>> $file
-80.000 32.900
END

#plot a text label for significant historical EQ's
pstext -R -JM -D0.25i/0.1i -Ggold -S2p -O -K <<END>> $file
-80.000 32.900 24 0 1 BL 1886 M7.3 Charleston
END

#plot a star at the location of downtown Charleston, SC
psxy -R -JM -Sc0.2i -Gwhite -W1.25p -O -K <<END>> $file
-79.931065 32.776473
END

#plot a text label for Charleston
pstext -R -JM -D-0.15i/-0.1i -Gwhite -S2p -O -K <<END>> $file
-79.931065 32.776473 18 0 1 TR Charleston, SC
END

#change axis font size
gmtset ANNOT_FONT_SIZE=16p
 
#plot a legend
pslegend -R -JM -Dx9.5i/0.5i/4.0i/2.25i/BR -F -Gwhite -O --FRAME_PEN=2p <<END>> $file
G 0.15i
H 18 1 Middleton Place-Summerville
H 18 1 Seismic Zone
G 0.35i
M 280 33.5 50k+u f
G 0.25i
M 280 33.5 50m+u f
END

#print the final line to the .eps file so that it will open properly
#echo showpage >> $file

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

#display the map
#gv $file &
illustrator $file &

#convert to pdf
echo -e "Converting to pdf...";
ps2raster $file -Tf
#Now open the pdf with acrobat
acrobat EastUS.pdf &
#convert to jpg
echo -e "Converting to jpg...";
#ps2raster $file -Tj
#convert to png
echo -e "Converting to png...";
#ps2raster $file -TG

echo -e "Done! \a"
exit;

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